Can Trends in Snowmelt Streamflow Be Explained by Terrain Characteristics or Trends in Vegetation, Precipitation and SWE? An Examination of Small High-Elevation Colorado (U.S.A.) Rivers

1. Introduction

Mountain snowmelt generates water for streams and rivers and is a major source for a substantially increasing portion of the Earth’s population [2]. Across the semi-arid western United States, a majority of the precipitation falls as snow [3,4,5]. The timing of the start of the melt season and when the snowmelt enters streams in the high elevation watersheds is crucial for estimating water availability [1,6]; this timing has shifted [7,8,9,10,11] due mostly to climate change [12,13,14,15,16].

1.1. Snowmelt Streamflow Timing Metrics

Peak flow date is a simple metric of streamflow timing but neglects the remaining data for a given year [17]. Court [17] introduced the half-flow or Center of Volume date, i.e., the day when 50% of total annual flow has passed a stream gauging station (t_Q50), to assess the characteristics of streamflow timing. This t_Q50 is used extensively as a streamflow timing metric [18,19,20], especially to assess the impacts of climate change [1,13,15,21,22,23]. Other percentages of annual flow passage have been used as proxies for the start of snowmelt contribution, i.e., the date of 20% (t_Q20) [15] or 25% [22] of flow, and the end of snowmelt, i.e., the date of 75% [22] or 80% (t_Q80) [15] of flow. To highlight the snowmelt period further, Dudley et al. [23] proposed the Center of Volume (COV) for 50% of the flow from January to July (t_QDudley). However, all these methods are static based on a specific quantity of the total annual (or winter [23]) streamflow. The use of t_Q20, t_Q50, and t_Q80 are not appropriate indicators specifically of snowmelt timing [24] and sometimes due to large precipitation events [20,24]. These metrics can be influenced by inter-annual variability in streamflow volume [24] and do not reflect how a changing climate impacts streamflow timing [1].

1.2. Objectives of the Paper

An analytical approach considering the change in streamflow using a departure from baseflow has previously been proposed to identify the start (t_Qstart) and end of the snowmelt contribution to streamflow (t_Qend), i.e., considering the characteristics of the hydrograph [1]. Here, we use that approach to quantify the start and end of snowmelt from the hydrograph and to determine if and how snowmelt timing is changing. Since snowmelt streamflow characteristics change for a variety of reasons [25,26], we evaluate these changes considering the terrain parameters (e.g., elevation, aspect, slope), canopy from the Normalized Difference Vegetation Index (NDVI), and winter precipitation from the Parameter-elevation Regression on Independent Slopes Model (PRISM) dataset for 39 watersheds, less than 900 km² in size. Temperature has had a strong correlation to trends in the specified percentage of flow that has passed [15,23]. However, there is a significant inhomogeneity in the middle (approximately 1998 to 2007) of the time series at the high elevation Snow Telemetry (SNOTEL) stations used to derived mountain temperatures across the Western U.S. [27,28,29]. Therefore, spatial temperature data were not used to investigate snowmelt streamflow changes in Colorado mountain watersheds over the 40-year study period. We used the SNOTEL Snow Water Equivalent (SWE) station closest to each streamflow gauge to identify annual peak SWE.

The objectives of this paper are as follows: 1) apply a new method of estimating snowmelt timing and volume for streamflow for the Southern Rocky Mountains of Colorado, 2) conduct a trend analysis of different snowmelt timing and volume variables, 3) determine possible explanations for these trends based on time trends in vegetation, winter precipitation, and peak SWE, as well as terrain parameters. We explored these high-elevation watersheds in Colorado, as the state of Colorado is a headwater state. These include the Colorado River and its tributaries (Yampa, Gunnison, Uncompahgre, San Miguel, Dolores, Animas, and San Juan) [20], the North and South Platte Rivers, the Arkansas River, and the Rio Grande [15]. The highest mountain peaks reach over 4,400 meters in elevation and snowcover persists from October through May [3]. Streamflow in these watersheds is snowmelt dominated, with 60 to 80% of the annual streamflow coming from snow [4,6,15,16].

2. Methodology

The t_Qstart, or timing (date) of the start of the snowmelt contribution to streamflow, was computed as the increase in streamflow from baseflow by a change in slope of at least 10 mm/day [1]. The t_Qend, or the date of the end of snowmelt contribution, was computed as the decrease in streamflow back to baseflow as the change in slope of at least 17.5 mm/day [1]. The timing of snowmelt (t_Qstart-end) was computed as the number of days between t_Qstart and t_Qend (Figure A1). This was used in lieu of t_Q50 or t_QDudley. We then determined volumes of flow, in particular the total annual runoff (Q₁₀₀), the volume that passed the gauge prior to the start and end of snowmelt (Q_start and Q_end, respectively), and the volume in between (Q_start-end) (Figure A1).

The rate of change for the trends were calculated as the Theil-Sen’s Slope [31,32], and the significance was calculated using the Mann-Kendall Test [33,34]. Since previous studies that have examined trends in timing of streamflow snowmelt have primarily relied on climactic indices to explain their observations [13,15,16,23], we used precipitation data from the Precipitation-Elevation Regressions on Independent Slopes Model (PRISM) dataset [35] to evaluate winter precipitation (October through March), starting in 1982.

We included mean incoming winter solar radiation, basin elevation, basin slope, and location (latitude and longitude) [36] to evaluate parameters that could influence trends. To address changes within the watershed from land use, beetle-kill, or wildfires, we collected NDVI data from the U.S. Geological Survey [37]. These data start in July of 1989 but would still capture major changes in vegetation because major fires and beetle-kill didn’t occur in the Southern Rocky Mountains until the late 1990s and early 2000s [38,39]. We calculated the correlation coefficients between the trends in snowmelt streamflow timing and flow volume versus terrain parameters, plus trends in vegetation and precipitation. Further, we evaluated multi-variate linear regressions using all the variables and the most highly correlated variables from the individual regression. The independent variables were standardized to between 0 and 1 so that the coefficients could be compared for each regression.

3. Study Domain

We examined 40-years of streamflow (1976 through 2015) for 39 United States Geological Survey (USGS) gauging stations across the Southern Rocky Mountains of Colorado, each with at least 30 years of record (Figure 1; Table A1). Streamflow data were obtained from the National Water Information System [30]. All were headwater streams gauged at an elevation higher than 2000 meters above sea level (Figure 1; Table A1). The mean basin elevation varied from 2494 to 3644 m.a.s.l. (Figure 2a), with the mean April clear sky solar radiation loading of 1407 to 1760 Wh/m² (Figure 2b). The basins had a mean slope from 17 to 26^o (Figure 2c) and ranged in size from 4 to 878 km² in size (Figure 2d). The stations are summarized in Pfohl and Fassnacht [1]. The SWE data were obtained from the Natural Resources Conservation Service [40].

4. Results

The canopy density, as per NDVI was increasing for 37 of the 39 all watersheds (Figure 2e), significantly at four (moderately significant at one). Both winter precipitation (Figure 2f) and the adjacent SNOTEL peak SWE (Figure 2g) were decreasing for all but one watershed. Most of the trends in winter precipitation and peak SWE were not significant.

The snowmelt characteristics of streamflow have changed across most of the watersheds over the study period (Figure 3). Most (34) see a trend of an earlier start of the snowmelt streamflow, while only three are later (Figure 3a), with a third being of the trends being significant (and five being moderately significant). Twenty-nine watersheds see an earlier end of the snowmelt contribution and seven are later (Figure 3b), with about 40% being significant (9 watersheds) or moderately significant (7 watersheds). The change in timing of the peak, denoted t_Qstart-end, is mixed, being earlier at 12 watersheds and later at 20 (1 significantly in each direction; Figure 3c). For 27 watersheds, both t_Qstart and t_Qend trends were earlier while for only one watershed both became later. Earlier trends were observed for all three metrics in nine watersheds.

Trends for the volume of flow that has passed the gauge were more mixed (Figure 3d–g), i.e., both increasing and decreasing. Total annual streamflow (Q₁₀₀) increased in 23 (1 significantly and 3 moderately significant) watersheds while it decreased at the (16) others (2 significantly and 1 moderately significant) (Figure 3d). Q_start changed by the smallest amount (Figure 3e). Trends for Q_end (Figure 3f) and Q_start-end (Figure 3g) were similar (16 with more and 23 with less) with 35 having the same sign (15 watersheds where both increased streamflow and 20 where both decreased). The trend was in the same direction for Q_start and Q_end at 20 watersheds (7 less, 13 more), and for 18 watersheds for all four metrics (6 less, 12 more). Trends were in the same direction and significant for three watersheds: Joe Wright Creek (more streamflow), Vasquez Creek (more streamflow), and Conejos River (less streamflow). The trends in Q₁₀₀, Q_end and Q_start-end illustrated a latitudinal pattern with most stations north of 39.7^o increasing in flow volume and most south decreasing (Figure 3d,f,g).

With the exception of Q_start, snowmelt streamflow trends are more correlated to winter precipitation or peak SWE than NDVI (Figure 4). Winter precipitation is more correlated with NDVI (R = 0.43) than with peak SWE (R = 0.29). Mean radiation and elevation are weakly correlated to streamflow. Mean basin slope is significantly correlated (negatively) to trends in t_Qend, Q₁₀₀, Q_end and Q_start-end. Latitude is positively correlated to all streamflow characteristics (3 significantly) while longitude is less correlated (except t_Qend and t_Qstart-end that are moderately significant) than latitude.

A linear multi-variate regression illustrates some significant correlation between trends in streamflow timing and volume metrics with terrain parameters (Table 1), but not with variables with trends, i.e., NDVI, winter precipitation, or peak SWE. The strongest correlations were for Q_end and Q_start-end including all variables (R² of 0.52 and 0.46, respectively). The regression between Q₁₀₀ and all variables was moderately significant (p < 0.1). Considering only some of the variables made the regressions for t_Qend, Q₁₀₀, Q_end, and Q_start-end significant (Table 1b–d). The variance explained decreased, as shown by R², but the individual regression variables were more consistently significant. Slope was negatively correlated, and latitude was positively correlated with snowmelt streamflow timing and volume trends. Trends in the start of snowmelt streamflow, i.e., t_Qstart and Q_start, as well as t_Qstart-end, were poorly explained by the regression, not significant, and R² was mostly less than 0.1 (Table 1).

5. Discussion

Snowmelt-driven streamflow is occurring earlier for most basins across the study domain (t_Qstart in Figure 2e and the t_Qend in Figure 2f), as has also been seen using time-constant streamflow metrics, i.e., t_Q20, and t_Q80 [7,8,9,10,11,13,15,16,21,22]. However, the trends in the t_Qstart-end or mean of start and end timing were mixed (Figure 3c, Figure 3, Table 1), reflecting t_Qstart and t_Qend trends (Table A2) and their differences (Figure 3a versus Figure 3b). The time-constant metrics that are meant to represent the middle of the snowmelt streamflow peak, i.e., t_Q50 [17] or t_QDudley [23], are getting earlier [15]. These time-constant streamflow metrics are relative to the water year, while t_Qstart-end is relative to the characteristics of the hydrograph [1], as recommended by Whitfield [24]. The change in slope for identifying the start and end of melt (here 10 and 17.5 mm/day) can influence the estimation of the timing metrics, possibly for other climate regions. For high-elevation Colorado streamflow gauges, the values used herein were shown to be acceptable minima [1].

The title of this paper and the third objective of this study is to determine if trends in snowmelt streamflow can be explained by watershed parameters or trends in canopy, precipitation and peak SWE (temperature was not assessed as described above). The simple answer is that slope, negatively, and latitude, positively, explained changes in total flow (Q₁₀₀), the end of snowmelt streamflow (t_Qend, Q_end), and Q_start-end (Figure 4 and Table 1). The negative correlation with mean basin slope could imply that gentler sloped watersheds are possibly melting out later and thus having a later t_Qend and larger Q_end [41]. Higher elevation watersheds tend to be steeper (Table A2). However, slopes usually vary substantially across mountain watersheds and the mean slope may not represent watershed processes well [42]. Winter precipitation across the state of Colorado is correlated with latitude (R = 0.59 in Table A3) [43], with southern stations seeing a larger in snowfall (Figure 2f) since about 2000 [39,44]. This correlation is also seen between NDVI and winter precipitation (R = 0.43 in Table A3) and thus NDVI and latitude (R = 0.40). Peak SWE trends were correlated with elevation (Table A3) [45,46], but here (Figure 2g), less correlated with winter precipitation (R = 0.29 in Table A3). Peak SWE was extracted from SNOTEL station data [40] and these may not be representative of the watershed [47]. These are mostly small watersheds (Figure 2d), so current SWE products [48] may not have the necessary resolution to assess changes. Snowpack and hydrological modeling could provide more insight into changes of processes that may dictate altering of streamflow timing [49].

There are some spatial patterns in the changes in snowmelt-driven streamflow, specifically latitude, and to a lesser degree longitude (Figure 3 and Figure 4). Others [15] have used the Regional Kendall test [50] to evaluate trends and their significance across an area; due to the limited spatial patterns observed here (Figure 3), it is recommended that Mann-Kendall test [33,34] and Theil-Sen slope [31,32] on individual stations. Using the Regional Kendall test can produce trends that are smaller in magnitude than observed trends at individual sites [51].

The method used herein presents the timing and volume of water at the start, end, and average of the peak (start-end) from snowmelt contribution [1]. This information may be helpful for water forecasters and managers making decisions about water storage and reservoirs for the future [52], especially if timing of peak flow is incorporated [53,54].

The approach used herein [1] identifies the start and end of the snowmelt contribution for snowmelt dominated systems as an improvement to the traditional statics approaches, such as t_Q20, t_Q50, and t_Q80 [15,17]. It does not specifically identify baseflow, although it has been used for that purpose [55]. Baseflow separation techniques could be used to identify when direct or non-baseflow starts to contribute to streamflow. This could be applied to a snowmelt dominated system to determine when snowmelt streamflow started and ended. There are analytical approaches [54] using only streamflow data. Snowmelt is often separated from baseflow using isotopes [55]. However, such measurements are labor and cost intensive. Specific conductance is measured as an in-situ water quality variable in a few locations and has been used with streamflow to separate baseflow from non-baseflow [56]. There are now some time series long enough to examine trends.

Temperature increases are a major indication and result of climate change [12,13,14,15,16]. Where temperature data are reliable, these data can be used to assess changes to snowmelt-driven streamflow. Across Colorado [29] and the western U.S. [27,28], the inconsistency in the temperature time series limited their use in this study. However, future investigations could use this time series, as the period of record of the new SNOTEL time series is now 20 or more years long [40].

Here we used NDVI [37] to assess changes in canopy. Change in land cover type and the nature of the canopy will influence snowmelt and thus streamflow timing [57,58]. There are other datasets that may be more useful than NDVI, such as OpenET [59].

6. Conclusions

We applied the snowmelt timing and streamflow volume metrics previously proposed [1] for 39 watersheds higher than 2,500 meters across the U.S. state of Colorado. We found that the onset and end of snowmelt-driven streamflow was occurring earlier in almost all of the watersheds. The total annual streamflow increased at a majority of the watersheds, as did the volume before the onset of snowmelt and the volume at the end of snowmelt. These trends were most correlated with winter precipitation, slope (negatively), and latitude. There was correlation with peak SWE for total runoff volume and the volume at the end of snowmelt; these two variables are highly correlated. Due to climatic differences across the domain, in particular drying trends in southern Colorado, winter precipitation was correlated with latitude. Multi-variate regressions illustrated the more highly correlated variables.