Changing snow water storage in natural snow reservoirs

This work introduces a novel snow metric, snow water storage (SwS), defined as the integrated area under the snow water equivalent (SWE) curve (units: length-time, e.g., m d). Unlike other widely used snow metrics that capture snow variables at a single point in time (e.g., maximum SWE) or describe temporal snow characteristics (e.g., length of snow season), SwS is applicable at numerous spatial and temporal scales. This flexibility in the SwS metric enables us to characterize the inherent reservoir function of snowpacks and quantify how this function has changed in recent decades. In this research, changes in the SwS metric are evaluated at point, gridded and aggregated scales across the conterminous United States (hereafter US), with a particular focus on 16 mountainous Environmental Protection Agency (EPA) Level III Ecoregions (ER3s). These ER3s account for 72 % of the annual SwS (SwS_A) in the US, despite these ER3s only covering 16 % of the US land area. Since 1982, spatially variable changes in SwS_A have been observed across the US with notable decreasing SwS_A trends in the western US and in the 16 mountainous ER3s. All mountainous ER3 (except for the Northeastern Highlands in New England) exhibit decreasing trends in SwS_A resulting in a 22 % overall decline in SwS_A across mountainous ER3s. The peak monthly SwS (SwS_M) occurs in March at all spatial scales, while the greatest percentage loss of SwS_M occurs early in the snow season, particularly in November. Unsurprisingly, the highest elevations contribute most to SwS_A in all mountain ranges, but the specific elevations that have experienced loss or gain in SwS_A over the 39-year study period vary between mountain ranges. Comparisons of SwS with other snow metrics underscore the utility of SwS, providing insights into the natural reservoir function of snowpacks, irrespective of SWE curve variability or type (e.g., ephemeral, mountain, permanent). As we anticipate a future marked by increased climate variability and greater variability in mountain snowpacks, the spatial and temporal flexibility of snow metrics such as SwS may become increasingly valuable for monitoring and predicting snow water resources.

1 Introduction

Seasonal snow is a critical resource in mountainous regions and at high latitudes across the United States (US), and many other countries, providing an important ecosystem service by functioning as a natural and spatially distributed reservoir (Barnett et al., 2005). These snow reservoirs play a key role in the water cycle by storing water during the cool season and releasing water gradually throughout the warm season when human and ecological demand is the highest (Li et al., 2017). The natural reservoir function of snowpacks is at risk due to anthropogenic climate change, which has been shown to decrease snowpack magnitude and persistence while increasing snowpack variability (Siirila-Woodburn et al., 2021; Scalzitti et al., 2016; Sospedra-Alfonso et al., 2015; Morán-Tejeda et al., 2013). The variability in climatic variables that drive snowpack variability, such as precipitation and temperature, has increased in the recent past and is projected to continue to increase as a result of climate change (Scalzitti et al., 2016; Sospedra-Alfonso et al., 2015; Morán-Tejeda et al., 2013; Pendergrass et al., 2017; Ohmura, 2012). Given the vulnerability of seasonal snow water storage to climate change and the importance of snow-derived water to municipalities, agriculture, ecosystems and hazard forecasters, it is vital to understand how water storage in our natural snow reservoirs is evolving in the context of a changing climate (Immerzeel et al., 2020; Sturm et al., 2017; Barnett et al., 2005; Li et al., 2017; Siirila-Woodburn et al., 2021).

Snow water equivalent (SWE) is a relevant snowpack characteristic for many water resources applications. SWE is the depth of water obtained upon melting a column of snow. Having an estimate of SWE across a watershed is analogous to knowing the stage (water elevation) in a surface reservoir; it quantifies the amount of water being stored for later use. Other directly measurable snow characteristics largely fall into two categories: (1) they can be temporal snapshots that give us information about snow magnitude at a certain point in time or (2) they can provide information about snow timing (Nolin et al., 2021). The 1 April SWE, snow-covered area (SCA) and peak SWE (SWE_max) are examples of temporal snapshot metrics. Snow metrics that give us information about the timing of snow include snow cover duration (SCD), date of snow onset (DSO) and date of snow disappearance (DSD).

Composite snow metrics, such as the snow storage index (SSI) (Hale et al., 2023) and the water tower index (WTI) (Viviroli et al., 2007), are not directly measurable, but they combine information streams in order to relate snowpack to water storage. The SSI indicates the degree to which snowpack delays the timing and magnitude of surface water inputs relative to when it falls as precipitation, whereas the WTI identifies locations where mountain runoff contributes disproportionately to lowland water supplies. Additionally, a global WTI was developed, which ranks all water towers in terms of their water-supplying role and downstream societal and ecological demand (Immerzeel et al., 2020).

A conceptual SWE curve is shown in Fig. 1a. The conceptual SWE curve referenced throughout this paper is for mountain snowpacks and is delineated by three points: the DSO, the peak SWE (SWE_max) and the DSD. SWE accumulation begins at the DSO and continues up to a SWE_max, which may or may not occur on 1 April (Northern Hemisphere). After SWE_max, the ablation phase of the snow season begins and the SWE depth declines until it reaches zero at the DSD. The SCD is captured by the width of the SWE curve. Multiple factors can result in systematic changes in the shape of the SWE curve, including climate change (Lute et al., 2015); natural land cover change, such as wildfire (Gleason et al., 2019) or beetle kill (Pugh and Small, 2012; Boon, 2007; Winkler et al., 2014); and anthropogenic land cover change, such as forest thinning (Krogh et al., 2020; Sun et al., 2022) or logging (Winkler et al., 2005; Troendle and Reuss, 1997).

Figure 1 Conceptual illustration of a SWE curve and various ways that it could change from past to present.

The shaded regions in the other panels of Fig. 1 provide plausible examples of how the SWE curve may have changed from the past to present day. For example, a current SWE curve could be a scaled (reduced) version of a past SWE curve (Fig. 1b). This would result in a later DSO, a lower SWE_max, an earlier DSD and a shorter SCD. Changes in SWE curves could also result from a temporal shift in the historic curve (Fig. 1c). This would not impact SWE_max or SCD, but the metrics including 1 April SWE, DSO and DSD would be affected. Figure 1d gives yet another example of a theoretical current scenario, compared to a historic one. In this case, the shape of the conceptual SWE curve is changed by repeated accumulation and melt events during the accumulation season. As shown in this graphic, metrics such as DSO, DSD, SCD, SWE_max and 1 April SWE could all remain unchanged, but it is clear that the snowpack is different than in the past. Previous literature has quantified increasing ablation during the accumulation period by defining a “melt fraction” (Musselman et al., 2021), which is the ratio of the melt that occurs during the accumulation phase to the total melt. Their metric helps to identify snowpacks that have considerable variance and vulnerability to warming and rain-on-snow events. Another example of changing snowpack is shown in Fig. 1e. Here, the SCD remains constant due to consistent DSO and DSD, but SWE_max decreases in magnitude, resulting in less snow overall. Finally, Fig. 1f shows a theoretical future in which DSO, DSD, SCD, 1 April SWE and SWE_max all remain constant, but it is clear that there is less snow present throughout the season.

The conceptual SWE curve and the above discussion is focused on a mountain snowpack, a snowpack with a distinct period of steady accumulation up to a SWE_max, followed by a similarly steady ablation season that persists throughout the winter. While mountain snowpacks play a key role in natural water storage, other types of snowpacks also have distinct characteristics and are important to the hydrological cycle. For example, ephemeral snowpacks, snowpacks that tend to be have a lower cold content than mountain snowpacks and come and go throughout the winter (Sturm et al., 1995; Hatchett, 2021), play a role in soil moisture and runoff regimes (Livneh and Badger, 2020; Hamlet and Lettenmaier, 2007). Ephemeral snowpacks tend to experience accumulation and ablation processes nearly in tandem (Liston and Elder, 2006). Alternatively, Greenland and Antarctic ice sheets primarily experience accumulation processes (Liston and Elder, 2006). While metrics such as 1 April SWE, SWE_max and SCD do a good job of characterizing mountain snowpacks, they are not as useful at capturing the transient nature of ephemeral snowpacks or the lack of an ablation season on ice sheets.

This work aims to characterize the extent to which snowpacks serve as natural reservoirs and evaluate spatial and temporal changes in snow water storage in a new, integrated way. As we move into a future of increased climate and snowpack variability, we need snow metrics that can capture diverse and dynamic snowpack regimes. As the majority of natural snow water storage in the US occurs in mountainous regions, it is important to understand how the natural reservoir function of snowpacks is changing across individual mountain ranges. When attempting to quantify snow water storage change, it can be difficult to merge the scale at which most in situ observations are available (at the point/station scale) and the scale at which snowpacks operate as natural reservoirs (at the mountain range scale). This study presents a new spatially and temporally flexible snow metric, snow water storage (SwS), in order to address the reality of changing snowpack regimes and the challenge of spatial variability between snow observations and decision-making. This study examines SwS trends in mountain snowpacks by addressing the following research questions:

1. What are the trends in monthly and annual SwS across the US at discrete point, gridded and aggregated scales?

2. What is the role of mountainous ER3s in US SwS and how has this changed in recent decades?

3. How does SwS relate to other common snow metrics?

2 Methods

2.1 Snow water storage metric (SwS)

SwS quantifies the depth of water stored in snow reservoirs over time and is calculated by integrating the area under the SWE curve:

$$\begin{array}{}\text{(1)}& \mathrm{SwS}=\int \mathrm{SWE}\left(t\right)\mathrm{d}t,\end{array}$$

where SWE has dimensions of length, and integration occurs over a time period (water year, a given month, etc.) of interest. If daily SWE data are used for this calculation at a given point, SwS will have dimensions of meter-days (or md).

As defined above, SwS is a quantity computed at a single point, e.g., a SNOTEL location. However, the SwS metric can also be aggregated across various spatial scales. There are numerous reanalysis products that provide spatially distributed SWE information on a regular grid. In this case, SwS can be computed for a horizontal area (say a particular watershed) of interest, resulting in SwS having dimensions of cubic meter days (m³ d).

SwS can also be computed for various integration periods. If the integration is done over the entire water year, this yields annual SwS (SwS_A). In the integration is for a particular month, this yields monthly SwS (SwS_M). Integrating daily SWE data over a single day produces the daily value of SwS (SwS_D), but this is simply is the same as daily SWE.

SwS is the integrated area under the SWE curve, indicating the cumulative meter-days of water stored as a snowpack. SwS quantifies the degree to which a snowpack functions as a water storage reservoir. Unlike 1 April SWE or SWE_max, SwS can be applied to mountain, ephemeral or permanent snowpacks. Unlike other storage-related metrics such as SSI and WTI, SwS is directly measurable and does not require the combination of multiple data streams to calculate. Ultimately, this integrated metric helps us to understand how much water is held in our snow reservoirs and for how long.

2.2 Data

Daily observations of SWE were obtained from Natural Resources Conservation Service (NRCS) Snow Telemetry (SNOTEL) stations (Serreze et al., 1999) and from Cooperator Snow Sensors (COOP). The SNOTEL network provides data at discrete scattered points across the western US, while the COOP stations used in this study provide data across California. This study used the 465 stations that have a period of record from at least water year 1982 to water year 2020 with less than 10 % of days missing during that period.

Figure 2 Map of ER3s in the US. Mountainous ER3s are colored and labeled.

This study also uses the University of Arizona SWE (UASWE) dataset (Zeng et al., 2018; Broxton et al., 2019), a daily 4 km gridded dataset that spans the US. The UASWE dataset assimilates SWE and snow depth observations into an empirical temperature index snow model that is forced with PRISM temperature and precipitation data (Daly et al., 2008). The primary value of this dataset is that it provides SWE estimates at locations other than the SNOTEL stations. This allows for the aggregation of SWE information over areas of interest (Zeng et al., 2018). The UASWE product has been shown to outperform (Dawson et al., 2018) other gridded SWE products such as the SWE estimates from the Snow Data Assimilation System (SNODAS; NSIDC, 2004). Additionally, a spatially continuous, gridded product allows us to build a more complete picture of spatial changes in SwS and how changes in SwS are occurring at aggregated scales. Despite the good performance of the UASWE product, there are limitations to using any modeled SWE product. Factors such as imperfect physics, inaccurate boundary conditions and rescaling errors can contribute to inaccuracies in the modeled field, SWE in this case (Sturm, 2015; Zhang et al., 2007).

The Environmental Protection Agency (EPA) Level III Ecoregions (ER3s) (McMahon et al., 2001; Omernik and Griffith, 2014), which are regions with similar ecosystems and environmental resources, were used to identify mountainous regions and to delineate the grid cells in the UASWE dataset that were associated with each ER3 (Fig. 2). Mountainous ER3s were included in this study if at least half of their area resided in the snow-covered mask (described in Sect. 2.3 below). As each ER3 has similarities in biotic, abiotic, terrestrial and aquatic ecosystem components, examining SwS change in any given ecoregion may help us understand ecosystem impacts that are related to changes in SwS. Numerous ER3s correspond to the major mountain regions in the western and eastern US that serve as the largest natural reservoirs in the country.

Finally, NASA Shuttle Radar Topography Mission (SRTM) digital elevation data (Farr et al., 2007) were regridded to create a digital elevation model (DEM) matching the grid of the UASWE product. Elevation data were used to calculate watershed hypsometry in each ER3. The procedure used to calculate a hypsometry grid is described in Sect. 2.4.2.

Although the station and ER3 datasets extend beyond the conterminous US, the UASWE dataset does not. All datasets were spatially constrained to the conterminous US in order to facilitate the comparison of results between spacial scales. A summary of all of these datasets is provided in Table 1.

Table 1 Summary of data used in this work. n/a stands for not applicable.

2.3 Study area

As noted above, this study considers both discrete station data that focus on the western US and spatially continuous gridded data that cover the conterminous US. Regarding the gridded data, many locations have little to no snow. Therefore, the analysis of the gridded product is restricted to locations that have a mean of at least 30 snow-covered days per year based on the 39-year climatology (1982–2021) of the UASWE dataset (Fig. 3). As expected, snow cover duration increases with latitude and elevation, with the longest snow cover duration found along mountain tops in the western US. In the ER3 SwS change analysis, all ER3s are considered that contain grid cells that meet the 30 d snow cover threshold, although the mountainous ER3s are more closely examined because they store the bulk of our winter water.

Figure 3 Study area for the UASWE dataset, indicated by color shading, showing the mean number of SCDs across the contiguous US in locations that have a minimum average of 30 snow-covered days per year over the period of study.

2.4 Analysis

2.4.1 SwS trends

To evaluate significant trends in SwS_A and SwS_M across the US, these quantities were computed over a 39-year period of study (water years 1982–2020) at stations and at UASWE grid cells, which have an area of 16 km². The grid-cell-based SWE from the UASWE product was additionally aggregated for each ER3 in order to assess trends at larger scales. To compute SwS at aggregated ER3 scales, the gridded SWE data within an individual ER3 were simply integrated spatially, resulting in SwS with units of cubic meter days (m³ d).

This study used the Hamed and Rao modified Mann–Kendall (MK) test from the pyMannKendall Python package to compute trends in SwS (Hussain and Mahmud, 2019). The MK test is a rank-based nonparametric test that is used to evaluate monotonic (increasing or decreasing) trends in temporally varying data (Hirsch et al., 1982). Thus, the null hypothesis is that the data are randomly and independently ordered, whereas the alternative hypothesis is that a monotonic trend exists in the data. Although the MK test is widely used in hydrological studies, it does not account for positive autocorrelation, which increases the probability of detecting trends when no trends exist. Because of this, many studies have turned to a modified MK test that does account for autocorrelation (Hamed and Rao, 1998).

2.4.2 Trends by elevation in mountainous ecoregions

This analysis focused on 16 ER3s corresponding to the mountain ranges that receive substantial snowfall relative to surrounding ecoregions. A total of 12 of these ER3s are located in the western US and 4 are located in the eastern US. The relative elevation of SwS_A change in each ER3 is examined in this study. In order to make trends in SwS_A comparable over the wide range of elevations across the US, the elevations of each ER3 are converted to hypsometry scores. Each ER3 boundary is used to select co-located elevation data from the regridded NASA SRTM digital elevation dataset. ER3 hypsometry is calculated by determining the percentage of the ER3 area that falls below a given elevation within that ER3. Thus, there is 0 % of the ER3 at the lowest elevation of the ER3 and 100 % of the ER3 is below the highest elevation. Each elevation grid cell in the DEM is turned into a value between zero and one based on where that grid cell lies relative to other elevation grid cells within the same ER3. Hypsometry scores in each mountainous ER3 are then binned into 10 % increments, from 0 % of an ER3 below to 100 % of an ER3 below, in order to compute the mean SwS_A and the percentage change in each hypsometry band from 1982 to 2020. The percentage change in the interquartile range (IQR) of SWE was also computed for each hypsometry band from 1982 to 2020. To calculate the percent change in the IQR, the IQR in each ER3 is calculated for each year in the study by subtracting the 25th percentile from the 75th percentile of SWE. The trend is evaluated in each hypsometry band following the trend analysis described in Sect. 2.4.1.

2.4.3 SwS_A compared with other snow metrics

SwS_A trends are compared with other commonly used snow metrics including 1 April SWE, SWE_max, day of SWE_max and SCD in order to evaluate what type of information the SwS_A metric provides that other metrics do not. This is done in four ways using the station data. First, the percentage of stations with positive, positive significant, negative and negative significant trends in each metric is computed. Second, the trend in the annual number of snow-free periods is calculated at each station from 1982 to 2020 to evaluate whether snowpacks are becoming more ephemeral using the Hamed and Rao modified MK test described in Sect. 2.4.1. The annual number of snow-free periods is defined as the number of times in a water year that there is no snow following a period of snow. Next, the utility of the SwS metric compared with other snow metrics is shown using a case study of SNOTEL station 706 (Quartz Mountain, Oregon), a station that has transitioned from a mountain snowpack to an ephemeral snowpack. Third, a regression is computed between the percentage change in SwS_A and each other metric above using empirical data from the stations. Finally, the relationship between the percentage changes in the empirical data is compared with the percentage changes that would be expected based on the conceptual SWE curve. For example, the empirical relationship between the percentage change in SwS_A and the percentage change in SWE_max is compared with what it would be if there was a uniform scaling in the conceptual SWE curve as shown in Fig. 1b.

3 Results

3.1 SwS change trends

3.1.1 SwS_A change trends

The average SwS_A across all stations in this analysis is 60 md. The lowest SwS_A found at a single station was 0 md, whereas the maximum SwS_A observed was 510 md. Changes in SwS_A range from a decrease of 122 md to an increase of 69 md over the period of study (Fig. 4). Of the 97 SNOTEL and COOP stations with increasing trends in SwS_A, only 10 had significant (p<0.1) increases (Fig. 4). Significant decreasing SwS_A trends were found at 123 of the 367 stations. Losses in SwS_A ranged from 2 to 122 md. Spatially, there are widespread decreasing SwS_A trends across most of the 11 western states that contain snow stations, with declines ranging from 17 % to 87 %. The 10 stations with significant increases in SwS_A range from a 6 % increase to a 78 % increase. The stations with increasing SwS_A trends are mostly located in the Northern and Middle Rockies and also include a few station in the Southern Rockies and in the Cascades.

Figure 4 Change in SwS_A (md) (a) and percentage change in SwS_A (b) across US stations from water years 1982 to 2020. Large outlined circles indicate stations with p<0.1.

Moving from discrete station data to the spatially continuous gridded UASWE data, there is a mean SwS_A of 1.8×10⁸ m³ d across grid cells (Fig. 5). Mountainous ER3s in the western US have an average SwS_A of over 1.4×10⁹ m³ d and the maximum SwS_A is 2.6×10¹⁰ m³ d. The average SwS_A in much of New England and the Upper Peninsula of Michigan ranges from 3.2×10⁸ to 6.4×10⁸ m³ d.

Figure 5 Change in SwS_A (md) (a) and percentage change in SwS_A (b) across the UASWE dataset. Stippling indicates locations with p<0.1.

The grid-cell-scale change analysis yields similar geographic patterns of significant changes in SwS_A in the western US to the station-scale analysis (Fig. 5). This is not surprising given that the UASWE product assimilates SNOTEL (and other) station data. The benefit of including a spatially distributed product such as UASWE in this analysis is that it adds detail and insight as to where changes in SwS_A are occurring beyond the western US and in between the locations where discrete stations are located. Significant increases in grid-cell SwS_A are primarily found in the north central and northeastern US. Only 5 % of US grid cells have significant increasing trends with a mean increase of 84 %. From 1986 to 2015, the north central and northeastern US experienced an increase in annual precipitation, particularly in spring and fall, although these regions also show spatially variable increases in precipitation during the winter (Easterling et al., 2017). These precipitation changes may partially explain the increases in SwS_A, although these regions have also experienced increases in winter temperatures over the same time period. Significant decreases in SwS_A are more widespread and are found across the western US, the Appalachian Mountains, the Blue Ridge Mountains and in the Ozarks. Of the 54 % of US grid cells that have decreasing trends in SwS_A, 11 % have significant decreasing trends. The mean decline in SwS_A for the grid cells with significant trends is 44 %.

Figure 6 indicates the raw change and percentage change in SwS_A across ER3s. Aggregating UASWE SwS_A at ER3 scales spatially filters (and thus mutes) some of the grid-cell-scale trends in SwS_A, as can be seen when comparing Figs. 4 and 5. Of the 51 ER3s that are evaluated in this study, 19 have increasing trends and 32 have decreasing trends. Only one ER3, the non-mountainous Lake Agassiz Plain, has a significant positive trend in SwS_A (86 % increase). None of the non-mountainous ER3s have significant decreasing SwS_A trends. The specific ways in which SwS_A has changed across the 16 mountainous ER3s and how these changes relate to other snow metrics will be discussed in Sect. 3.3.

Figure 6 Change in SwS_A (md) (a) and percentage change in SwS_A (b) aggregated across ER3s over water years 1982–2020. Stippling indicates ER3s with p<0.1.

3.1.2 SwS_M change trends

The highest monthly mean SwS_M occurs in March at stations (12.8 md), grid cells (2.7×10⁷ m³ d) and in mountainous ER3s (7.4 m³ d). Figures 7–9 summarize trends in SwS_M change evaluated at stations, UASWE grid cells and ER3s, respectively. Significant decreases in SwS_M occur across all months at all spatial scales examined. November experienced the highest number of stations, grid cells and ER3s with significant SwS_M losses compared with any other month. The greatest monthly median percentage loss of SwS_M occurred in November at stations (56 %) and grid cells (44 %), while it occurred in March at mountainous ER3s (61 %). Looking at raw change values, the months with the greatest decrease in SwS_M are not the same as the months with the largest percentage changes. March, December and January are the months with largest decrease in median SwS_M at stations (4.6 md), grid cells (6.4×10⁶ m³ d) and mountainous ER3s (1.1×10¹¹ m³ d), respectively.

Figure 7 Changes (significance p<0.1) in SwS_M across US stations, in dimensional units (a) and in terms of percentage change (b). The rectangle indicates the interquartile range, with the middle bar indicating the median. The blue pluses are outlier points. Panel (c) shows the fraction of stations that had significant increases or decreases in SwS.

Figure 8 As in Fig. 7 but for the gridded UASWE dataset.

Figure 9 As in Fig. 7 but for changes in SwS_M aggregated across ER3s from water years 1982 to 2020. Red components of this figure indicate the results only considering mountainous ER3s.

Although there is an overall negative median percentage change in SwS_M in all winter months at the station and grid-cell scale, February and March have higher occurrences of significant SwS_M increase than any other months. At the grid-cell scale, October, March–May and July–August all have a 0 % median change in SwS_M because most grid cells within the snow-covered mask are snow-free during these times. At the ER3 scale, the median percentage change in SwS_M is negative in all months. Most data points that indicate significant positive increases in monthly storage are considered outliers at all spatial scales.

3.2 SwS_A change trends in mountainous ER3s

Analysis of mountainous ecoregions illuminates the large role that mountains play in storing winter snow water resources as snowpack, particularly in the western US. An average of 72 % of the annual SwS_A in the US (3.5×10¹³ m³ d) is held in the 16 mountainous ER3s, despite these ER3s only covering 16 % of the US land area (Figure 6). Western mountainous ecoregions cover 12 % of the US land surface and store an average of 65 % of the annual SwS_A. Across all mountainous ER3s, there has been a 22 % decline (6.1×10¹² m³ d) in SwS_A over the 39-year period of study. Over the same time span, there has been a 24 % decline in SwS_A in western mountainous ER3s, indicating that western snow reservoirs are shrinking faster than eastern snow reservoirs.

Of the 16 mountainous ER3s (outlined in red), only the Northeastern Highlands has a (nonsignificant) increasing SwS_A trend, whereas the other 15 mountainous ER3s have decreasing SwS_A trends. This means that the snow water storage in 94 % of mountainous ER3s has declined from 1982 to 2020. Five of the mountainous ER3s, the Cascades, the Eastern Cascade Slopes and Foothills, the Southern Rockies, the Idaho Batholith, and the Arizona/New Mexico Mountains, have significant decreasing SwS_A trends, with a mean decrease of 38 %.

Table 2 summarizes the fraction of US SwS_A in each mountainous ER3, the percentage change in SwS_A from water years 1982 to 2020 and the p value associated with the change. In the western US, the Middle Rockies are responsible for the greatest fraction (12 %) of SwS_A in the country, followed by the Southern Rockies (10 %) and the Idaho Batholith (8 %). SwS_A has declined in all mountainous ER3s in the western US over the last 39 years. The greatest declines in western SwS were in the Arizona/New Mexico Mountains (56 % decline), the Eastern Cascade Slopes and Foothills (40 % decline), and the Cascades (39 % decline). All eastern mountainous ER3s showed declines in SwS_A over the last 39 years with the exception of one, the Northeastern Highlands. The Northeastern Highlands are responsible for the greatest fraction (4 %) of SwS_A in the eastern US; SwS_A has increased 13 % in this region over the last 39 years. The greatest decline in SwS_A in the eastern US was in the Ridge and Valley (11 % decline), which holds an average of 0.2 % of US SwS_A.

Table 2 Overview of SwS_A in mountainous ER3s. Bold p values are significant at p<0.1.

The greatest SwS_A is found in the highest 10 % of mountainous ER3s elevations (Fig. 10). Most mountainous ER3s have decreasing trends in SwS_A across all hypsometry bins, although the Sierra Nevada (5), the Wasatch and Uinta Mountains (19), the Southern Rockies (21), the North Central Appalachians (62), and the North Cascades (77) show increasing trends in SwS_A at low elevations and decreasing trends at higher elevations. The Northeastern Highlands is the only ER3 that shows increasing trends in SwS_A at all elevations. Increasing trends in SwS_A at low elevations in some ER3s may partially be a result of very low SwS_A to begin with; thus, small changes in SwS_A may suggest large percentage changes.

Figure 10 SwS_A in each mountainous ER3 as a function of ER3 hypsometry. The teal line indicates the average SwS_A in each hypsometry bin from water years 1982 to 2020 (left axis). The solid red line indicates the percentage change in subscript A as a function of hypsometry for each mountainous ER3 over the time period of study (right axis). Black “×” symbols indicate where the percentage change in SwS_A is significant (p<0.1). The dashed red line indicates the percentage change in the IQR of daily SWE as a function of hypsometry in each mountainous ER3 over the time period of study (right axis). Black “o” symbols indicate where the percentage change in the IQR is significant (p<0.1). Refer to Table 2 for ER3 names.

By looking at the percentage change trends in the IQR of SWE, this study gives an idea of how interannual SWE variability has changed from 1982 to 2020 (Fig. 10). Several ER3s have increases in the SWE IQR of the lowest hypsometry bands, which correspond to the lowest parts of ER3s. This could be a result of increasing snow variability as freezing levels move to higher elevations, resulting in increased irregularity in precipitation form. In the middle and upper hypsometry bands of most ER3s, there is largely a decrease in the IQR. This may be a result of declining snowpacks, which would allow for less variability in the range of SWE values overall. The Northeastern Highlands, Ridge and Valley, Central Appalachians, and North Cascades stand out in that they have increasing tends in IQRs across most hypsometry bands.

3.3 SwS_A compared with other snow metrics

Figure 11 demonstrates how the conceptual SWE curve is changing in each of the ER3s based on the trend analysis of four common snow metrics: DSO, SWE_max, SWE_max day of water year (D_max) and DSD. Trends in DSO, SWE_max, D_max and DSD were evaluated because these metrics serve as anchor points that define the boundaries of the conceptual SWE curve. These conceptual SWE curves are superimposed on the observed mean SWE curves for the first and last 20 years of study in each ER3. With the exception of the Northeastern Highlands and the North Cascades, the 2020 SWE curve (solid red line) delineates a smaller conceptual SWE curve than in 1982 (dotted red line). The specific anchor points that cause the shrinkage of the conceptual SWE curve are variable across ER3s. Of the ER3s that experienced significant (p<0.1) decreases in SwS_A, the Southern Rockies and the Arizona/New Mexico mountains also had significant changes in the DSO, DSD and SWE_max. The Cascades, the Eastern Cascade Slopes and Foothills, and the Idaho Batholith also experienced significant declines in SwS_A. These three ER3s had significantly earlier DSDs. Although evaluating change at the anchor points tells us how the conceptual SWE curve is changing, the actual SWE curve is not a triangle and is subject to complex patterns of change including notable accumulation- and ablation-season SWE variability. A comparison of the mean SWE curves from the beginning and end of the study period yields a smaller SWE curve in the last decade of study in all western ER3s, although the eastern ER3s (Northeastern Highlands, North Central Appalachians, Ridge and Valley, and Central Appalachians) have more nuanced change in the SWE curve. These subtler shifts involve certain aspects of the snow season having higher SWE values in the first decade of study, while other parts of the snow season have higher SWE values in the later period of study. The SwS metric accounts for these complex patterns and is able to quantify natural reservoir storage and change.

Figure 11 Observed SWE curves for mountainous ER3s shown with historic and current conceptual SWE curves. The black dotted line indicates the mean SWE curve from the first 20 years of our study period in each ER3 (1982–2002). The solid black line indicates the mean SWE curve from the last 20 years of our study period for each ER3 (2000–2020). Light blue shading indicates the IQR of observed daily SWE and gray shading indicates the minimum and maximum range of observed daily SWE. Red ER3 labels signify ER3s that had a significant decrease in SwS_A over the period of study. The red dotted line indicates the conceptual SWE curve for each ER3 at the start of our study period (1982) based on the trend analysis of DSO, SWE_max, D_max and DSD. The solid red line indicates the conceptual SWE curve for each ER3 at the end of our study period (2020). Red stars on the x axis indicate significant decreases in DSO, D_max and DSD. Red stars on the y axis indicate significant decreases in SWE_max. There were no significant increases in any of the snow metrics. Refer to Table 2 for ER3 names.

Trends in snow metrics (SwS_A, 1 April SWE, SWE_max and SCD) are not changing in the same direction at all stations (Fig. 12). SNOTEL and COOP stations were placed to capture mountain snowpack regimes where snow increases up to a maximum value throughout the accumulation season and then disappears across the ablation season. The 1 April SWE, SWE_max, SwS and SCD measures are snow metrics that are useful for characterizing and monitoring change in mountain snowpack regimes and describe points relevant to a conceptual SWE curve. While trends in these metrics have changed in the same direction at 286 of the 465 SNOTEL and COOP stations, 38.5 % of stations have nonuniform inter-metric trend directions. This means that if trends in only one snow metric were to be examined, it could paint an incomplete picture of change.

Figure 12 Trend directions in SwS_A (a), 1 April SWE (b), SWE_max (c) and SCD (d) across stations that do not have the same trend directions across all metrics.

The inability of common one-dimensional snow metrics to reflect snow storage change is particularly apparent when snowpacks transition from one snow regime to another, such as a permanent snowpack transitioning to a mountain snowpack or a mountain snowpack transitioning to an ephemeral snowpack. In our observational record, this study finds that snow regimes are becoming increasingly ephemeral at many station locations. This study found the number of annual snow-free periods has significantly increased at 23 % of stations over the last 38 years. The Quartz Mountain (OR) SNOTEL station (706) provides an example of where a snowpack has increased in ephemerality (Fig. 13). This station went from an average of 1.7 snow-free periods per water year over the first decade of study to an average of 6.3 snow-free periods per water year over the last decade of study. A side-by-side comparison of trends in SwS_A, 1 April SWE, SWE_max and SCD illustrates how the SwS metric is able to capture complex change, such as that associated with an ephemeral snowpack. At the Quartz Mountain station, there are significant negative trends in SwS and SWE_max, no trend in 1 April SWE, and a positive trend in SCD (Fig. 13). In this example, 1 April SWE is largely not relevant as a snow monitoring metric because the majority of years are snow-free on 1 April at this location. This station is interesting because it has opposite significant trends in SWE_max, which is decreasing, and SCDs, which are increasing. As SwS_A is the integral of the SWE curve, both the magnitude and duration of snow cover are incorporated into its calculation. This allows the SwS metric to provide a robust picture of the degree to which the snowpack is serving as a reservoir for water storage and how that reservoir function may be changing.

Figure 13 SWE curves from every second water year from 1982 to 2020 at the Quartz Mountain SNOTEL station.

There are other ways to demonstrate the utility of SwS as an additional tool for gaining insight into our changing snowpack. The 1 April SWE is overwhelmingly the metric cited by water resource managers as the singular measure of the season's snow. However, how do changes in that measure correspond to changes in others? First, Fig. 14 shows the relationships between percentage changes in various snow metrics of interest. It can next be useful to return to the idea of a conceptual SWE curve (Fig. 1a). Based on this simple geometry, the SwS_A is given by

$$\begin{array}{}\text{(2)}& {\mathrm{SwS}}_{\mathrm{A}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{SWE}}_{\max }\mathrm{SCD}.\end{array}$$

If the geometry of the conceptual SWE curve were to be preserved over time, the SWE curve would be uniformly scaled (Fig. 1b). In this scenario, one would expect the percentage change in SWE_max and SCD to each be half the percentage change in SwS. However, the regression plots in Fig. 14 reveal that the change in SCD is roughly 26 % that of the change in SwS. Furthermore, it is found that the change in SWE_max is 86 % that of the change in SwS_A. What this means is that the conceptual SWE curve has been flattening over the period of study of the data at stations. Therefore, relying on a single metric like 1 April SWE gives an incomplete assessment of the storage of snow throughout a full season, and a more holistic metric like SwS may be more informative when considering a full snow season.

Figure 14 Regression of percentage change in SwS_A with percentage change in 1 April SWE (a), SWE_max (b) and SCD (c). Regression of SWE_max with SCD (d).

4 Discussion

SwS is a unique snow metric because it essentially has unlimited degrees of freedom – any change in the SWE curve (changes in SWE magnitude, timing, variability, etc.) will be captured in its calculation. Thus, regardless of how the SWE curve changes, the SwS metric is able to provide information about water storage in natural reservoirs at a given location. SwS is different from the other common snow metrics discussed in this paper (1 April SWE, SCD and SWE_max) because the other metrics have 1 degree of freedom: they provide data on one dimension of the SWE curve. The flexibility of the SwS metric is particularly useful when attempting to quantify storage in snowpacks that have fundamental shifts in their SWE curve. As demonstrated above, the SwS metric can still quantify natural reservoir storage when a snowpack has transitioned from a mountain-type snowpack to an ephemeral type. If a snowpack were to transition from a permanent snowpack to a mountain-type snowpack, the SwS metric would also be able to provide information on the storage.

Of the snow metrics discussed in this paper, SwS is uniquely positioned to capture storage change at aggregated scales, across the full SWE curve. Metrics such as DSO, SWE_max, DSD and SCD are essentially anchor points for the conceptual SWE curve. Although it can be beneficial to note changes at any of these points, change can also happen in between these anchor points, such as an increase in accumulation-season ablation or an increase in ablation-season accumulation, as is seen across mountainous ER3s in Fig. 11. This type of change has already been documented as changing melt fractions (Musselman et al., 2021). At aggregated scales, such as across a watershed, snowpack variability (due to landscape features such as aspect or elevation) influences the SWE curve between the anchor point metrics. Change that occurs between the anchor point metrics is not inherently captured by these metrics. As the SwS metric can accommodate various spatial scales, it is able to capture natural reservoir storage regardless of variability in snow change or snowpack fluctuations that occur in between anchor point metrics.

The widespread losses of SwS_A over the last 39 years reported in this study are consistent with the broader narrative of snowpack change literature, which has established declines in snow-covered area, snow cover duration, 1 April SWE and SWE_max, among others (Rupp et al., 2013; Mote et al., 2018; Notarnicola, 2020; Marshall et al., 2019; Bormann et al., 2018; Choi et al., 2010; Huning and AghaKouchak, 2018). Losses of winter snowpack are largely attributed to increasing global temperatures (Hamlet et al., 2005), which have resulted from a combination of natural variability and anthropogenic climate warming (Rupp et al., 2013; Pederson et al., 2013). The declining trends in SwS_A are a reflection of declining trends in SwS_M in nearly every month, at every scale. The greatest percentage losses in SwS_M occur early in the snow season, particularly in November. The loss of early-season SwS is consistent with previous work that used satellite imagery and reported that DSO is occurring later (Notarnicola, 2020). While future work could explore the exact mechanistic drivers of predominantly decreasing SwS_A trends, these findings are reasonable in the context of mechanistic drivers explored in other snow change literature. From an energy budget standpoint, snow falling at warmer temperatures (as a result of climate warming) and overall shallower snowpacks (due to reduced snowfall fractions) contribute to reduced cold content and more readily ripening snowpacks (Jennings and Molotch, 2020). Additionally, shallower snowpacks are susceptible to enhanced snowmelt from the albedo feedback as vegetation and soil are exposed (Kapnick and Hall, 2010). Although the majority of SwS_A trends are declining, the North Central Plains and the Northeastern Highlands show increasing SwS_A trends. Snowmelt and rain-on-snow are known to be flood-generating mechanisms in New England, Minnesota, and along the Mississippi and Missouri rivers (Collins, 2009; Novotny and Stefan, 2007; Wiel et al., 2018; Olsen et al., 1999). The increasing SwS trends in these regions may, therefore, have implications for flood hazards.

Spatial scale has long been a topic of conversation in snow hydrology, as certain processes that occur at very small scales contribute to considerable within-grid-cell heterogeneity as one scales up from point to grid cell to regional scales (Blöschl and Sivapalan, 1995; Molotch and Bales, 2005). This works finds differences in the magnitude and timing of significant changes in SwS_A and SwS_M when different spatial scales are compared. For example, less of the US landscape shows significant changes in SwS_A in the ER3 analysis compared with the grid-cell analysis. Thus, the aggregation of SwS_A into ER3s filters some of the grid-cell-scale spatial SwS_A trends. Temporally, there is a higher fraction of sites with a significant positive increase in SwS_M from October to March in the grid-cell analysis compared with the ER3 analysis. This indicates that local significant increases in SwS_M at grid-cell scales are offset by smaller-magnitude increases in SwS_M or decreases in SwS_M at many locations once the SwS_M is aggregated to ER3 scales. From a water resources perspective, these findings underscore the importance of choosing an appropriate aggregation scale in order to accomplish management goals.

In the western US, where snowmelt is vital to supplement warm-season water supplies, about 70 % of runoff in mountainous regions originates as snow (Li et al., 2017). Snowpacks also play an important role in climate, ecological processes and recreation in both eastern and western mountains. Across all mountainous ER3s, there has been a 22 % decline in SwS_A over the 39-year period of study. The ER3 mountain ranges considered in this work include the headwaters of 13 of the 18 water basins located in the US, underscoring the importance of these natural reservoirs to water resources. The loss of SwS in these regions is of further concern, as the warm season is projected to increase in length due to anthropogenic climate warming (Mallakpour et al., 2018; Padrón et al., 2020; Siler et al., 2019). Siler et al. (2019) also suggest that declining snow trends may accelerate once the current natural climate mode changes because natural variability has slowed the decline in western snowpacks since the 1980s. The capacity of natural snow reservoirs is declining in most of the western US and across most mountain ranges in the US. As declining trends are expected to continue into the future, monitoring our natural snow reservoirs is essential. Metrics that are highly flexible in space and time (like SwS) can be used to monitor change and evaluate future projections.

SwS_A, the changes in SwS_A and the variability in SWE are all influenced by elevation. The greatest amount of snow water storage occurs on a disproportionately small fraction of our landscape – at the highest elevations of mountainous ER3s. Almost all mountainous ER3s are losing SwS_A at all elevations. In the majority of mountainous ER3s, the highest elevations have experienced the greatest SwS_A losses over the last 39 years. The elevation-dependent changes in our natural snow reservoirs are likely associated with documented elevation-dependent changes in temperature and precipitation (Wang et al., 2014; Harpold et al., 2012; Pepin et al., 2015, 2022; Qixiang et al., 2018). Winter temperatures have increased significantly in the recent past (Vose et al., 2017), which increases the vapor pressure deficit in the atmosphere and may enhance sublimation and vapor fluxes (Harpold et al., 2012). Higher elevations have also warmed at faster rates than their low-elevation counterparts, where there have been increasing trends in precipitation (Wang et al., 2014; Pepin et al., 2022). Wang et al. (2014) suggest that elevational warming amplification is likely associated with effective moisture convection. These mechanistic drivers are a plausible explanation for finding the greatest SwS_A loss at the highest elevations.

This work also finds elevation-dependent changes in SWE variability. SWE variability has likely increased as a result of winter freezing levels moving to higher elevations (Catalano et al., 2019), an increased fraction of precipitation falling as rain instead of snow and more rain falling on snow (McCabe et al., 2007), all of which are related to increasing winter temperatures. Decreases in SWE variability at higher elevations, where there are declining trends in SwS_A, may be a result of shallower snowpacks overall.

5 Conclusions

A new snow metric, SwS, is defined and used to identify where and to what extent water storage in natural snow reservoirs has already changed in the observational record. Mountains, especially western mountains, play a disproportionate role in natural water storage relative to the surrounding landscape. High-elevation natural snow reservoirs are responsible for the greatest SwS_A and have generally experienced the greatest declines in SwS_A. Declines in SwS_A are associated with a fundamental shift in the shape of the conceptual SWE curve, as it appears to be flattening across stations. As we move into a future of increased snow variability and diminished snowpacks and as more of the winter snow landscape transitions to ephemeral regimes, temporally static metrics such as 1 April SWE and SWE_max may become less representative of our snowpacks. Concurrently, it may be useful to have metrics such as SwS that can adapt to a wide range of circumstances. Spatially and temporally flexible metrics such as SwS may become increasingly valuable, particularly when it comes to monitoring change.

Declining storage in our natural snow reservoirs has broad implications for human and ecological systems. Natural snow reservoirs help to increase water storage far beyond the capacity of artificial reservoirs in the western US, supporting their roll in linking cool-season precipitation to warm-season water demand. As one of the most robust projected impacts of climate change is a continued increase in air temperatures, it is likely that declining trends in SwS_A will continue. Water managers, planners and decision-makers will need to account for these declines in natural snow water storage as they relate to streamflows for fish migration and recreation, municipal and agricultural water supplies, and flood hazards. Although this paper does not focus on future predictions of snowpack, SwS could be a useful tool for understanding how our natural snow reservoirs change in the future. Change in our natural snow reservoirs is multidimensional and already happening. Metrics are needed that can capture this complexity of change.