1 Introduction
Snow is a first-order control on ecosystems, influencing energy balance, hydrology, geomorphic processes, flora, and fauna (
Billings and Bliss, 1959;
Totland and Alatalo, 2002;
Homan et al., 2011;
Tahir et al., 2011;
Moore et al., 2015). The US mountain west has seen an earlier arrival of spring weather since the 1970s, leading to snow melting 14‒45 days earlier than average (
Stewart, 2009), changes in plant phenology (
Totland and Alatalo, 2002;
Steltzer et al., 2009;
Cornelius et al., 2013), altered water resource availability (
Sproles et al., 2013), and a larger area burned by wildfires (
Westerling et al., 2006;
O’Leary III et al., 2016). Furthermore, model-based projections suggest that snowmelt may arrive an average of 10‒40 days earlier across the continental US by the year 2100 (
Stewart, 2009). The snow-dominated watersheds of the Cascade Mountains may be particularly vulnerable to climate change, as the mid-elevations of these mountains spend much of the winter near 0°C, and a change of only a few degrees could therefore alter the snowpack substantially (
Nolin and Daly, 2006;
Sproles et al., 2013) as an increase in freezing elevation level and possibility of rain-on-snow events deplete snowpack rapidly (
McCabe et al., 2007). Perhaps as a harbinger of change, the years 2014 and especially 2015, experienced an unusually low snow volume and early snowmelt in the Cascade Mountains (
Margulis et al., 2016a;
Sproles et al., 2017). Events like these and the widely-anticipated future reduction of snowpack and earlier arrival of spring snowmelt have created an increased interest in monitoring the timing of snowmelt.
While snowmelt timing is often defined as a single value for a watershed using streamflow (
USGS, 2018) or measured in a single location at meteorological stations (
NRCS, 2016), snowmelt timing is in reality a highly spatially heterogeneous process due to fine-scale variability in snow accumulation and spring heating patterns. Recent research has therefore started to explore the use of remotely sensed data to map snowmelt timing across the landscape as decreasing snow-covered area (
Hall et al., 2012,
2015;
O’Leary III et al., 2016;
Nolin et al., 2017). Snowmelt timing has profound ecological implications at the landscape level (
Billings and Bliss, 1959;
Westerling et al., 2006;
Inouye, 2008) and many snowmelt timing investigations have been motivated by questions of snow’s control on ecological processes with the goal of providing a Day Of Year (DOY) for snowmelt timing to compare with similar DOY products for plant phenology (
Jönsson and Eklundh, 2004;
Zhang et al., 2006;
Ganguly et al., 2010;
Tahir et al., 2011). Previous attempts to define the timing of snow cover loss using optical remote sensing imagery have been largely unsuccessful due to cloud interference and validation issues (
Narasimhan and Stow, 2010). However, studies from the wildfire literature have successfully identified the influence of snowmelt timing on annual area burned using cloud-free datasets (
Semmens and Ramage, 2012;
O’Leary III et al., 2016). With an increasing interest on climate’s influence on plant phenology (
University of Wyoming, 2017;
Sherwood et al., 2017) and wildfire (
Westerling, 2016), there is an obvious need for a gridded snowmelt timing dataset for use in ecological investigations (
Littell et al., 2009;
Abatzoglou and Kolden, 2013;
O’Leary III et al., 2018).
In this paper we describe the development of annual Snowmelt Timing Maps (STMs) (
O’Leary III et al., 2017), which identify the DOY that a given location becomes snow-free for the years 2001‒2015, based on the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra satellite standard snow-cover product MOD10A2 (
Hall and Riggs, 2016;
Riggs et al., 2017). We address the question: what was the spatial extent and temporal magnitude of the snowmelt timing anomaly of 2015 in the Cascade Mountains? In answering that question, we identify the areas having the greatest snowmelt timing anomaly in terms of their elevation, relative topography, and geographic location. We focus on the Cascade Mountains, and report specifically on results within three National Park Service (NPS) unit focus areas. Finally, we discuss the successes and limitations of this approach, with a particular focus on applications within the fields of hydrology, ecology, and land management.
2 Methods
The Cascade Mountains of the western conterminous US are part of a volcanic cordillera stretching from central Washington, through Oregon, and into northern California (Fig. 1). This region is known for its deep seasonal snowpack, expansive protected and recreational lands, and towering glacier-flanked volcanoes (
USGS, 2014;
NPS, 2016a). Here, we consider the Cascade Mountains along with case studies that focus on three National Park Service units: Mt Rainier (MORA), Crater Lake (CRLA), and Lassen Volcanic (LAVO) National Parks. These protected areas are the sites of active and long running ecological investigations and monitoring programs sponsored by federal agencies, non-profits, and academic institutions (
NPS, 2016a,
b,
c). Additionally, these high-elevation areas are particularly vulnerable to climate change (
Rangwala and Miller, 2012), and many at-risk species compete with invasive species for resources and survival (
Wolf et al., 2016). Finally, these three study areas contain distinct environments offering a variety of snowmelt patterns-- the massive glaciers of Mt Rainier contrasting with the subalpine meadows of Crater Lake, and the drier California climate of Mt Lassen volcano.
We developed the STMs (
O’Leary III et al., 2017) spanning the area of North America, Greenland, and Iceland using a time-series analysis of the MOD10A2 Collection 5 data product (
Hall et al., 2006). The MOD10A2 product is derived from a fractional snow cover product with the snow/no-snow threshold defined by 50% snow covered area. In this composite product, eight consecutive days of observations are compiled into a single image representing the ‘maximum’ snow cover for that period, with each pixel defined as snow-covered if snow was observed for one or more days during the composite period.
To define snowmelt timing, we recorded the DOY for the first snow-free image following snow cover, or a mixture of clouds and snow cover. To account for cloud-obscured MOD10A2 images we interpolated the dates between the last observed snow cover and the first observed snow-free date (Table 1). We recorded cloud interference as the number of consecutive 8-day composite images that were cloud-covered+1. Pixels with more than four consecutive cloud-obscured images (32 days) between snow and no snow were omitted because of their poor validation results. We also calculated the count of STM values (‘Count’) identifying the number of years on record where the STM has a value (i.e., snow was present at some point that year, and snowmelt timing was successfully calculated). While late-season snowstorms contribute appreciable water to the ecosystem, they are not the seasonal indicator and hydrologic input of the primary snowpack departure that we sought to capture in the STMs. To remove late-season snowstorms we disregarded any snow readings that occur following a 48-day period without snow (six consecutive 8-day composite images) beginning on February 26 or later. To validate our STMs, we compared the snowmelt timing values with the SNOw TELemetry (SNOTEL) station network spanning the western United States (
NRCS, 2016) (
n=854). For a detailed explanation of STM development and validation see Appendix 1.
After validating against SNOTEL data, we divided our study area into four sections: the Cascade Mountains as a whole, and MORA, CRLA, and LAVO National Parks. We calculated Snow-Cover Depletion Curves (SCDCs) (
Hall et al., 2015) for each region in software R version 3.2.3 by calculating the percentage of pixels (
Y-axis) in each annual STM image with snowmelt DOY value greater than or equal to the daily time step (
X-axis). The snowpack of the Cascade Mountains experiences substantially different dynamics at low, medium, and high elevations (
Gleason et al., 2017), though comparing a particular elevation between the northern and southern extents of the Cascade Mountains may find stark differences (e.g., at an elevation of 2000 m in MORA, one may encounter glaciers, while the same elevation may be hot, bare, and dry in LAVO). To capture these differences, we separated each study area into three elevation quantiles of equal area relative to the elevation range of each park (Table 2). We then found the mean snowmelt DOY for 2001‒2015 and calculated the 2015 snowmelt anomaly for each area and elevation quantile.
3 Results
Validation results of snowmelt timing values between the STMs and SNOTEL stations reveals a high level of agreement. Differences between the STMs and all available SNOTEL stations have a mean error of 4.31 days. When we include only those stations within the Cascade Mountains we find a mean error of –4.76 days. These errors are less than the 8-day temporal resolution of the MOD10A2 product and therefore indicate that our STM algorithm captures on-the-ground snowmelt timing accurately. For a detailed description of this validation see Appendix 1.
The STMs show snowmelt timing gradually changing across the landscape, with low elevation areas and the arid lands to the east melting earlier than the land west of and including the crest of the Cascade Mountains (Fig. 2). When we investigate the 2015 snowmelt anomaly (relative to the 2001‒2015 average), we find widespread earlier melt across the study area, with the vast majority of pixels experiencing a negative anomaly (Fig. 3). The average snowmelt timing anomaly was ‒41.45 days for the Cascade Mountains as a whole, and ‒32.94, ‒40.12, and ‒47.49 days for MORA, CRLA, and LAVO, respectively (Fig. 4, Table 3). Looking at the relative elevation of the Cascade Mountains as a whole, the lowest elevation areas experience less of an anomaly (‒25.40 days) than do the middle (‒45.05 days) and highest (‒50.46 days) areas. This pattern is reversed for each national park, where the highest elevations experienced a lesser anomaly than the lower elevations within each park (Table 3). The year 2015 shows a generalized early melt, with much of the Cascade Mountains melting weeks to months earlier than the 15-year mean (Figs. 2 and 3). For the Cascade Mountains as a whole, and each national park studied here, there were locations where snowmelt timing occurred more than 100 days earlier than the 2001‒2015 average (Fig. 3). 93.3% of the land area of the Cascade Mountains experienced an earlier than average snowmelt timing in 2015, along with 95.1% of MORA, 98.6% of CRLA, and 99.8% of LAVO. Ten percent of the land area of the Cascade Mountains had a snowmelt timing anomaly of ≥83.73 days earlier than average, and LAVO similarly had ten percent of its land area melt ≥80.53 or more days earlier than average. However, MORA and CRLA were somewhat insulated from the greatest anomaly magnitudes, with ten percent of their land area melting>57 days earlier than average (Table 4).
The cloud interference output shows that the Cascade Mountains are subject to cloud interference during the snowmelt period (Fig. 5). In particular, cloud interference occurs most in the northern reaches, the low-elevation areas of the Columbia River gorge, and the western flanks of the Cascade Mountains crest. When we consider the number of years on record with a STM value, it is evident that the areas with the lowest Count of STM values often coincide with the highest cloud interference values (Fig. 5). While the mountainous areas of western Washington, Oregon, and northwestern California show high Count values approaching 15 years (i.e., they experience snow accumulation and snowmelt every year on record), the lower-elevation valleys and some parts of the central Oregon and the northeastern California altiplano show lower Count values, indicating an inconsistent inter-annual snowpack. Notably, MORA experiences high cloud interference values along with high Count values, suggesting that the snowmelt season is particularly cloudy at this location, leading to ambiguous detection of snowmelt timing. SCDCs for the Cascade Mountains and national parks show that in 2015 snow melts earlier than all other years throughout the melt season (Fig. 6). For each region 2015 has an earlier initial melt, earlier complete melt, and a lower SCA than all other years for the entire melt season (with the exception of LAVO melting out completely in 2001 several days earlier than in 2015).
4 Discussion
4.1 STMs and the 2015 snowmelt timing anomaly
The STM algorithm was successful in creating intuitive, ecologically sensible, and easy-to-employ maps for environmental analysis, which compare favourably with in-situ snow measurements. For large regions, visual inspection of the STMs is straightforward (Fig. 2), showing landscape patterns that relate to topography and regional climate. Viewing smaller areas (such as the national parks of interest, Fig. 4) we clearly see pixilation from the 500 m resolution of the MODIS snow-cover product. Still, subtle differences between elevations, slope gradients, and aspects highlight the landscape-level spatial heterogeneity of snowmelt and may yield insight for land managers or other stakeholders.
The year 2015 had a widespread early snowmelt event across the Cascade Mountains. Our SCDCs show that 2015 had substantially lower snow-covered area than all other years on record for the Cascade Mountains as a whole and each national park of interest (Fig. 6). These results are consistent with other analyses of the 2015 snow year throughout the Cascade and Sierra Nevada ranges (
Margulis et al., 2016a,
b;
Sproles et al., 2017). Compared with the spatial mean of snowmelt DOY values, we see that the 2015 early snowmelt event was widespread, though non-uniform, throughout the study area. Stratifying study areas by elevation we see that for the Cascade Mountains as a whole, the high elevations had the greatest anomaly (Table 3). However, the individual parks are much higher in elevation, with LAVO existing completely within the Cascade Mountains highest quantile. When we look at the quantiles within these smaller study areas we see that the lowest elevations experienced the greatest anomalies (LAVO Q1 is particularly devastated at ‒62.35 days), while the highest elevations experienced the smallest anomalies within the parks (~‒30 days) (Fig. 3, Table 3). Some high-elevation mountain ecosystems experience greater warming than the global average (
Rangwala and Miller, 2012), and while we see this effect when looking at the Cascade Mountains as a whole, our results show that the highest elevation areas within the national parks experienced less of a snowmelt anomaly than the lower-elevation slopes surrounding the volcanoes (Table 3). This may be driven by many different factors, including land cover, stochastic weather events, and vegetation. Forest composition influences snowpack dynamics by changing canopy interception and radiation balance, and forested and non-forested areas retain and ablate snow in different ways (
Roth and Nolin, 2017), though it is unclear how these forces have contributed to this event of 2015. Perhaps the areas near treeline experience the greatest anomaly because they are typically slightly below freezing, needing an increase of only a few °C to cause major snowmelt timing changes, while the highest elevations were still well below freezing for much of the melt season (
Mote et al., 2005). While there is no single driver of the elevational pattern that we see in our results, further investigations may consider these elements when isolating snowmelt events in particular locations.
Cloud cover occasionally prevents an ideal snowmelt detection, leading to a temporally interpolated snowmelt DOY value (Table 1, Case 3+), and a reduced temporal precision. In our validation, increasing cloud cover corresponded with greater errors between the STMs and SNOTEL stations, with the STMs estimating a later date of snow loss than the SNOTEL stations (Appendix 1). Generally, cloud interference is greater in the northern end of the Cascade Mountains, as Washington State, and MORA in particular, show higher cloud interference values than the Oregon and California parts of the ecoregion (Fig. 5). Throughout the range, cloud cover is also greater along the western foothills, owing to the orographic cloud formation and precipitation driven by the uplift of the Cascade Mountains. Predominant weather patterns direct humid coastal air from west to east, which frequently forms clouds on the western side of the mountains as topography-driven uplift causes rapid cooling and precipitation. This also creates a well-known “rain shadow” east of the Cascade Mountains where sunny, dry days are common. Ultimately, the patterns shown in Fig. 5 are sensible and expected for this ecoregion.
There are many areas that experience snow cover for some, but not all, years on record (Fig. 5). Many of the low-elevation valleys on the western flanks of the Cascade Mountains, as well as the low-elevation areas separating LAVO from the rest of the Cascade Mountains, have a Count value less than 15. We calculated the mean cloud interference values by ignoring years without a snowmelt value, so within areas that have low Count values a single year with a high cloud interference value may result in a higher mean cloud interference value than an area with 15 years of coverage. Keeping this effect in mind, we see the areas within Fig. 5 with a low Count often have a high mean cloud interference value, particularly within the Puget Sound and Willamette Valley. Interestingly, within the Columbia River Valley we see a high Count value with a high cloud interference value, likely owing to persistent fog and low clouds that funnel into the Columbia River Valley from the Pacific coast during the melt season. For our visualizations we mask pixels with Count values less than 8 to highlight regions with a consistent perceptible snowmelt (>50% occurrence).
4.2 Limitations
There are numerous limitations involved in the application of our STMs. Of primary concern is the 8-day resolution of this product. While the MOD10A2 does improve upon the MOD10A1 daily product in terms of cloud-free data, it is still limited by some cloud cover, its 8-day resolution, and inherent snow reporting biases (
Hall and Riggs, 2007). This 8-day resolution differs from the daily or hourly resolution available using meteorological stations (SNOTEL) and limits the precision of any snowmelt observation to eight days. Still, when investigating large areas, particularly if averaging snowmelt DOY values, this fairly coarse temporal resolution is smoothed by the many points available from the moderate spatial resolution, and the 8-day composites streamline data processing compared to daily resolution products. Alternative snow datasets are available, though they have their own limitations. The SNOw Data Assimilation System (SNODAS) offers 1 km snow products including extent, SWE, temperature, and more, from late 2003‒present (
Barrett, 2003). This product is limited by its methods for this kind of application (computer modelled from several
in-situ and remotely-sensed inputs), and therefore likely has a reduced accuracy at the margins of snow cover compared to the MOD10 products. The widely-employed Ice Mapping System (IMS) (
Ramsay, 1998) is also derived from multiple inputs, though spatial resolution for the time period of this study is poor at 24 km, for the earlier part of our study period, with the 4 km and 1 km resolution IMS products becoming available more recently. Further research should employ the cloud-gap-filled daily MOD10A1F product in Version 6 (
Hall and Riggs, 2016), though this would require a non-trivial revision of the Python scripts used for the snowmelt timing algorithm, and substantially increase storage and computational demand. These methods will also be robust in the future using the Visible Infrared Imaging Radiometer Suite (VIIRS) (
Riggs et al., 2017).
When looking closely at the national parks of interest we see areas with low Count values around the glaciers of MORA, and the seasonal lake ice of CRLA. In MORA, some years show extensive snowmelt, and possibly glacier retreat, exposing crags and rocky slabs that are typically snow and ice covered throughout the year. This leads to snowmelt values in locations that do not usually experience snowmelt, resulting in a low Count values on and around this volcano’s glaciers. These areas also experience extensive cloud cover for much of the year, also leading to low Count values. The possible signal from glacier retreat should not be confounded with high cloud interference, and deserves further investigation, though we recommend the use of higher spatial resolution products for glacier monitoring in temperate regions. CRLA experiences different lake ice extents depending on the winter. When this lake ice is snow-covered it is classified as snow in the MOD10A2 data product, leading to snowmelt timing values where there typically are none. These effects lead to subtle changes in the SCA by increasing the total snow-covered area of the study site (Fig. 6), and may be of interest for particular applications, though for our purposes these data points have little influence on the bigger picture in terms of ecology and hydrology.
Validating the STMs with SNOTEL measurements is useful, though the differences between the data collections methods complicate their comparison. Cloud interference is not a problem for SNOTEL locations, as SWE observations are collected using a pressure transducer. This difference may lead to snowmelt timing disagreement where the STM for that year is cloud-interfered. SNOTEL locations are reported by latitude and longitude to the hundredth of a decimal degree, and this low spatial precision may lead to registration mismatch between SNOTEL location and the STM pixel that it exists within. Additionally, the forest (or lack thereof) surrounding any SNOTEL station may influence this reading as different forest types have substantially different interception and radiation properties (
Lundquist et al., 2013;
Dickerson-Lange et al., 2017;
Roth and Nolin, 2017). Finally, the 30 cm SWE threshold may not be appropriate for all ecosystems and years, and adjusting this value may improve the comparison results (Appendix 1).
4.3 Further investigation
Continued development of the snowmelt detection algorithm will yield further insight into quantifying snowmelt timing. The STMs are currently limited to an 8-day resolution stemming from the MOD10A2 parent data. Creating new STMs using the Collection 006 cloud-gap-filled product, MOD10A1F, could improve temporal resolution, possibly improving the SNOTEL comparison and delivering a higher temporal precision. Improving explicit reporting of details such as cloud coverage and late-season snow events may benefit certain users who demand specific criteria for ecological modelling. Expanding the spatial domain of this dataset into a global data product would be a substantial improvement but would require experts from around the world for validation and harmonization with local snow cycles. Research into the specific meteorological conditions leading up to and including the snowmelt season of 2015, including air temperature and Pacific Ocean surface temperature anomalies, could help to pair some of the anomaly characteristics with particular weather events such as shifts in freezing level or precipitation. This would be particularly insightful to explore why the low elevations of the national parks experienced the greatest anomalies. Finally, collaboration with researchers in many fields will yield insights into improvements that would benefit further iterations of the STMs. These STMs have been found to have significant relationships with plant phenology data based on the normalized-difference vegetation index, within CRLA (
O’Leary III et al., 2018) and similar snowmelt timing maps are known to have ecologically significant relationships with wildfire (
Semmens and Ramage, 2012;
O’Leary III et al., 2016), demonstrating that many other ecological questions can be addressed using the STMs.
5 Conclusions
We developed an algorithm for time series analysis of the MOD10A2 MODIS snow-cover data product and used this algorithm to quantify snowmelt timing and create snowmelt timing maps (STMs) for North America from 2001‒2015 (
O’Leary III et al., 2017). We validated the STMs against
in-situ SNOTEL stations located throughout the western United States including Alaska, which proved to have strong agreement (Appendix 1). The derived snowmelt timing values follow topography, with elevation, slope, and aspect all contributing to spatial heterogeneity. Comparing snow years from 2001‒2015 we present snow-cover depletion curves (SCDCs) for the Cascade Mountains as a whole, and for Mt Rainier, Crater Lake, and Lassen Volcanic National Parks individually, and identify spatial and temporal snowmelt trends.
Within the 2001–2015 study period, the year 2015 experienced the earliest initial melt, the earliest complete melt, and the lowest snow-covered area throughout the entire melt period for each region of interest, with a single exception (Fig. 6). The 2015 snowmelt averaged 41.45 days earlier snowmelt in the Cascade Mountains as a whole, with the highest elevations having the greatest snowmelt timing anomaly. Interestingly, the national parks in isolation had the largest snowmelt anomaly at their lowest elevations, suggesting that this mountain range is particularly vulnerable to early snowmelt events near and below treeline. Our new method of quantifying snowmelt timing is especially useful to identify locations where seasonal snowpack is highly sensitive to changes in climate, such as the Cascade Mountains (
Nolin and Daly, 2006;
Sproles et al., 2013). These results are relevant in hydrology (
Musselman et al., 2017), cryosciences, earth energy balance, and phenology (
Sherwood et al., 2017 ; O’Leary III et al., 2018), as well as to land managers, and the general public as they keep a close eye on snowpack dynamics of the Cascade Mountains and other lands for the decades to come.
Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature