Next Article in Journal
Estimating the Precipitation Amount at Regional Scale Using a New Tool, Climate Analyzer
Previous Article in Journal
Nutrient Atmospheric Deposition on Utah Lake: A Comparison of Sampling and Analytical Methods
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Snow Water Equivalent Accumulation Patterns from a Trajectory Approach over the U.S. Southern Rocky Mountains

by
Isaac J. Y. Schrock
1,
Steven R. Fassnacht
2,3,4,*,
Antonio-Juan Collados-Lara
2,5,
William E. Sanford
6,
Anna K. D. Pfohl
2 and
Enrique Morán-Tejeda
2,7
1
Civil Engineering, Colorado State University, Fort Collins, CO 80523-1372, USA
2
ESS-Watershed Science, Colorado State University, Fort Collins, CO 80523-1476, USA
3
Cooperative Institute for Research in the Atmosphere, Fort Collins, CO 80521-1375, USA
4
Natural Resources Ecology Laboratory, Fort Collins, CO 80523-1499, USA
5
Instituto Geológico y Minero de España (IGME, CSIC), 18006 Granada, Spain
6
Geosciences, Colorado State University, Fort Collins, CO 80523-1482, USA
7
Departamento de Geografía, Universidad de las Islas Baleares, 07122 Palma, Spain
*
Author to whom correspondence should be addressed.
Hydrology 2021, 8(3), 124; https://doi.org/10.3390/hydrology8030124
Submission received: 16 July 2021 / Revised: 13 August 2021 / Accepted: 16 August 2021 / Published: 18 August 2021

Abstract

:
The spatial characteristics and patterns of snow accumulation and ablation inform the amount of water stored and subsequently available for runoff and the timing of snowmelt. This paper characterizes the snow accumulation phase to investigate the spatiotemporal snow water equivalent (SWE) distribution by fitting a function to the trajectory plot of the standard deviation versus mean SWE across a domain. Data were used from 90 snow stations for a 34-year period across the Southern Rocky Mountains in the western United States. The stations were divided into sub-sets based on elevation, latitude, and the mean annual maximum SWE. The best function was a linear fit, excluding the first 35 mm of SWE. There was less variability with SWE data compared to snow depth data. The trajectory of the accumulation phase was consistent for most years, with limited correlation to the amount of accumulation. These trajectories are more similar for the northern portion of the domain and for below average snow years. This work could inform where to locate new stations, or be applied to other earth system variables.

1. Introduction

Across the western United States, human existence and economic activity are defined by the availability or scarcity of water [1]. With persisting climate change, water resources are becoming more in one of the two extreme states: drought or flood [2]. Urban areas are especially susceptible to both drought and flood [3]. Many of these areas are downstream of mountains, where snow is persistent in the winter [4], with about 22% of the world’s population living downstream of the mountains [5].
In the Southern Rocky Mountains (SRM; Southern Wyoming, Colorado, and Northern New Mexico), the snowpack is the dominant contributor to surface water flows, and in Colorado, surface water was estimated to comprise 85% of all water withdrawals [6]. Drought conditions have been identified over the southern portions of the SRM (Figure 1) since the early 2000s [7], and these have persisted for the past two decades [8]. As well, a decline in the snowpack has occurred over the past few decades across most of the SRM [9,10,11]. The changing climate may not influence the snowpack accumulation, melt and subsequent runoff uniformly [12].
Since these drier conditions are leading to a more uncertain future [13,14], an understanding of the spatial distribution of the snowpack is vital for management of water resources [15]. Spatial patterns of snow distribution have been considered at various scales [16] from local [17] to regional [18]. The distribution of snow dictates the rate and timing of snowmelt and the subsequent streamflow. Evaluation of the distribution of snow and terrain and land cover/canopy data has often focused on small domains (km) at fine resolution (m) [19,20,21,22,23,24], with fewer studies examining medium [25] or larger scales [26]. A variety of approaches exist to estimate SWE spatially [27].
Egli and Jonas [18] characterized the dynamics of spatiotemporal snow depth distribution for both accumulation and ablation phases based on six years of data from 77 weather stations across the Swiss Alps. They applied analytical methods originally developed in the field of statistical physics by Barabási and Stanley [28] to describe snowpack depth as a growing then diminishing surface over time. To visualize and analyze snow depth and variability across the domain, a statistical technique was adapted from Crow and Wood [29] and soil analyses by Famiglietti et al. [30] to present snow depth distribution, as the standard deviation, as a function of mean snow depth. The technique allowed for the identification and parameterization of mean snow depth distributions for accumulation and ablation phases.
This study applied the assumptions and methods developed by Egli and Jonas [18] for the Swiss Alps to the headwaters region in the Southern Rocky Mountains (SRM) in the western United States, considering several adaptations. Across the SRM region (Figure 1), there are different snow climatologies [31], with large spatial and elevation variability [32], but a lower spatial density of snow measurement stations compared to the Alps study. While Egli and Jonas [18] used snow depth, this work focuses on snow water equivalent (SWE) data due to their length of record across the SRM. Specifically, the inter-annual variability in accumulation characteristics could be characterized due to the length of record; as well, the latitudinal and elevational differences could be assessed. The objectives are as follows: (1) compare the statistical fit of snow depth versus SWE trajectories, (2) identify regression models and parameters to approximate the accumulation phase of SWE trajectories over the period of analysis, and (3) identify characteristics and patterns across the SRM by sub-dividing the domain by latitude and elevation. A method is evaluated to automatically differentiate between accumulation at all locations and mixed accumulation–ablation. The maximum peak SWE of each snow year was also considered. The focus in this paper is the accumulation phase across all stations; the inter-annual differences in the mean rate of accumulation as a function of snowpack variability provide insight into snowfall patterns and how established stations measure those.

2. Study Area and Data

This work, in the SRM region, builds upon previous snow hydrology studies using the same dataset, in particular the derivation of snow-cover depletion curves [33], the spatiotemporal variability of snowmelt factors [34], the evaluation of snowmelt rates compared to precipitation as potential for flood risk [32], and snowmelt modeling [35].
This paper examined 34 water years (1 October 1981–30 September 2015) of SWE data across the SRM. This time series contains 90 Snow Telemetry (SNOTEL) operated by the Natural Resources Conservation Service (NRCS) [32] (Figure 1). Stations within the study domain stretch from approximately 36 to 42.5 degrees North and 105 to 108.5 degrees West within the states of Colorado, northern New Mexico and southern Wyoming and across the elevation range of 2268 to 3536 m (Figure 1). Mean daily SWE values recorded between water years 1982 and 2015 were obtained from the NRCS <wcc.nrcs.usda.gov> (accessed on 16 August 2021). For these years, the available period of record for the stations was between 28 and 34 years. Quality-controlled SWE data were obtained from Fassnacht and Records [32]. Snow depth data have been available at most of these stations since the late 1990s to late 2000s; thus, there is only about half of the period of record. A 30 m digital elevation model (DEM) and the 2011 land use and land cover (LULC) map were obtained from the U.S. Geological Survey National Map <nationalmap.gov> (accessed on 16 August 2021). For each station, the elevation and land cover type were extracted from the DEM and LULC map, respectively, using the station coordinates obtained from the NRCS. Due to the potential for snow cover canopy interception, land cover was classified as either evergreen forest (50 stations) or non-evergreen forest [34]. The non-evergreen land cover types were predominantly deciduous forest (12 stations), mixed forest (7), or grassland/herbaceous (15).

3. Methodology

3.1. Station Sub-Classification by Latitude, Elevation and Snow Year

Quality-controlled water year data yielded a daily SWE value for each calendar day for each of the 34 water years [32]. The mean SWE and standard deviation of SWE were computed for the entire domain and the two groups divided by latitude (north and south). Each of the two latitude groups was further divided by elevation yielding north-high, north-low, south-high, and south-low (Figure 1a), resulting in a total of seven unique mean and standard deviation SWE values for each calendar day. The latitude-based group divided the SNOTEL stations across the SRM at 38.75 degrees (see Figure 1a) north latitude into north and south sub-sets, based on regions of homogeneity identified with self-organizing map analyses [31] and correlogram analyses of accumulation slopes [36]. The latitude- and elevation-based group further subdivided the north–south SNOTEL station sub-sets by an elevation threshold of 2800 m in the north and 2900 m in the south; the “high” group was above this threshold elevation. Elevation thresholds were selected to the nearest 100 m and to have a similar proportional distribution of low and high elevation stations for the north and south sub-sets, to ensure that there were enough stations in each sub-domain so that statistics could be adequately computed. For each spatial sub-set, all 34 snow years were classified according to the annual maximum daily SWE measurement into three groups: above average (AA, more than ½ standard deviation deeper than the mean), average (AVG, within (±) ½ standard deviation of the mean) or below average (BA, less than ½ standard deviation shallower than the mean) (Figure 2). For assessing the regression model, seven characteristic snow years were selected (above average: 1993, 1995, 2011; average: 2010; below average: 2002, 2012, 2013; highlighted in Figure 2). Since the main purpose is to evaluate the inter-annual variability of accumulation slopes across the entire domain and sub-domains, the seven years were assumed to be representative of snowpack conditions; only one average snow year was selected for the evaluation, as above and below the average year represent the outliers, which tend to be more difficult to adequately fit with statistics.

3.2. Derivation of SWE Standard Deviation versus Mean Trajectories

For each station (i) in the study domain, and day (t) in the snow year, the daily SWE value was represented by xi(t). Averaging these values across a single sub-set yielded x i ( t ) ¯ , or the mean daily SWE. For each sub-set of mean SWE values, the daily standard deviation σ(xi,t) was computed, as per Egli and Jonas [18]. Thus, for each calendar day (t) from 1982 to 2015, and each sub-set of SNOTEL stations, the mean and standard deviation SWE were plotted, yielding 238 unique instances over the 34-year period of analysis.
Egli and Jonas [18] utilized this methodology “to describe the evolution of the snow cover as a growing surface”, based on work by Barabási and Stanley [28] and Löwe et al. [37]. When the snow cover ( x i ( t ) ¯ , SWE) is considered to be a growing/diminishing surface over the snow year, the standard deviation σ(xi,t) within the data set reflects the magnitude of spatial variability. A unique standard deviation SWE value was correlated to each mean SWE value. Each plot of σ(xi,t) versus x i ( t ) ¯ has three components: All stations are accumulating (growing surface), all stations are melting (diminishing surface), and a mixture of accumulation and melting stations, with inflection points in between each set (see Figure 3 as an example). The approach enabled the seasonal SWE distribution, or snow depth, to be modeled across multiple spatiotemporal scales, independent of time.

3.3. SWE versus Snow Depth Trajectories

Egli and Jonas [18] used snow depth to model the trajectories. Here, SWE data were used since those data were available for a longer period of record and at a higher spatial density than snow depth for the SRM study domain. SWE and snow depth data were evaluated across the full domain for three years (2010 as an average snow year, 2011 as a high snow year, and 2012 as a low snow year) to compare accumulation phase regression statistics and slope coefficients between SWE and snow depth trajectories. For all three snow years, the daily standard deviation data over 90 (87) SWE (snow depth) stations across the full domain were plotted against the daily mean over all 90 (87) stations to develop both SWE and snow depth trajectories. The inflection point between accumulation and mixed periods (called hysteresis by Egli and Jonas [18]) was identified visually and used as an upper bound for regression models. Linear regression equations were applied to the accumulation phase, regression statistics (slope coefficients, coefficients of determination) were extracted, and differences between the regression data were evaluated.

3.4. Accumulation Trajectory Models

To characterize and describe snowpack as a singular growing surface, both the regression model (e.g., linear, power) and inflection point parameters (e.g., slope threshold) were investigated. Standard deviation versus mean trajectories from seven characteristic snow years in the full domain sub-set were evaluated to identify the optimal regression model and inflection point parameters. Of the seven characteristic snow years, one to three years of data were obtained from each snow year magnitude category to represent BA (2002, 2012, 2013), AVG (2010) and AA (1993, 1995, 2011) snow years (Figure 2). Snow years were selected based on the presence or absence of distinct trajectory characteristics, such as the curvature and shape of the transition between accumulation, mixed, and ablation phases, and the relationship between the maximum standard deviation SWE and maximum mean SWE. Special consideration was given to select trajectories with distinct behavior.
A sensitivity analysis (Appendix A and Appendix B) was conducted to identify parameters that enable automated detection of the accumulation phase inflection point, i.e., the point between accumulation and mixed accumulation (Figure 3; Appendix A). For the seven sample years (Figure 2), this inflection point was manually identified for comparison to the automated method. A moving-slope difference approach was used, and three different variables were tested: indicator position, magnitude difference in slope, and last feasible calendar date that all stations could be in accumulation phase. For every data point (day) in the domain, linear regression variables were calculated based on the adjacent data points to determine the 10-day or 11-day average accumulation slope. A 10-day linear regression was conducted for the leading (9 days prior, present day) and lagging (present day, future 9 days) indicator positions, and an 11-day linear regression was conducted for the central indicator position (5 days prior, present day and future 5 days). The difference in slopes between consecutive days was calculated. Large differences in slope were assumed to correlate with the end of the quasi-linear accumulation phase. Two thresholds (0.3 and 0.5) were tested to calculate the difference in slope, compared to the manually extracted inflection. Multiple dates for the feasible calendar date for all stations to be in the accumulation phase were also tested. The feasible calendar date was included to maintain data quality.
Inflection points obtained from all six permutations of slope difference threshold and indicator position variables were compared to the ideal inflection point, which was identified by visual inspection for the seven characteristic snow years (Appendix A). The Nash–Sutcliffe efficiency coefficient (NSE) was calculated to measure the ability of each set of model parameters to replicate the ideal inflection point. Linear regression was applied to the observed ideal inflection point versus the modeled inflection point, and the coefficient of determination (R2) was calculated.
After the optimal inflection point parameters were identified, three regression methods were applied to the same seven years’ standard deviation versus mean SWE domain data. Linear, truncated linear, and power regression models were considered (Appendix B), based on prior work by Egli and Jonas [18] and Pomeroy et al. [38]. For both the linear and power regressions, the entire accumulation data set from initial non-zero values through the inflection point were considered, and coefficients M and B and α and β were identified by the following equations:
σ ( x i , t ) = M   ×   x i ( t ) ¯ + B
σ ( x i , t ) = α   ×   x i ( t ) ¯ β
A truncated linear regression model was developed specifically to exclude SWE data from the beginning of the accumulation phase, characterized by non-linear patterns and significant moving-average slope variability. The non-zero initial point (first occurrence of mean SWE ≥ 35 mm) was selected by reviewing the standard deviation versus mean SWE trajectories trends visible on all seven characteristic years, with a focus to exclude the initial non-linear snowfall accumulation phase (Figure 3).

3.5. Evaluation of Accumulation Trajectories and Best-Fit Equations

After the ideal regression model and parameters were identified, both were applied to the accumulation phase of 238 standard deviation versus mean SWE trajectories (one trajectory per year in seven sub-sets). The iterative analysis was conducted with Python code, leveraging SciPy [39] and pandas [40] packages. For each trajectory, slope and coefficient of determination regression statistics were obtained. Within each sub-set, maximum, minimum, average and standard deviation statistics based on the 34 slope values were calculated and each snow year was identified as BA, AVG and AA. The upper and lower boundaries used to categorize the snow year were developed uniquely for each sub-set as described above in Section 3.2 based on the mean and standard deviation of the sub-set.

4. Results

4.1. Evaluation of SWE versus Snow Depth Trajectories

The standard deviation versus mean snow depth and SWE trajectories for the three characteristic snow years (2010, 2011, and 2012) have similar shapes (Figure 4). Within the accumulation phase, the trajectory can be broken into two distinct sections for both SWE and snow depth: a non-linear, hysteretic section for low mean snow depth/SWE values (typically less than 35 mm mean SWE), and a linear section with a positive average slope (Figure 4). Snow depth values have a higher degree of scatter along the trajectory than SWE (Figure 4). The degree of scatter is reflected in the regression statistics of the accumulation phase; coefficient of determination values for SWE are consistently close to 1 (linear), while coefficient of determination values for snow depth values are consistently lower (Figure 4). Snow depth trajectories exhibit abrupt changes in the magnitude and direction of slope that are not seen for SWE during the accumulation phase. The linear accumulation slopes for SWE were greater than snow depth by 24% in 2010 and 2012, and 31% in 2011 (Figure 4).

4.2. Optimal Accumulation Trajectory Model

The linear, truncated linear, and power regression models each fit the measured accumulation phase data (three characteristic years across the full domain shown in Figure 4; Table A1) effectively. Coefficient of determination values are reported to three decimal places to differentiate the fit (Table 1). The R2 values for the two linear models are closer to a perfect fit (1) than for the power function. By excluding SWE data below 35 mm, the linear model fit improves, and the average R2 is 0.994 for the seven characteristic years.
The slope difference threshold value of 0.5 did not lead to detection of an inflection point for the 2010 snow year (AVG), influencing the negative NSE and low R2 values for linear regression between observed and modeled inflection points (Table 2 and Table A2). In contrast, the slope difference threshold of 0.3 led to the detection of an inflection point for all seven years (Table A3). With a slope difference threshold value of 0.3, both the central 11-day and lagging 10-day indicator position parameters reflected the observed inflection point (NSE of 0.85, 0.81, respectively), and exhibited a more linear relationship between observed and modeled mean SWE (R2 of 0.86, 0.87, respectively), than the leading 10-day indicator position parameter. The central 11-day (0.3 slope difference threshold) was utilized in regression model identification with the truncated linear equation, since this combination of parameters yielded the most accurate results for the seven characteristic snow years.

4.3. Accumulation Trajectories and Best-Fit Models

Out of the 238 trajectories, the model did not yield an inflection point or a realistic inflection point for nine trajectories (3.7%): full domain for 2000, south for 2000 and 2012, north-high for 1984 and 2001, south-high for 1987 and 2000, and south-low for 2000 and 2008 (Appendix C). For these cases, the inflection point was obtained by visual inspection. Realistic accumulation inflection points were assumed to occur prior to May 15th.
The 34-year average of SWE accumulation slopes across all seven sub-sets range from 0.36 (south-low) to 0.40 (full domain) (Table 3). The maximum slope was observed in the full domain sub-set in 2006 (slope = 0.59), and the minimum slope was observed in the south-low elevation sub-set in 2007 (slope = 0.12). Across the four latitude and elevation based sub-sets, a slight decreasing trend in average slope was observed (0.01 change per sub-set), from highest and northernmost (north-high = 0.39) to lowest and southernmost (south-low = 0.36). However, average slopes were fairly consistent across all seven sub-sets. The largest range between maximum and minimum slopes was observed in the south-low sub-set (0.40), while the smallest range was observed in the north sub-set (0.18).
The 34-year average R2 and standard deviation values of slope were employed to assess the degree of inter-annual variability (or scatter) on the accumulation phase of the standard deviation versus mean SWE trajectories (Table 3). The south-low sub-set exhibits a high degree of scatter (Figure 5) and had the lowest coefficient of determination (0.95), and highest slope standard deviation (0.097) of the slopes across all sub-sets. Conversely, the north and north-low sub-sets exhibited the highest coefficient of determination (0.99), but the north sub-set exhibited the lowest variation (slope standard deviation of 0.042).
Little to no correlation between scatter and elevation was observed, since separating the south and north sub-sets by elevation into high and low sub-sets yielded equal or slightly lower (0.01–0.02) coefficient of determination values. Conversely, the magnitude of scatter and elevation appear to be correlated, because separating the south and north sub-sets by elevation into high and low sub-sets increased the standard deviation for all four sub-sets. The largest increase in standard deviation was observed in the south-low sub-set, as standard deviation increased 0.020 (from 0.077 (south) to 0.097 (south-low)), while the increase in standard deviation for the other three sub-sets was observed between 0.006 (south to south-high) and 0.011 (north to north-high).
When the accumulation phase data are discretized by snow year magnitude (BA, AVG, AA) and sub-set, several characteristics become apparent (see Figure 2 and Figure 5a for type of snow year). Generally, below average (BA) snow years exhibited lower average accumulation slopes than average (AVG) or above average (AA) snow years (six out of seven sub-sets); the lowest average accumulation slope—regardless of snow year—was observed in the south-high region in the BA snow year category. The standard deviation values were fairly similar within each sub-set, except the south-low sub-set where accumulation slopes were much more variable in BA snow years than AVG or AA snow years.
Between the north and south sub-sets, the average accumulation slopes for BA and AA snow years were fairly consistent, but a large difference in average accumulation slope was observed for AVG snow years. When elevation was considered in the north sub-set (splitting into north-high and north-low sub-sets), no significant impact from snow year category was observed, since the average slope and standard deviation values were fairly consistent. However, when elevation was considered in the south sub-set (south-high and south-low), significant fluctuations in average slope and standard deviation were observed. For example, average slopes for BA snow years decreased, AVG snow years increased, and AA snow years both increased (high elevation) and decreased (low elevation). Discretizing the south sub-set data by elevation revealed standard deviation (scatter) increased for all snow years and both elevation sub-sets, except BA snow years in the south-high sub-set (Figure 5 and Table 3).
The largest difference in average accumulation slope within a sub-set was observed in the south-high sub-set, where the range between the largest (AA snow year = 0.42) and smallest (BA snow year = 0.32) average slopes was 0.10 (Table 3). The smallest difference in average accumulation slope within a sub-set was observed in the south-low sub-set, with a range of 0.02, between AVG (0.37) and BA (0.35) snow years. The range between the largest and smallest accumulation slopes in south-low and south-high sub-sets were in contrast to the ranges observed in the other five sub-sets, between 0.04 and 0.06. Collectively, the results suggest that average accumulation slope values are correlated strongly, moderately, or weakly with the annual peak SWE; south-low has little to no correlation, south-high has strong correlation, and the other five sub-sets have moderate correlation. Stated differently, average accumulation values are fairly consistent for the south-low sub-set, regardless of snow year. Conversely, for the south-high sub-set, average accumulation slopes vary greatly, based on the snow year. For the remainder of the sub-sets, accumulation slopes vary some with the amount of snow each year.

5. Discussion

5.1. Applicability of SWE and Snow Depth for Trajectories

An initial objective of this study was to assess whether SWE or snow depth was better for the trajectory evaluation. The use of SWE instead of snow depth data was motivated by the robustness of the available data. In the SRM, SWE data are available for a longer period of record and higher spatial extent than snow depth data. A second motivation was variability in the trajectory data. Most importantly, from a hydrology and water resources perspective, we ultimately care about SWE [27].
SWE data are a good substitute for and provide several advantages over snow depth data (Figure 4). For the three snow years where regression was applied to snow depth and SWE data (2010, 2011, 2012), the R2 values were closer to 1 for SWE data (Figure 4). The regression fit parameters confirm visual observations from the trajectories (Figure 4); the standard deviation varies less with SWE data than with snow depth data. In the accumulation phase, snow depth data lead to a trajectory with multiple standard deviation values for a given SWE value, which was not observed for trajectories based on SWE data (Figure 4). Egli and Jonas [18] observed the same behavior—high variability and frequent changes in snow depth standard deviation; this behavior was attributed to the physical processes of settling and/or densification, which are measured by snow depth. The SWE trajectory does not exhibit abrupt changes in magnitude/direction or multiple values of standard deviation for a given SWE value, because SWE measures the mass of accumulated water rather than snow depth and decreases in SWE are smaller due to sublimation. The SWE trajectories do not decrease, i.e., the value of the mean SWE, until the start of ablation phase at some stations, because mid-winter (accumulation phase) melt events are rare across the SRM [41]. The use of SWE data eliminated variability due to the physical processes of settling and densification (Figure 3), resulting in an accumulation trajectory to which regression can be applied with greater accuracy than snow depth data (Figure 4).

5.2. Regression Model and Parameters Applied to SWE Trajectories

Although the original work on growing surfaces by Barabási and Stanley [28] and Egli and Jonas [18] elected to model the accumulation phase of the trajectory with a power function, the truncated linear model most accurately fit the accumulation phase in the SRM. For the seven characteristic years in the full domain sub-set, the truncated linear model yielded a coefficient of determination equal to or greater than the linear or power models (Table 1). However, all three models yielded high coefficient of determination values, confirming prior findings by Egli and Jonas [18] and Pomeroy et al. [38] that both power and linear models can be applied to the accumulation phase across different snow measurement variables (snow depth/SWE) and snow hydrology regions (SRM, Swiss Alps).
The truncated linear regression model improves the fit over the non-truncated linear model by excluding mean SWE values below 35 mm from the regression. For all seven years analyzed in depth, the trajectory data below 35 mm were highly variable and had multiple standard deviation values for a given SWE value. The 35 mm threshold was selected to eliminate the majority of the variability and duplicate standard deviation values by visual analysis from the seven snow year trajectories studied in depth (Figure 3). The driver of this behavior is likely caused by both individual SNOTEL stations shifting from accumulation to ablation, and the 90-station aggregate including stations in both accumulation and ablation phases at the onset of the winter season. Similar behavior can be observed during the mixed period between purely accumulation and ablation phases (Figure 4).
Whether generated from snow depth or SWE data, average accumulation slope values are a measure of spatial variability; larger slope values indicate snowpack data are more heterogenous, while lower slopes indicate snowpack data are more homogenous. Both SWE and snow depth accumulation slopes in the SRM were generally lower than the linear regression slopes for snow depth observed by Egli and Jonas [18] in the Swiss Alps. In the SRM, snow depth accumulation slopes ranged from 0.29 to 0.38 (the average over 2010, 2011, 2012 was 0.32), and 34-year average SWE accumulation slopes ranged 0.36 (south-low sub-set) to 0.40 (full domain) (Table 3). Egli and Jonas [18] reported six-year average accumulation slopes between 0.41 and 0.54. Physical processes, such as snow settling, account for the difference in slope between snow depth and SWE data, but do not account for the difference in snow depth slope observed between the SRM and Swiss Alps.
The smaller accumulation slopes in the SRM could be due to differences in data density, period of record analyzed, or snowpack climatology. Differences in snow hydrology between the Swiss Alps and SRM may contribute as well. For example, the annual peak station-averaged snow depth observed by Egli and Jonas [18] in the Swiss Alps was between 1.53 and 2.42 m, while the annual peak station-average snow depth in the SRM was between 0.75 and 3.0 m <wcc.nrcs.usda.gov> (last accessed on 16 August 2021). Egli and Jonas [18] used six snow years that had some variability. In this paper, three years of snow depth data illustrated larger inter-annual variability (Figure 4). In total, 34 years of SWE were considered for the SRM, illustrating much variability (Table 3, Figure 5). Finally, the SRM domain in this paper is an order of magnitude larger than the Swiss Alps region examined by Egli and Jonas [18]. Egli and Jonas [18] utilized data from 77 stations located over approximately 30,000 km2 (390 km2/station), while this study utilized data from 90 stations spread out over approximately 300,000 km2 (3330 km2/station). As such, there is much difference in snowpack accumulation patterns across the SRM [36] that are not seen in the Swiss Alps study domain. These differences across the SRM [31] are especially pronounced between the northern and southern portions of the area (Figure 1). The SRM domain is large enough (>700 km north to south) that snowfall usually arrives from different storm systems across the domain.

5.3. Accumulation Slope Dynamics across the Southern Rocky Mountains

The 34-year average accumulation slopes of standard deviation versus mean trajectories illustrate that similar snowpack accumulation patterns exist across the SRM, which has previously been detailed by Sturm and Wagner [17]. Unique snowpack accumulation patterns were observed across the distinct snow climatologies of the SRM previously identified and characterized by Fassnacht and Derry [31]. The snowpack across the south region was more spatially homogenous than in the north region for individual years, since the 34-year average accumulation slope is lower in the south region than in the north (Table 3 and Figure 5). However, the higher standard deviation for the 34-year average accumulation slope in the south region indicated more inter-annual variability of the snow surface than in the north region (Table 3).
The most numerically significant impact of elevation was observed when the north and south sub-sets were split into high and low sub-sets, because the 34-year average standard deviation values increased for all elevation- and latitude-based sub-sets. While the relatively small number of data points within the twice-divided sub-sets may account for a portion of this change, the results may also suggest substantial snowpack variability exists across sites within similar elevation and latitude ranges. However, the accumulation slopes from the same elevation- and latitude-based sub-sets are very similar to the latitude-based sub-sets. Consistency in the average accumulation slope between elevation- and non-elevation-based sub-sets suggests that for larger spatial scales and data sets, elevation does not dramatically impact spatial or inter-annual snowpack accumulation patterns.
When annual maximum daily SWE for each snow year is discretized into BA, AVG and AA snow years, clear correlations between slope, peak annual SWE, elevation and latitude were observed across the 34-year accumulation phase trajectories. Average accumulation slopes differ somewhat for BA, AVG and AA snow years for the north sub-set, and those differences are amplified in the north-high and north-low sub-sets (Table 3). The consistency suggests that distinct snowpack accumulation patterns exist for different annual maximum SWE categories, and that elevation does not impact snowpack accumulation behavior. BA snow years exhibit the lowest average slope, corresponding to the lowest degree of spatial heterogeneity, followed by AA snow years, and then, AVG snow years with the highest degree of spatial snowpack variability. Within the south-low region, there is more variability in the southern stations and thus, less consistency (Figure 5c). This also suggests peak annual SWE has no impact on the spatial distribution of SWE. However, the large range of standard deviation values in the south-low sub-set reflects a unique characteristic of BA peak SWE snow years—a high degree of spatial variability in snowpack.

5.4. Other Applications

This paper examined 34 years of SWE data (Figure 2) and the difference between high and low snow years can be seen (Figure 5). In the future, if there is less snow in a warmer climate [42], as has been the case in the southern domain since about 2000 [7,8], the application of the rate of snow accumulation trajectory shown herein (Figure 3) can provide insight into accumulation variability across large domains (1000s of km2). It can also be used at finer scales [37], provided the proper balance of spatial and temporal measurements is available. This methodology can easily be extended beyond the United States [43] and used anywhere time series of snowpack (SWE or depth; Figure 4) are available.
Most earth system properties have systematic spatiotemporal variability that is difficult to assess due to a lack of measurements, as is seen with snow [27,44]. We use remote sensing, but in most cases, the resolution is too coarse spatially or temporally to fully assess the dynamic nature of an earth system property [16]. New approaches are combining multiple data sources, such as field measurements, remote sensing, and modeling [27,44]. To extend the approach shown here (Figure 3), more data are needed at finer spatial and temporal resolutions, as well as better models. We must also start to consider how the earth system has been changing [45], how humans interact with the earth [46], how our actions influence it [47], and how policy will impact future changes [48]. For water resource systems, this could provide insight into the occurrence and persistence of droughts and floods [2]. Examining the trajectory models of systematic spatiotemporal hydro-climatic variables could help assess where more monitoring is necessary and where models need to be improved.

6. Conclusions

Linear accumulation slopes for standard deviation versus mean for SWE were greater than snow depth by 24% in 2010 and 2012, and 31% in 2011. Multiple values of standard deviation for a given mean value were more likely on the snow depth trajectory. This hysteretic behavior was attributed to snow depth capturing the physical processes of settling. Snow water equivalent data did not exhibit the same behavior and do not measure these physical processes. While the linear fit to depth was quite good, it was better for SWE.
From seven characteristic snow years, the average coefficient of determination for linear and power models was equal (0.98). Eliminating the initial period (mean SWE < 35 mm) from the linear model yielded an R2 of 0.99. Using the slope difference method, an 11-day moving average and slope difference threshold of 0.3 most accurately predicted the end of the accumulation phase, compared to visual observation.
Across all seven sub-sets, the 34-year average accumulation slopes ranged from 0.36 for the south-low sub-domain to 0.40 for the full domain, suggesting some snowpack variability is different based on latitude. There is less inter-annual variability north of 38.75 degrees (north subsets) and more variability in the south sub-sets, especially below 2900 m in elevation (south-low). For all sub-sets, below-average snow years exhibited lower average accumulation slopes than above average snow years, illustrating that snowpack variability is correlated with maximum annual peak SWE. The differences in the slopes for two regions (north versus south), two elevation subsets (high versus low), and in three magnitudes of snow years (AA, AVG, versus BA) illustrate the nature of spatiotemporal variability in accumulation across the study domain; this can help indicate the transferability of accumulation patterns from one location to another. The approach presented in this paper can be applied to other regions to determine coarser resolution snow accumulation patterns, as well as possible applications to other related datasets.

Author Contributions

Conceptualization, S.R.F. and I.J.Y.S.; methodology, S.R.F. and I.J.Y.S.; software, I.J.Y.S.; formal analysis, I.J.Y.S. and S.R.F.; data curation, S.R.F.; writing—original draft preparation, I.J.Y.S. and S.R.F.; writing—review and editing, S.R.F., A.-J.C.-L., W.E.S., A.K.D.P. and E.M.-T.; visualization, S.R.F. and I.J.Y.S.; supervision, S.R.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The snow data are available from the Water and Climate Center of the Natural Resources Conservation Service at https://www.wcc.nrcs.usda.gov/ (date accessed on 16 August 2021). The elevation data are the U.S. Geological Survey National Map at https://nationalmap.gov/ (date accessed on 16 August 2021). Data were taken from Fassnacht and Records [11].

Acknowledgments

This work came about after a discussion with Luca Egli. The authors appreciate his insight. One anonymous reviewer provided a thorough review of the paper that helped correct various grammatical errors.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Regression Methodology

Three best-fit regression methods were fit to the accumulation data: linear, truncated linear (excluding mean SWE < 35 mm), and power for seven snow years with unique accumulation phase characteristics. The 11-day moving-average breakpoint method identified in the Methods section was used to identify the end of the accumulation phase. The truncated linear model yielded the largest coefficient of determination values in 6 of 7 years (Table A1).
Table A1. Regression parameters and coefficients of determination for three regression methods based on data from seven characteristic snow years.
Table A1. Regression parameters and coefficients of determination for three regression methods based on data from seven characteristic snow years.
Snow YearLinear (Equation (1))Linear Excluding
SWE < 35 mm
Power Function
(Equation (2))
R2bmR2bmR2αβ
19930.9983.480.4320.9976.580.4250.9830.8090.715
19950.9809.260.4590.99123.430.3930.9590.6800.753
20020.9867.260.4310.9899.130.4190.9861.2600.608
20100.97716.520.4050.99428.010.3640.9841.4120.614
20110.98711.230.4950.99729.900.4350.9951.0860.702
20120.9808.360.3590.99714.790.3260.9711.2290.592
20130.9877.560.4500.99516.640.4060.9911.2570.623

Appendix B. Inflection Point Identification

A variety of approaches were tested to determine a replicable method to identify the inflection point between the purely accumulation and mixed accumulation/ablation phases, independent of the regression model applied to the accumulation phase. Three time periods were tested to determine the inflection point: leading, lagging and centralized daily linear regression methods were tested, with date ranges varying between 10 (leading/lagging) and 11 days (central). Preliminary attempts to implement this method resulted in many false positives—changes in slope above the threshold, but much earlier than the measured accumulation phase terminated. To mitigate false positives, a minimum date threshold of 15 February was included to disallow pre-emptive terminations of the accumulation phase.
After a series of preliminary analyses, two slope difference thresholds were identified for larger-scale analyses. The two slope difference thresholds (0.3 and 0.5) were tested with leading 10-day, lagging 10-day, and central 11-day methods for seven snow years. The seven years selected for testing exhibited unique characteristics on the plot of standard deviation versus mean SWE. Results were compared to the observed inflection point, which was determined through visual inspection. Based on the observed inflection point, the Nash–Sutcliffe model coefficient of efficiency was calculated to assess how accurately the predicted inflection point aligned with the observed inflection point (Table A2 and Table A3).
Table A2. Comparison of breakpoint prediction values of mean SWE for slope difference threshold of 0.5.
Table A2. Comparison of breakpoint prediction values of mean SWE for slope difference threshold of 0.5.
Snow YearObserved Mean SWE [mm]Modeled
Leading 10-DayCentral 11-DayLagging 10-Day
1993629235628526
1995411348347342
2002262259261256
2010343001
2011445455454445
2012317313316316
2013295294302297
NSE0.7060.8470.810
R20.7570.8550.870
Table A3. Comparison of breakpoint prediction values of mean SWE for slope difference threshold of 0.3.
Table A3. Comparison of breakpoint prediction values of mean SWE for slope difference threshold of 0.3.
Snow YearObserved Mean SWE [mm]Modeled
Leading 10-DayCentral 11-DayLagging 10-Day
1993629527624526
1995411348347342
2002262174261253
2010343414444391
2011445455454445
2012317313316316
2013295296303296
NSE0.7060.8470.810
R20.7570.8550.870

Appendix C. Evaluation of Accumulation Slopes

The process of identifying inflection points, and executing a linear regression for each trajectory (34 years of data, 7 sub-sets) was simplified with development of a number of Python scripts. However, the method and parameters were not able to identify a point of inflection for 1 one instance out of 238 (0.42%). For this trajectory (north-low sub-set, 1987, Figure 1), the inflection point was easily identified visually, but was not detected numerically.
After reviewing the raw data obtained through Python scripts, nine specific instances with false inflection points were discovered. False inflection points were defined as occurring after day 250 of the water year (i.e., June, well after melt phase should have begun), and often were characterized by a coefficient of determination substantially lower than 0.7. This occurred throughout all sub-sets, over a variety of years (1 in the summary sub-set (2000), 2 in the south sub-set (2000, 2012), 2 in the north-high sub-set (1984, 2001), 2 in the south-high sub-set (1987, 2000), 2 in the south-low sub-set (2000, 2008)). The local slope near the visually observed breakpoints ranged from 0.286 to 0.465 (0.391, 0.384, 0.292, 0.322, 0.371, 0.286, 0.391, 0.465, 0.402, respectively). The issue appeared to be more prevalent in non-average snow years (5 in low, 3 in high and 1 in mid).

References

  1. Kearney, M.S.; Harris, B.H.; Hershbein, B.; Jácome, E.; Nantz, G. In times of Drought: Nine Economic Facts about Water in the United States. The Brookings Institution Policy Memo. Hamilt. Proj. 2014. Available online: https://www.brookings.edu/research/in-times-of-drought-nine-economic-facts-about-water-in-the-united-states/ (accessed on 16 August 2021).
  2. He, X.; Pan, M.; Wei, Z.; Wood, E.F.; Sheffield, J. A Global Drought and Flood Catalogue from 1950 to 2016. Bull. Am. Meteorol. Soc. 2020, 101, E508–E535. [Google Scholar] [CrossRef] [Green Version]
  3. Güneralp, B.; Güneralp, I.; Liu, Y. Changing global patterns of urban exposure to flood and drought hazards. Glob. Environ. Chang. 2015, 31, 217–225. [Google Scholar] [CrossRef]
  4. Hammond, J.C.; Saavedra, F.A.; Kampf, S.K. Global snow zone maps and trends in snow persistence 2001–2016. Int. J. Clim. 2018, 38, 4369–4383. [Google Scholar] [CrossRef]
  5. Immerzeel, W.W.; Lutz, A.F.; Andrade, M.; Bahl, A.; Biemans, H.; Bolch, T.; Hyde, S.; Brumby, S.; Davies, B.; Elmore, A.C.; et al. Importance and vulnerability of the world’s water towers. Nat. Cell Biol. 2019, 577, 364–369. [Google Scholar] [CrossRef] [PubMed]
  6. Maupin, M.A.; Kenny, J.F.; Hutson, S.S.; Lovelace, J.K.; Barber, N.L.; Linsey, K.S. Estimated use of water in the United States in 2010. U.S. Geol. Surv. Circ. 2014. [Google Scholar] [CrossRef]
  7. Piechota, T.C.; Timilsena, J.; Tootle, G.; Hidalgo, H. The western U.S. drought: How bad is it? Eos Trans. Am. Geophys. Union 2004, 85, 301–304. [Google Scholar] [CrossRef]
  8. Williams, A.P.; Cook, E.R.; Smerdon, J.E.; Cook, B.I.; Abatzoglou, J.T.; Bolles, K.; Baek, S.H.; Badger, A.M.; Livneh, B. Large contribution from anthropogenic warming to an emerging North American megadrought. Science 2020, 368, 314–318. [Google Scholar] [CrossRef]
  9. Clow, D. Changes in the Timing of Snowmelt and Streamflow in Colorado: A Response to Recent Warming. J. Clim. 2010, 23, 2293–2306. [Google Scholar] [CrossRef]
  10. Fassnacht, S.R.; Venable, N.B.; McGrath, D.; Patterson, G.G. Sub-Seasonal Snowpack Trends in the Rocky Mountain National Park Area, Colorado, USA. Water 2018, 10, 562. [Google Scholar] [CrossRef] [Green Version]
  11. Fassnacht, S.; Patterson, G.; Venable, N.; Cherry, M.; Pfohl, A.; Sanow, J.; Tedesche, M. How Do We Define Climate Change? Considering the Temporal Resolution of Niveo-Meteorological Data. Hydrology 2020, 7, 38. [Google Scholar] [CrossRef]
  12. Thackeray, C.W.; Derksen, C.; Fletcher, C.G.; Hall, A. Snow and Climate: Feedbacks, Drivers, and Indices of Change. Curr. Clim. Chang. Rep. 2019, 5, 322–333. [Google Scholar] [CrossRef]
  13. Cayan, D.R.; Das, T.; Pierce, D.W.; Barnett, T.P.; Tyree, M.; Gershunov, A. Future dryness in the southwest US and the hydrology of the early 21st century drought. Proc. Natl. Acad. Sci. USA 2010, 107, 21271–21276. [Google Scholar] [CrossRef] [Green Version]
  14. Apurv, T.; Cai, X. Regional Drought Risk in the Contiguous United States. Geophys. Res. Lett. 2021, 48. [Google Scholar] [CrossRef]
  15. Bales, R.C.; Molotch, N.P.; Painter, T.H.; Dettinger, M.D.; Rice, R.; Dozier, J. Mountain hydrology of the western United States. Water Resour. Res. 2006, 42. [Google Scholar] [CrossRef]
  16. Blöschl, G. Scaling issues in snow hydrology. Hydrol. Process. 1999, 13, 2149–2175. [Google Scholar] [CrossRef]
  17. Sturm, M.; Wagner, A.M. Using repeated patterns in snow distribution modeling: An Arctic example. Water Resour. Res. 2010, 46. [Google Scholar] [CrossRef] [Green Version]
  18. Egli, L.; Jonas, T. Hysteretic dynamics of seasonal snow depth distribution in the Swiss Alps. Geophys. Res. Lett. 2009, 36. [Google Scholar] [CrossRef] [Green Version]
  19. Meiman, J.R. Snow accumulation related to elevation, aspect, and forest canopy. In Proceedings of the Snow Hydrology, Snow Hydrology Workshop Seminar, Fredericton, NB, Canada, 28–29 February 1968; Canadian National Committee of the International Hydrological Decade: Fredericton, NB, Canada, 1968; pp. 35–47. Available online: https://www.cabdirect.org/cabdirect/abstract/19730606963 (accessed on 16 August 2021).
  20. Elder, K.; Dozier, J.; Michaelsen, J. Snow accumulation and distribution in an Alpine Watershed. Water Resour. Res. 1991, 27, 1541–1552. [Google Scholar] [CrossRef] [Green Version]
  21. Erxleben, J.; Elder, K.; Davis, R. Comparison of spatial interpolation methods for estimating snow distribution in the Colorado Rocky Mountains. Hydrol. Process. 2002, 16, 3627–3649. [Google Scholar] [CrossRef]
  22. Meromy, L.; Molotch, N.P.; Link, T.E.; Fassnacht, S.; Rice, R.H. Subgrid variability of snow water equivalent at operational snow stations in the western USA. Hydrol. Process. 2012, 27, 2383–2400. [Google Scholar] [CrossRef]
  23. McCreight, J.L.; Slater, A.; Marshall, H.P.; Rajagopalan, B. Inference and uncertainty of snow depth spatial distribution at the kilometre scale in the Colorado Rocky Mountains: The effects of sample size, random sampling, predictor quality, and validation procedures. Hydrol. Process. 2012, 28, 933–957. [Google Scholar] [CrossRef]
  24. Fassnacht, S.R.; Brown, K.S.J.; Blumberg, E.J.; Lopez-Moreno, I.; Covino, T.P.; Kappas, M.; Huang, Y.; Leone, V.; Kashipazha, A.H. Distribution of snow depth variability. Front. Earth Sci. 2018, 12, 683–692. [Google Scholar] [CrossRef]
  25. Sexstone, G.A.; Fassnacht, S.R. What drives basin scale spatial variability of snowpack properties in northern Colorado? Cryosphere 2014, 8, 329–344. [Google Scholar] [CrossRef] [Green Version]
  26. Fassnacht, S.R.; Dressler, K.A.; Bales, R.C. Snow water equivalent interpolation for the Colorado River Basin from snow telemetry (SNOTEL) data. Water Resour. Res. 2003, 39. [Google Scholar] [CrossRef] [Green Version]
  27. Dozier, J.; Bair, E.H.; Davis, R.E. Estimating the spatial distribution of snow water equivalent in the world’s mountains. Wiley Interdiscip. Rev. Water 2016, 3, 461–474. [Google Scholar] [CrossRef]
  28. Barabási, A.L.; Stanley, H.E. Fractal Concepts in Surface Growth; Cambridge University Press: Cambridge, UK, 1995. [Google Scholar]
  29. Crow, W.T.; Wood, E.F. Multi-scale dynamics of soil moisture variability observed during SGP’97. Geophys. Res. Lett. 1999, 26, 3485–3488. [Google Scholar] [CrossRef]
  30. Famiglietti, J.; Ryu, D.; Berg, A.; Rodell, M.; Jackson, T.J. Field observations of soil moisture variability across scales. Water Resour. Res. 2008, 44. [Google Scholar] [CrossRef] [Green Version]
  31. Fassnacht, S.R.; Derry, J.E. Defining similar regions of snow in the Colorado River Basin using self-organizing maps. Water Resour. Res. 2010, 46. [Google Scholar] [CrossRef]
  32. Fassnacht, S.R.; Records, R.M. Large snowmelt versus rainfall events in the mountains. J. Geophys. Res. Atmos. 2015, 120, 2375–2381. [Google Scholar] [CrossRef]
  33. Fassnacht, S.R.; Sexstone, G.A.; Kashipazha, A.H.; Lopez-Moreno, I.; Jasinski, M.F.; Kampf, S.K.; Von Thaden, B.C. Deriving snow-cover depletion curves for different spatial scales from remote sensing and snow telemetry data. Hydrol. Process. 2015, 30, 1708–1717. [Google Scholar] [CrossRef]
  34. Fassnacht, S.R.; Lopez-Moreno, I.; Ma, C.; Weber, A.N.; Pfohl, A.K.D.; Kampf, S.K.; Kappas, M. Spatio-temporal snowmelt variability across the headwaters of the Southern Rocky Mountains. Front. Earth Sci. 2017, 11, 505–514. [Google Scholar] [CrossRef]
  35. Ma, C.; Fassnacht, S.; Kampf, S. How Temperature Sensor Change Affects Warming Trends and Modeling: An Evaluation across the State of Colorado. Water Resour. Res. 2019, 55, 9748–9764. [Google Scholar] [CrossRef]
  36. von Thaden, B.C. Spatial Accumulation Patterns of Snow Water Equivalent in the Southern Rocky Mountains. Master’s Thesis, Colorado State University, Fort Collins, CO, USA, 2016. [Google Scholar]
  37. Löwe, H.; Egli, L.; Bartlett, S.; Guala, M.; Manes, C. On the evolution of the snow surface during snowfall. Geophys. Res. Lett. 2007, 34. [Google Scholar] [CrossRef] [Green Version]
  38. Pomeroy, J.; Essery, R.; Toth, B. Implications of spatial distributions of snow mass and melt rate for snow-cover depletion: Observations in a subarctic mountain catchment. Ann. Glaciol. 2004, 38, 195–201. [Google Scholar] [CrossRef] [Green Version]
  39. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. McKinney, W. Data Structures for Statistical Computing in Python. In Proceedings of the 9th Python in Science Conference, Austin, TX, USA, 28 June–3 July 2010; pp. 51–56. [Google Scholar]
  41. Fassnacht, S.R.; Patterson, G.G. Niveograph Interpolation to Estimate Peak Accumulation at Two Mountain Sites. Cold and Mountain Region Hydrological Systems under Climate Change: Towards Improved Projections. In Proceedings of the Symposium H02, IAHS-IAPSO-IASPEI Assembly, Gothenburg, Sweden, 22–26 July 2013; IAHS: Oxfordshire, UK, 2013; Volume 360, pp. 59–64. [Google Scholar]
  42. Räisänen, J. Warmer climate: Less or more snow? Clim. Dyn. 2007, 30, 307–319. [Google Scholar] [CrossRef]
  43. Guest, C.; Shearer, H.; McKean, M.; Reiner, R. America. Track 5 on This Is Spinal Tap; Polydor: London, UK, 1984. [Google Scholar]
  44. Mudryk, L.R.; Kushner, P.J.; Derksen, C. Interpreting observed northern hemisphere snow trends with large ensembles of climate simulations. Clim. Dyn. 2013, 43, 345–359. [Google Scholar] [CrossRef]
  45. Steffen, W.; Richardson, K.; Rockström, J.; Schellnhuber, H.J.; Dube, O.P.; Dutreuil, S.; Lenton, T.M.; Lubchenco, J. The emergence and evolution of Earth System Science. Nat. Rev. Earth Environ. 2020, 1, 54–63. [Google Scholar] [CrossRef] [Green Version]
  46. Montanari, A.; Young, G.; Savenije, H.; Hughes, D.; Wagener, T.; Ren, L.; Koutsoyiannis, D.; Cudennec, C.; Toth, E.; Grimaldi, S.; et al. “Panta Rhei—Everything Flows”: Change in hydrology and society—The IAHS Scientific Decade 2013–2022. Hydrol. Sci. J. 2013, 58, 1256–1275. [Google Scholar] [CrossRef]
  47. Steffen, W.; Rockström, J.; Richardson, K.; Lenton, T.M.; Folke, C.; Liverman, D.; Summerhayes, C.P.; Barnosky, A.D.; Cornell, S.E.; Crucifix, M.; et al. Trajectories of the Earth System in the Anthropocene. Proc. Natl. Acad. Sci. USA 2018, 115, 8252–8259. [Google Scholar] [CrossRef] [Green Version]
  48. Kotzé, L.J.; Kim, R.E. Earth system law: The juridical dimensions of earth system governance. Earth Syst. Gov. 2019, 1, 100003. [Google Scholar] [CrossRef]
Figure 1. Location of the stations (a) as the study map illustrating the elevation group (high is 2800 m in the north and low is 2900 m in the south), land cover (evergreen/non-evergreen) and mean peak SWE (mm) of 90 SNOTEL stations across the Southern Rocky Mountains. The latitudinal division is noted by the dashed line. The study area is (b) located in the Western United States. The number of stations in each group (Section 3.1) is presented in (a).
Figure 1. Location of the stations (a) as the study map illustrating the elevation group (high is 2800 m in the north and low is 2900 m in the south), land cover (evergreen/non-evergreen) and mean peak SWE (mm) of 90 SNOTEL stations across the Southern Rocky Mountains. The latitudinal division is noted by the dashed line. The study area is (b) located in the Western United States. The number of stations in each group (Section 3.1) is presented in (a).
Hydrology 08 00124 g001
Figure 2. Mean annual peak SWE across the entire domain highlighting the above average (blue), average (green), and below average (red) snow years. The seven characteristic snow years that were used to assess the regression model are identified by the bold outline and peak SWE is labeled (above average: 1993, 1995, 2011; average: 2010; below average: 2002, 2012, 2013).
Figure 2. Mean annual peak SWE across the entire domain highlighting the above average (blue), average (green), and below average (red) snow years. The seven characteristic snow years that were used to assess the regression model are identified by the bold outline and peak SWE is labeled (above average: 1993, 1995, 2011; average: 2010; below average: 2002, 2012, 2013).
Hydrology 08 00124 g002
Figure 3. Standard deviation versus mean SWE trajectory across the full domain for an average snow year 2010, high snow year 2011, and low snow year 2012, with accumulation start (non-linear) as a dashed arrow, and the accumulation/melt inflection points identified. The approximate accumulation phase at all stations (square symbol), accumulation at some stations & melt at others (triangles symbol) and melt at all stations (round symbol). The three solid arrows show evolution over time, illustrating the direction of accumulation, accumulation & melt, and melt periods.
Figure 3. Standard deviation versus mean SWE trajectory across the full domain for an average snow year 2010, high snow year 2011, and low snow year 2012, with accumulation start (non-linear) as a dashed arrow, and the accumulation/melt inflection points identified. The approximate accumulation phase at all stations (square symbol), accumulation at some stations & melt at others (triangles symbol) and melt at all stations (round symbol). The three solid arrows show evolution over time, illustrating the direction of accumulation, accumulation & melt, and melt periods.
Hydrology 08 00124 g003
Figure 4. Mean and standard deviation of mean daily data for snow depth (87 SNOTEL stations) and SWE (90 SNOTEL stations) for (a) 2010 (average snow year), (b) 2011 (above average snow year) and (c) 2012 (below average snow year). The slope and R2 values listed in the lower right of each figure are for the accumulation phase.
Figure 4. Mean and standard deviation of mean daily data for snow depth (87 SNOTEL stations) and SWE (90 SNOTEL stations) for (a) 2010 (average snow year), (b) 2011 (above average snow year) and (c) 2012 (below average snow year). The slope and R2 values listed in the lower right of each figure are for the accumulation phase.
Hydrology 08 00124 g004
Figure 5. Mean accumulation slopes per year from the standard deviation versus mean trajectories for (a) the entire domain, (b) the north station sub-sets, and (c) the south station sub-sets. The station subsets are further divided into high and low elevation.
Figure 5. Mean accumulation slopes per year from the standard deviation versus mean trajectories for (a) the entire domain, (b) the north station sub-sets, and (c) the south station sub-sets. The station subsets are further divided into high and low elevation.
Hydrology 08 00124 g005
Table 1. Comparison of coefficient of determination values from three accumulation phase regression models for the seven study years (full domain) (see Figure 4). The type of snow year (Figure 2) is listed as above average (AA), average (AVG) or below average (BA).
Table 1. Comparison of coefficient of determination values from three accumulation phase regression models for the seven study years (full domain) (see Figure 4). The type of snow year (Figure 2) is listed as above average (AA), average (AVG) or below average (BA).
Snow YearLinearLinear Truncated at 35 mmPower
1993 (AA)0.9980.9970.983
1995 (AA)0.9800.9910.959
2002 (BA)0.9860.9890.986
2010 (AVG)0.9770.9940.984
2011 (AA)0.9870.9970.995
2012 (BA)0.9800.9970.971
2013 (BA)0.9870.9950.991
average0.9850.9940.981
Table 2. Accumulation phase inflection point model parameter comparison.
Table 2. Accumulation phase inflection point model parameter comparison.
Slope Difference ThresholdIndicator PositionNSER2
0.3leading 10-day0.710.76
central 10-day0.850.86
lagging 10-day0.810.87
0.5leading 10-day−1.970.01
central 10-day−0.310.55
lagging 10-day−0.420.44
Table 3. Full domain and six sub-sets summary of accumulation phase slope statistics (maximum, mean, minimum, standard deviation) plus the mean R2 of all truncated linear fits over the entire 34-year period of analysis, and mean slope for above average (AA), average (AVG), and below average (BA) snow years. See Figure 5 for specific years.
Table 3. Full domain and six sub-sets summary of accumulation phase slope statistics (maximum, mean, minimum, standard deviation) plus the mean R2 of all truncated linear fits over the entire 34-year period of analysis, and mean slope for above average (AA), average (AVG), and below average (BA) snow years. See Figure 5 for specific years.
GroupSub-SetSlopeMean R2AAAVGBA
Max.MeanMin.Std. Dev.
summaryfull domain0.590.400.260.0620.980.370.390.33
latitude-basedNorth0.470.390.290.0420.990.360.370.33
South0.500.370.230.0770.970.370.370.33
latitude and elevation-basednorth-high0.490.390.280.0530.980.340.370.34
north-low0.470.380.260.0490.990.360.360.28
south-high0.570.370.230.0830.970.390.340.29
south-low0.520.360.120.0970.950.350.360.25
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Schrock, I.J.Y.; Fassnacht, S.R.; Collados-Lara, A.-J.; Sanford, W.E.; Pfohl, A.K.D.; Morán-Tejeda, E. Snow Water Equivalent Accumulation Patterns from a Trajectory Approach over the U.S. Southern Rocky Mountains. Hydrology 2021, 8, 124. https://doi.org/10.3390/hydrology8030124

AMA Style

Schrock IJY, Fassnacht SR, Collados-Lara A-J, Sanford WE, Pfohl AKD, Morán-Tejeda E. Snow Water Equivalent Accumulation Patterns from a Trajectory Approach over the U.S. Southern Rocky Mountains. Hydrology. 2021; 8(3):124. https://doi.org/10.3390/hydrology8030124

Chicago/Turabian Style

Schrock, Isaac J. Y., Steven R. Fassnacht, Antonio-Juan Collados-Lara, William E. Sanford, Anna K. D. Pfohl, and Enrique Morán-Tejeda. 2021. "Snow Water Equivalent Accumulation Patterns from a Trajectory Approach over the U.S. Southern Rocky Mountains" Hydrology 8, no. 3: 124. https://doi.org/10.3390/hydrology8030124

APA Style

Schrock, I. J. Y., Fassnacht, S. R., Collados-Lara, A. -J., Sanford, W. E., Pfohl, A. K. D., & Morán-Tejeda, E. (2021). Snow Water Equivalent Accumulation Patterns from a Trajectory Approach over the U.S. Southern Rocky Mountains. Hydrology, 8(3), 124. https://doi.org/10.3390/hydrology8030124

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop