Across the western United States, over half the water that annually replenishes rivers and reservoirs originates as snowmelt (Li et al., 2017). The timing and amount of snowmelt that becomes streamflow determines how much water will be available for irrigation, hydroelectric power generation and instream fish habitat (Barnett et al., 2005;Barnhart et al., 2016Barnhart et al., , 2020)). Each year, irrigators and water managers rely on seasonal streamflow predictions issued in the spring that forecast summer water availability (Garen, 1992;Pagano et al., 2004). Many of these regression-based seasonal predictive models are founded on the relationship between April 1st snow water equivalent (SWE) and summer streamflow. Seasonal streamflow predictions are notoriously difficult as they are limited by the information known at the time of prediction and uncertainty in future conditions (Pagano et al., 2004).
To combat this uncertainty, forecasters should perform crossvalidation, be wary over-fitting, use novel metrics and search for the optimal combination of predictors (Garen, 1992;Mendoza et al., 2017). Despite these best practices, studies across the western United States have shown that summer runoff has become more variable, even in places that have not experienced recent shifts in climate (He et al., 2016), and that traditional SWE metrics may become less accurate at predicting seasonal droughts (Livneh & Badger, 2020;Vano, 2020). Recent work in snow-dominated watersheds has addressed the declining accuracy of SWE metrics to predict summer runoff by showing the benefit of adding information about antecedent moisture conditions. These antecedent metrics can improve estimates of spring runoff efficiency or how much snow becomes streamflow as opposed to evapotranspiration (ET) (Castillo et al., 2003;Hammond et al., 2019;Lapides et al., 2022). Specifically, it has been recommended to use mean January runoff as a regression predictor of annual streamflow as a way to represent antecedent groundwater conditions (Brooks et al., 2021) or to use soil moisture data (Harpold et al., 2017;Koster et al., 2010). However, many areas still do not have a long enough soil moisture period of record to integrate into regression-based models; furthermore, ice cover can affect the accuracy of winter runoff measurements (Melcher & Walker, 1992). With the declining efficacy of traditional SWE metrics, it is vital to understand the processes that drive changes in runoff efficiency in seasonally snow-covered areas and develop mechanistic models and metrics that represent these processes at a timescale useful for seasonal predictions.
High-elevation mountainous regions are projected to warm faster than low-elevation areas (Pepin et al., 2015). With warming, much of the western United States may experience a transition from snow to rain (Klos et al., 2014), declining snowpacks (Mote et al., 2018), and earlier and potentially slower snowmelt (Musselman et al., 2017), which could have a profound impact on streamflow, and contribute to the declining ability of April 1st SWE to accurately predict runoff (Barnett et al., 2005;Goulden & Bales, 2014;Livneh & Badger, 2020). With changing melt rates, it is unclear how changes in the timing and rate of snowmelt will affect the amount of melt that goes to ET. Previous work suggests that lower snow accumulation and earlier melt may lead to earlier increases in ET and the start of the growing season (Cooper et al., 2020;Hamlet et al., 2007;Kraft & McNamara, 2022). If this is the case, more snowmelt may go to ET and less to groundwater recharge and annual runoff (Barnhart et al., 2016;Christensen et al., 2021;Goulden & Bales, 2014).
At the same time, other studies have found no change in streamflow with earlier melt (Hammond & Kampf, 2020), suggesting that snowmelt occurring when there is lower evaporative demand offsets the impact of a transition from snow to rain on runoff efficiency (Barnhart et al., 2020;Robles et al., 2021). These opposing results suggest that the hydrologic responses of changes in snow accumulation and melt may be dynamic and vary both regionally and temporally. Three primary mechanisms have been proposed to explain variable responses of streamflow to changes in snow accumulation and melt (Gordon et al., 2022). The first mechanism considers the magnitude of winter vapour fluxes such as sublimation in reducing how much water enters the ground; the second emphasizes the importance of rain and snowmelt input intensity (referred to here as surface water inputs); and the third is based on snow-energy synchrony, which controls how much snowmelt goes to ET. Another way to consider the third mechanism is to explicitly consider snow-vegetation interactions and the synchrony of vegetation greening or ET with snowmelt.
Previous work in the forested southeastern United States found that longer growing seasons and changing vegetation dynamics are driving nonstationary changes in catchment storage and ET (Hwang et al., 2018). In a global analysis, Knighton et al. (2020) found that the alignment of climate and vegetation phenology determines how streamflow responds to forest disturbance. In snow-dominated regions with seasonal offsets in the timing of water inputs and the growing season, they found a dampened streamflow response to changes in forest cover. This period between snowmelt and vegetation green-up, referred to as the vernal window, plays an important role in modulating spring runoff and dictating how snowmelt is partitioned into ET (Barnett et al., 2005;Grogan et al., 2020;Stewart et al., 2005). Many studies have found that vernal windows are lengthening due to earlier snowmelt outpacing changes in the start of the growing season (Grogan et al., 2020); however, this asymmetrical shift largely depends on the climate and vegetative characteristics of a region. These studies highlight the important role that the timing of snowmelt and green-up play on runoff generation processes; however, there is currently no consensus on how to account for these dynamics in regression-based predictive streamflow models. Most new metrics of antecedent moisture focus on absolute water deficits as opposed to timing, and very few explicitly account for changing vegetation dynamics.
Forests impact the magnitude and timing of streamflow generation, yet there is a substantial gap in relating forest processes to downstream water supplies, particularly in irrigated areas (Barnard et al., 2023). In seasonally snow-covered areas, forest structure, such as stem density and canopy cover, affects the amount of precipitation that is intercepted and sublimated (Montesi et al., 2004;Sexstone et al., 2018) and the energy balance of the surrounding snowpack (Kraft et al., 2022;Lawler & Link, 2011;Lundquist et al., 2013;Musselman et al., 2012). Forests also stabilize soils and provide clean water, making them crucial for watershed health and function (Bonan, 2008).
As precipitation regimes change and drought becomes more frequent and severe (Dai, 2013), many regions have started experiencing widespread forest disturbance due to drought-induced mortality, beetle die-off and wildfire (Abatzoglou & Williams, 2016;Adams et al., 2012;Seidl et al., 2017). Each of these disturbances may impact runoff efficiency, though the magnitude and direction of this impact will depend on local geology, topography, climate and vegetation successional classes (Goeking & Tarboton, 2020). Even in regions that have yet to experience large-scale disturbances, changes in local climate may exacerbate or alter how forest structure impacts local hydrology. For example, warmer winters and springs may lead to more canopy interception, early snowmelt and a diminished benefit of canopy shading for snow retention (Dickerson-Lange et al., 2021;Lundquist et al., 2013). Additionally, the amount of forest ET depends on atmospheric evaporative demand and each tree's physiological response to increasing water deficits (Martínez-Vilalta & Garcia-Forner, 2017;Massmann et al., 2019), making it difficult to predict how forests with diverse vegetation will respond to drought. The response of vegetation and ET to warmer and drier conditions has been shown to determine how watersheds and streamflow respond to and recover from drought (Avanzi et al., 2020;Maurer et al., 2022), or, in some cases, explain why watersheds experience a hydrologic regime shift following severe drought (Peterson et al., 2021).
The complex interplay of climate and forest structure, and the effects of vegetation-climate interactions on the hydrologic cycle, make it challenging to capture these dynamics in real time, let alone ahead of time, as needed for seasonal predictions. Thus, few seasonal models integrate information about forest dynamics. Those rare exceptions often use a relatively static metric of percent forest cover, which fails to capture intra-annual changes in vegetation dynamics (Hernandez et al., 2018;Sun et al., 2014).
Here, we address this knowledge gap by introducing a new metric, the Vegetation-Water Alignment Index (VWA), to investigate whether a metric that leverages dynamic vegetation behaviour can improve runoff predictions following a multi-year drought. The VWA quantifies the alignment between vegetation greenness and water availability, thereby both capturing changes in forest structure and providing an inference into water-energy synchrony. Summer runoff and low flows can retain a memory of the previous year's snowpack and storage deficits (Brooks et al., 2021;Godsey et al., 2014;Lapides et al., 2022). Therefore, we hypothesize that when annual snowmelt and rain inputs are more aligned with seasonal vegetation greenness, less snowmelt will go to runoff and recharge, which will reduce the following spring's runoff efficiency. Using a case study from the Greater Yellowstone Ecosystem, we select and validate new predictive models and show how including metrics of evaporative demand and vegetation dynamics improve seasonal predictions of summer runoff. These results emphasize the importance of explicitly accounting for vegetation dynamics as the mechanisms and relationships used to predict streamflow change with droughts and declining snowpacks.
The Henrys Fork of the Snake River (hereafter referred to as the Henrys Fork) is a major headwater tributary of the Snake River in eastern Idaho, USA (Figure 1b,c). The Henrys Fork is located in the Greater Yellowstone Ecosystem and sourced by three subwatersheds (from here on referred to as watersheds)-the Upper Henrys Fork, Fall River and Teton River-that collectively supply $25% of the water supply in the upper Snake River basin and the nearly $10 billion agricultural economy that exists in the basin (Idaho Water Resource Board, 2009;Van Kirk et al., 2019). Like many river basins across the western United States, the Henrys Fork is a seasonally snow-covered, forested headwater system that supplies water to downstream users and is expected to see shifts in climate and runoff reliability (Hostetler et al., 2021). As such, this region serves as a case study to investigate the effects of changing vegetation, snow and energy alignment on seasonal runoff predictions and runoff efficiency.
The hydrology of both the upper Snake River basin and the Henrys Fork reflects the complex volcanic geology of the region that formed in the wake of the migration of the Yellowstone Hotspot (Pierce et al., 2007). Due to the unique combination of geology, topography and hydrology, each watershed exhibits distinct streamflow generation processes. As seen in Figure 1a, the Upper Henrys Fork (HF) has the greatest groundwater contribution with higher winter runoff, lower and earlier spring runoff and higher summer baseflows. A majority of the groundwater that feeds the Upper Henrys Fork discharges from a series of springs that emerge at the foot of the Yellowstone Plateau, a nearly 400-m thick rhyolite tuff from the recent Lava Creek eruption ($600 Ka; Pierce & Morgan, 1992).
Snowmelt on the plateau recharges the springs' aquifer, which has a response time of approximately 3 years (Benjamin, 2000). The spring discharge joins the rest of the Upper Henrys Fork, which flows through a legacy caldera before descending onto the Eastern Snake River Plain. The Upper Henrys Fork has the lowest average elevation, followed by Fall River and Teton River, as seen in Figure S1.
In contrast to the strong groundwater influence in the Upper Henrys, Teton River is strongly snowmelt-dominated, with most of the water supply coming from snow that falls in the Teton Range.
Figure 1a shows the high spring flows and lower summer baseflows indicative of snowmelt-driven runoff. The remaining watershed, Fall River, is sourced by a mix of groundwater springs and snowmelt, resulting in a more intermediate hydrograph shape (Figure 1a).
The climate of the study area transitions from warm continental at low elevations to dry-summer subarctic at high elevations (Kottek et al., 2006). Across all three watersheds, there is a sharp precipitation gradient with around 300 mm of annual precipitation in the valleys and up to 2000 mm in the high-elevation mountains, with anywhere from 50% to 80% falling as snow, depending on elevation. The Greater Yellowstone Climate Assessment projects that this region will receive more annual precipitation, but a lower snow fraction, in the coming decades (Hostetler et al., 2021).
A large portion of each watershed is forested, as seen with the green shading in Figure 1b. The composition of these forests varies by elevation and aspect but are largely dominated by lodgepole pine (Pinus contorta), Engelmann spruce (Picea engelmannii), Douglas-fir (Pseudotsuga menziesii), whitebark pine (Pinus albicaulis) and subalpine fir (Abies lasiocarpa). These forests experienced substantial disturbance before the 1989 start of this study due to a combination of extensive timber harvesting and the 1988 Yellowstone Fires (Turner et al., 2003). Though most of the fire perimeter was outside of the watershed, with less than 5% of both the Upper Henrys and Fall River watersheds affected (Figure S2), beetle infestation following the fires led to additional delayed mortality in the early 1990s (Ryan & Amman, 1996). While there has been a recent resurgence of beetle mortality in the eastern Greater Yellowstone Ecosystem, the Henrys S1 gives description of each SNOTEL site. Green shading indicates land classified as evergreen, mixed or deciduous forests Dewitz, (2023). The colour of each SNOTEL site indicates the watershed that site is associated with and the shading indicates the relative elevation of the site where a darker colour indicates a lower elevation. (c) Regional site map of the Henrys Fork, a major headwater tributary in the Columbia-Snake basin.
Fork has experienced only small-scale disturbances in the past 20 years.
Water resources in the Henrys Fork are highly regulated by irrigation storage, delivery and return flow (Morrisett et al., 2023). In this study, we use natural streamflow (Q nat ) to refer to unregulated runoff in each watershed (Van Kirk, 2020), calculated as:
where Q reg is regulated streamflow at the USGS gage, diversions are surface water diversions recorded daily by Idaho Water District 01, ΔS is change in upstream reservoir storage, E res is upstream reservoir evaporation and injection is input to streamflow upstream of the gage via delivery of water from other watersheds through canals or via exchange wells, which pump groundwater directly into rivers (U.S.
Bureau of Reclamation and Idaho Water Resource Board, 2015).
Direct precipitation onto the reservoir is included in the reservoir evaporation term, making E res negative when precipitation exceeds evaporation. As each watershed has different infrastructure, the application of Equation ( 1) varies slightly between each watershed. A full description of these differences can be found in Supporting Information Text S1. The normalized discharge shown in Figure 1a is the mean daily natural streamflow on any given day of the year divided by the 1989-2023 annual average daily natural streamflow.
The precipitation, SWE and snow depth data used in this study come from all of the Natural Resources Conservation (NRCS) Snow Telemetry (SNOTEL) sites in and near the study watersheds (Figure 1b). Key topographic and vegetative characteristics of each site are given in Table S1. As advised by NRCS, the temperature data were corrected to resolve a known bias associated with a change in sensor (Atwood et al., 2023). Since hourly instantaneous data was not available for the entire period of record for all sites, we applied the correction equation to daily summary statistics. Over the relevant range of temperatures for a given day, the correction is nearly linear, indicating that applying the correction directly to the daily mean summary statistic is close to the true mean. We performed the analyses outlined below with both the daily mean and the mid-point of the range in daily minimum and maximum temperatures. We found no difference, and therefore, present the results using mean daily temperatures.
Using the SNOTEL SWE data, we calculated daily melt rates as the decrease in SWE. We also partitioned daily precipitation into rain and snow by classifying days with precipitation and an increase in snow depth or SWE as 'snowy' and days with new precipitation but no change or a decrease in snow depth or SWE as 'rainy' (Jennings & Molotch, 2019). Using this melt and rain data, we calculated daily surface water inputs (SWI) as melt plus rain to represent when liquid water enters the system (Kormos et al., 2014).
To compare sites, we will use the term 'lower' to refer to the SNOTEL stations that have an elevation between 1900 and 2100 m, 'mid' to refer to those between 2100 and 2300 m and 'high' to refer to the stations with an elevation greater than 2300 m.
The normalized difference vegetation index (NDVI) is a metric of vegetation greenness commonly used to describe growing season dynamics (White et al., 2009), leaf area index (Tewari et al., 2012), forest density and disturbance (Jin & Sader, 2005), and in some areas, NDVI has been shown to have a strong relationship with ET (Goulden et al., 2012). As a broad measure of vegetation dynamics, we use the 250-m MODIS (Moderate Resolution Imaging Spectroradiometer)
16-day NDVI data products from both the Terra and Aqua satellites (MOD13Q1 and MYD13Q1, collectively referenced as M*D13Q1).
Together, these products produce a 20-year time series that are widely used and publicly accessible through open-source software.
Before release, the 16-day NDVI products are internally processed for atmospheric correction and cloud cover, with the highest quality image and highest NDVI value selected for each pixel with the corresponding day of acquisition recorded (Didan et al., 2015). To further process this data, we used the M*D13Q1 reliability index to filter out any remaining low-quality pixels associated with cloud cover or interference other than snow cover. After filtering out low-quality pixels, we used the NDVI and the day of image acquisition for each pixel to generate a time series of average watershed NDVI, average watershed annual maximum NDVI and time series of NDVI extracted at the pixel containing each SNOTEL site location.
As a spectral index, NDVI is affected by snow cover extent (Wang et al., 2013), suggesting that in seasonally snow-covered areas, winter NDVI values will represent both vegetation greenness and snow cover. To retain an annual time series and incorporate insight associated with snow-covered extent and frequency (Hall & Riggs, 2007), we deliberately retained values flagged as snow-covered if the image was otherwise high quality. Once processed as described above, NDVI values show some variability related to early winter snow accumulation and depth, but are fairly consistent during midwinter months (Figure S3).
To investigate whether predictors of summer runoff in this region have changed, we used modified time trend analysis.
Time trend analysis can be used to isolate the impacts of changes in climate or vegetation on streamflow, such as the effects of forest disturbance on watershed yield (Goeking & Tarboton, 2022; Manning et al., 2022; Zhao et al., 2010). This is done by fitting a regression model to a predisturbance calibration period and testing whether the difference between the observed and predicted values for the post-disturbance period is significantly non-zero. If they are non-zero, this implies that the predictive relationships have changed between the calibration and test periods. For example, this method might use annual precipitation and temperature to predict annual runoff before and after a disturbance. Here, we use a modified approach to quantitatively investigate whether the relationships used to predict summer streamflow during and after drought differ from pre-drought relationships. In this modified approach, we use seasonal predictive models of summer runoff and test whether a multi-year drought changes runoff predictability.
To delineate pre-drought, drought and post-drought periods, we performed linear breakpoint analysis on a time series of the Palmer Drought Severity Index (PDSI) averaged across all three watersheds (accessed via ClimateEngine.org; Huntington et al., 2017). We used the resulting breakpoints to define pre-drought, drought and postdrought periods and investigate non-stationarity in predictive relationships as detailed below.
For the pre-drought predictive model, we use April 1st SWE and mean January runoff to predict summer runoff for each watershed, as given in Equation ( 2) where log represents the base-e logarithm.
Here, runoff refers to area-normalized Q nat (mm). Summer runoff is total runoff from April 1st-September 30th and Jan. runoff is mean daily January runoff (mm). April 1st SWE is the average from all SNO-TEL sites in or near each watershed (Figure 1b). We define summer runoff as runoff between April and September due to peak irrigation demand in the region occurring in July and August: this acknowledges the importance of predicting runoff through the irrigation season (past the commonly used April-July timeline). We include January runoff as it has been demonstrated to represent antecedent groundwater conditions and improve streamflow predictions in snow-dominated watersheds (Brooks et al., 2021). Although other factors inevitably influence summer runoff, we emphasize that this is a seasonal predictive model, and therefore, we can only include variables known by April 1st, when seasonal predictions are made. We used the predrought water years of 1989-1998 to calibrate the models using Equation ( 2). The Teton River model has a slight modification where we used April 1 SWE without a log transformation to improve the normality and leverage of the residuals. The normality of the residuals was confirmed with a Lilliefors test (Lilliefors, 1967) and the homoscedasticity of the residuals with scatterplots of residuals versus summer runoff and water year.
We then applied the pre-drought model to the drought and postdrought periods and calculated ΔQ nat as observed natural runoffpredicted natural runoff, so ΔQ nat is negative when the model overpredicts runoff. We used a one-sample t-test to test whether ΔQ nat for the water years in each period is non-zero. If ΔQ nat is non-zero, this suggests the accuracy of summer runoff predictions has changed due to something other than April 1st SWE and antecedent groundwater. To facilitate comparison of ΔQ nat across watersheds, we calculated the percent deviation of observed from predicted runoff using Equation (3).
Next, we investigated changes in climate and vegetation that may impact runoff efficiency. To do this, we ran Mann-Kendall trend analysis on a variety of temperature, precipitation and vegetation metrics.
The nonparametric Mann-Kendall trend test detects monotonic increasing or decreasing trends over time (Helsel et al., 2020). In this study, we performed this analysis over the entire period of record
To test whether adding information about vegetation dynamics improves runoff predictions, we introduce a metric that captures changes in both vegetation and water availability, or the water- growing season alignment. To characterize alignment, we use a similarity index that quantifies how in-phase calendar year NDVI is with SWI, which we will refer to as the VWA. As described in Section 2.2.2, SWI is the combination of snowmelt and rain, which represents when liquid water enters the ground. As with previous studies that calculate similarity indices using different climatic variables (Apurv & Cai, 2020;Hale et al., 2023;Woods, 2009), we first fit a sine curve to calendar year SWI at each SNOTEL site (Equation 4). For the purpose of capturing key characteristics of the seasonal periods, we assume that a single sine curve can reasonably fit SWI data, an assumption often used for precipitation, temperature and ET data, particularly in snowcovered regions (Apurv & Cai, 2020;Berghuijs & Woods, 2016;Hale et al., 2023;Milly, 1994;Potter et al., 2005;Woods, 2009). Figure 2a graphically summarizes Equations (4-6) and their interpretation.
In Equation ( 4), SWI is mean annual SWI, δ SWI is the dimensionless amplitude, t is time in days and s SWI is the phase shift in days. Text S2 in the Supporting Information provides an overview of using a sine function in this context and a discussion of Equations (4 and 5) within the context of the general sine form.
To integrate the annual variability of vegetation dynamics into this metric, we also fit a sine curve to calendar year MODIS NDVI extracted to each SNOTEL site (Equation 5), where δ NDVI is the dimensionless amplitude of the sine curve and s NDVI is the phase shift in days.
Here, we account for interannual variability in NDVI by multiplying by half of NDVI MAX instead of mean annual NDVI. Previous work has shown that maximum NDVI is a more reliable way of compositing NDVI given the heavy influence of cloud and atmospheric contamination on mean NDVI (Holben, 1986;Zhang et al., 2017). For both SWI and NDVI, sine functions were fit to annual data using the iterative, non-linear least-squares function 'nlsLM' in the R package minpack.lm (Elzhov et al., 2023). Figure 2a shows an example of how sine curves fit annual cycles of SWI and NDVI.
The phase and amplitude of each annual curve can be related using a similarity index, as given with Equation ( 6), which is weighted by the amplitude of NDVI and quantifies whether NDVI is in phase with SWI (Berghuijs et al., 2014;Woods, 2009). The VWA, or the similarity index of SWI and NDVI, ranges from close to À1 to 1, depending on δ NDVI . A value closer to À1 indicates that NDVI is out-of-phase with SWI and +1 that NDVI and SWI are in-phase.
Here, sgn refers to the sign (positive or negative) of δ SWI , which is used to assess whether NDVI and SWI periods exhibit the same seasonality (i.e., if the peaks both fall in the first or second 6 months of the year or not). To capture dynamic changes in vegetation, Equation ( 6) scales the distance of the phase shifts by δ NDVI . This is done so that when there is little annual variability in NDVI (a δ NDVI closer to 0), VWA will also be closer to 0. When NDVI shows strong annual variability, and therefore, has a higher δ NDVI , this scaling will move the value closer to ±1, with the sign dependent on the difference in phase shifts. Given the seasonal snow cover and energy limitations on the growing season, we do not observe the full range of possible VWA values in this system. Figure 2b gives an example using the NDVI and SWI curves associated with the lowest and highest observed VWA at the Island Park (IP) SNOTEL site. This example demonstrates how a shift in VWA can reflect changes in timing of peak SWI and greenness as well as the amplitude of NDVI. For a conceptual overview of the components of SWI and NDVI sines curves and a visualization of how shifts in each variable's amplitude and phase shift affects VWA, see Figures S4 and S5. Using Equations (4-6), we calculated VWA for each SNOTEL site in or near the watershed boundaries, except for Crab Creek (CC). We use the CC SNOTEL site for Upper Henrys meteorologic averages, but did not calculate VWA for that station as it did not have snow depth data available for the entire record.
We propose using VWA as a metric that integrates information about water-growing season alignment, changes in peak NDVI and changes in the annual distribution of SWI, all of which may influence annual runoff efficiency in seasonally snow-covered watersheds (Gordon et al., 2022). Given the memory that many watersheds have of antecedent moisture conditions (Brooks et al., 2021;Castillo et al., 2003), we propose that the alignment of NDVI with SWI for a given calendar year can be used to represent the relative storage state of the catchment at the end of the year and therefore be used to scale the following spring's seasonal runoff predictions. Figure 2b shows our hypothesis that when NDVI and SWI are more aligned, less melt and rain will go towards runoff and recharge. Instead, when peak SWI is more closely aligned with the start of the growing season, more of the snow meltwater entering the soil will be taken up by plants and lost to ET, creating a storage deficit. This is visualized in Figure 2b using arrow size to represent the relative flux of SWI, plant water uptake and runoff and recharge. Following a closely aligned year, we hypothesize that the memory of this reduced recharge will be carried over through a storage deficit that leads to lower runoff efficiency for a given snowpack the following spring. We propose that this lower runoff efficiency is due to the storage deficit needing to be refilled before runoff is produced (McDonnell et al., 2021).
To investigate whether VWA, or some other new metric, adds predictive power compared with traditional snowmelt-based models of streamflow (e.g., Equation 2) in the post-drought period, we performed statistical model selection using Akaike's information criterion (AIC) to rank models and select a new predictive model for the post-drought period in each watershed. We created models for the post-drought period to ensure consistent periods of record for all satellite and sensor data and to ensure the model is trained with years with similar forest cover and vegetative state. Here, we refer to the selected regression models as the 'new' models and the predictive models fit with April 1 SWE and January runoff (Equation 2) as the 'traditional' models.
To select the best new models for each site, we provided a set of candidate predictor variables to generate potential models and then ranked them according to AIC, which measures the variation explained by each model and the complexity of the model relative to the other candidate models (Claeskens & Hjort, 2008). All candidate predictors are listed and defined in Table S2. Given the short postdrought period of record, we used the small-sample-corrected AIC, referred to as AICc. Model selection was performed using the MuMIn:
Multi-Model Inference package in R version 4.2.3 (Bartoń, 2023).
After model selection, we performed leave-one-out cross-validation to test the predictive performance and bias of the new models. To do this, we iterated through each predicted year of the post-drought period (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023), fit a model using all but the selected year, and then used that model to predict runoff for the year not included in the model fit. To assess model performance, we calculated Nash-Sutcliffe efficiency (NSE), model percent bias (PBIAS) and root mean square error (RMSE) using the observed and predicted values. To facilitate comparison across watersheds, we scaled RMSE by the standard deviation of observed summer runoff, referred to here as RMSE sd (Moriasi et al., 2007). To compare the improvement offered by these new models, we also performed cross-validation on the traditional model using April 1 SWE and January runoff. By comparing performance metrics between the models, we can quantify the predictive improvement offered by the new model over the traditional one. Metrics of the pre-drought traditional model fits are given in Table S3.
Figure 3c,d shows that the time series of percent deviation of observed from expected runoff is below zero in these two watersheds for most of the drought period, meaning that observed runoff was consistently lower than expected in the drought period. Figure 3e shows the average ΔQ nat for each watershed, with significance noted for periods when ΔQ nat was significantly non-zero.
The average ΔQ nat improved from the drought to post-drought period in Fall and Teton Rivers. While the model over-predicted postdrought runoff, the over-predictions were insignificant. In the Upper Henrys, we see significant over-predictions and declining postdrought model performance ( p < 0.1). The observed over-prediction in the drought and post-drought periods suggests that the system has experienced mechanistic changes affecting runoff efficiency not captured by the pre-drought relationships between April 1 SWE, January runoff and summer runoff.
To investigate other potential drivers affecting the reliability of the historic stream-snowpack relationship, Figure 4 Comparison of precipitation metrics, including water year total annual precipitation, day of peak SWE, annual max SWE and snow fraction (peak SWE/annual precipitation), reveal no significant trends in any watershed (Figure 4).
Analysis of vegetation metrics shows that NDVI MAX is significantly increasing over both satellite periods of record in the Upper Henrys Fork. We also see a significant increase over the AVHRR period of record in Fall River. Trend analysis of the 1989-2020 fractional tree cover record shows significant increases in percent tree cover in all watersheds in all three periods. Both fractional tree cover and NDVI MAX suggest that forest greenness and density have changed over the course of the study period in the Upper Henrys and Fall River watersheds. In the Upper Henrys Fork and sometimes in Teton River, the lower elevation SNOTEL point measurements have a lower (more negative) VWA than the higher elevation stations (Figure 5a,c). A more negative VWA signifies that vegetation greenness is more out of alignment with snowmelt and rain entering the system at lower elevations than higher elevations. This may reflect more mid-winter melt events and earlier spring melt due to a lower snow fraction and warmer winter temperatures at lower elevations. This elevational relationship does not occur in Fall River (Figure 5b); however, it is worth noting that the lowest elevation SNOTEL site in Fall River is at the same elevation as the mid-elevation sites in the Upper Henrys Fork (Figure S1).
The dashed watershed-average VWA is often higher than what is captured at the SNOTEL sites, reflecting a closer growing season-water input alignment. However, in each watershed, the spatially averaged VWA is lower than the SNOTEL VWA in 2007 and 2015. Both of these are low snow years where contrasting patterns in the amplitude and phase shift of annual NDVI at the watershed-scale likely reflect early melt and water-limited phenological responses of other vegetation types (shrubs, grasses, crops, etc.) that are present throughout the watershed but not dominant at SNOTEL sites. Underlying all the VWA patterns, we found that annual SWI is far more variable than NDVI, but in many areas, peak NDVI is shifting earlier in the year (Figures S6-S8).
Post-drought AICc model selection shows that VWA is an important predictor of summer runoff. In Teton River and the Upper Henrys Fork, we selected the model with the lowest AICc, both with four fitted parameters. In Fall River, the top three models had similar AICc and we chose the most parsimonious, with only four parameters.
Figure 6 shows the new predictive models selected for each water- To test the predictive performance of these models, we report results from leave-one-out cross-validation using the 2007-2023 period of record (Figure 7), highlighting predicted versus observed runoff for each year using the new and traditional predictive models.
To facilitate comparison across watersheds, we scaled RMSE by the standard deviation of observed summer runoff, referred to here as RMSE sd . An RMSE sd of less than 0.5 indicates good model performance (Moriasi et al., 2007). The summary metrics are given in Table 1.
Overall, we see the best improvement and model performance for each year in the post-drought period in the Upper Henrys with an NSE of 0.79. We also see improvement in Fall River, with an NSE of 0.70 for the new model compared with 0.56 for the traditional model.
Teton River shows moderate improvements and performance with an NSE increase from 0.43 to 0.6. In all watersheds, we see low model
Table showing results of Mann-Kendall trend analysis of temperature, precipitation and vegetation metrics for all three watersheds. Each box is coloured by the Kendall tau of the trend where a darker green indicates a stronger increasing trend, and brown, decreasing. Black boxes group the metrics categorically. Significant trends are labelled with *** at p < 0.001, ** at p < 0.01, * at p < 0.05 and variables with periods of record shorter than the 1989-2023 periods of record are indicated with a diamond (◊). Advanced very-high-resolution radiometer (AVHRR) and MODIS refer to the corresponding satellite for each normalized difference vegetation index (NDVI) dataset. biases with slight model under-prediction with the new models. Finally, we see improvements in RMSE sd in Fall River and the Upper Henrys Fork, dropping RMSE sd for both sites close to or under 0.5. 4 | DISCUSSION 4.1 | Non-stationarity of predictive relationships, climate and vegetation dynamics requires new models
In this study, we investigated whether historic relationships between April 1 SWE, winter and summer streamflow still accurately predict summer runoff following the early 2000s severe drought. We analysed model performance before, during and after drought since it is well documented that drought reduces runoff efficiency and decreases the accuracy of seasonal predictive models (Saft et al., 2015;Tian et al., 2018). Furthermore, recent work looking at recovery following Australia's Millennium Drought suggests that some watersheds shift to a new state following severe drought, and that drought recovery is not just linked to catchment wetness, but also changes in vegetation that cause more precipitation to go to transpiration (Peterson et al., 2021).
Time trend analysis results show that the pre-drought relationships over-predict summer runoff during the drought (FR and TR) and post-drought (HF) periods. Before and after the drought, this region has not experienced any widespread changes in annual precipitation, snow accumulation or timing of melt (Figure 4), suggesting that the worsening performance in the Upper Henrys Fork is not due to changes in water inputs. What has changed is warmer late spring and summer temperatures (May through September) in the Upper Henrys Fork and Teton River watersheds as well as vegetation greenness and fractional tree cover. The latter two vegetation metrics increased the most in the Upper Henrys Fork and to lesser degrees in Fall River and Teton River (Figure 4). Our analysis shows that these changes in temperature and vegetation dynamics impact runoff efficiency as all three new models selected for the post-drought period include a metric of the previous summer's evaporative demand (ETr) and vegetation dynamics (VWA).
Collectively, these results affirm that accounting for changing snow dynamics alone may fail to capture important mechanisms underlying changes in runoff efficiency and the relationships used to predict streamflow. These findings also suggest that our understanding of future runoff dynamics hinges on our understanding of climatevegetation interactions.
The VWA proposed here is a metric that integrates information about changes in annual SWI, snow fraction, snowmelt rate, seasonality of precipitation, snow cover, timing of vegetation green-up, peak greenness, length of growing season and the synchrony of SWI and NDVI.
Many of these elements are predicted to change with warmer temperatures and a snow-to-rain transition (Cayan et al., 2001;Klos et al., 2014) and may impact runoff efficiency across multiple timescales (Barnett et al., 2005). Each new model includes a SNOTEL VWA as a negative predictor of summer runoff, with VWA showing up in nearly all the top-ranked models (Figure 6 and Tables S4-S6).
The importance of VWA as a negative predictor supports our hypothesis that the more closely aligned NDVI and SWI are for a given calendar year, the lower the following spring's runoff efficiency will be (Figure 2b). This decreased efficiency may be the result of more of the following year's snowmelt going towards ET and less going to recharge groundwater and deep soil moisture storage, leading to drier conditions going into the following spring and a larger storage deficit that must be overcome before runoff is produced (Castillo et al., 2003;Harpold et al., 2017;McDonnell et al., 2021;Milly, 1994). S1. runoff occurs after storage deficits is refilled (McDonnell et al., 2021).
The significance of this work goes beyond just a mechanistic understanding of streamflow generation; it also guides how we can integrate these processes into regression-based predictive models at the time scale needed for seasonal predictions.
In the context of the proposed mechanisms driving variable streamflow response to changes in snowpack, our results support the importance of capturing changes in energy-water synchrony (Gordon et al., 2022). While not explicitly modelling energy synchrony, the VWA indirectly accounts for spring and summer temperatures by using NDVI, which reflects the timing of snow disappearance, length of the growing season and provides inferences of ET without any of the assumptions often used in ET models, all of which reflect energy changes. Additionally, this metric also captures changes in the timing and magnitude of SWI inputs related to changes in melt rates, an important insight given that faster melt rates have been shown to produce larger streamflow responses (Barnhart et al., 2016(Barnhart et al., , 2020)).
The new models improved streamflow predictions the most in the Upper Henrys Fork (Figure 7). This improvement may be due to a combination of factors, including the groundwater-dominated hydrologic regime, greater forest regrowth following disturbance, and the watershed's topography. The Upper Henrys Fork is the only one of the three watersheds that has experienced both significant warming and higher NDVI MAX over both satellite periods of record (Figure 4).
This suggests that, as expected, the more change a watershed has experienced, the greater the need to revise traditional predictive relationships. Compounding this change in the Upper Henrys Fork is the memory of past deficits retained in the groundwater system (Benjamin, 2000).
Groundwater-dominated systems experience a larger absolute reduction in summer streamflow following changes in snowmelt timing and warming than their snow-dominated counterparts (Mayer & Naman, 2011;Tague & Grant, 2009). This is especially true in porous, young volcanic landscapes (Tague et al., 2008), such as the calderas and tuffs that underly the Upper Henrys. Not only do deficits accrue due to the slower drainage rates but also these porous substrates may have a large subsurface storage capacity. This stored water can modulate the hydrologic impacts of climate-ET interactions if this water is accessible to plants during dry periods (Garcia & Tague, 2015;Klos et al., 2018). This increased storage allows plants to deplete deeper water stores during multi-year dry periods, resulting in lower streamflow generation and greater deficits (Hahm et al., 2020(Hahm et al., , 2022)).
Improved predictive performance in Fall River, a watershed underlain by a mix of young volcanic and glacial deposits with a more mixed groundwater-snowmelt system, also supports the important role underlying geology and subsurface storage capacity may play in modulating runoff efficiency. Our analysis shows more moderate improvements in the snowmelt-dominated Teton River, perhaps indicative of the larger role of real-time melt dynamics in snowmelt-driven systems. The importance of capturing snowpack evolution and melt in near-real time has motivated efforts to operationalize physically based models to provide more accurate short-term forecasts (Meyer et al., 2023). However, these forecasts have much shorter lead times and depend on accurate meteorological forecasts, suggesting the value of also continuing to improve our understanding of snowmelt runoff efficiency on a seasonal scale.
In each new seasonal runoff model, SNOTEL sites in the low to intermediate elevation zones (1900-2100 m) were selected as the best predictors of streamflow (Figure 7). Previous work has shown that ET and vegetation greenness in intermediate elevation zones are more responsive to changes in annual precipitation and snowpack (Christensen et al., 2008;Kraft et al., 2022;Tague & Peng, 2013;Trujillo et al., 2012). This finding suggests that mid-elevation forests, which are able to optimally use available water, may exert a larger influence on storage deficits and runoff efficiency than forests that are more strongly water-or energy-limited. As warming shifts the timing between snowmelt and the growing season, annual and interannual hydrologic response will depend on how mid-elevation forests affect subsurface storage as they transition from energy-to waterlimited growing seasons.
This study provides a case study and proof of concept for the utility T A B L E 1 Leave-one-out performance metrics of traditional and new models for the post-drought period by watershed.
Watershed Model NSE PBIAS RMSE sd Upper Henrys Traditional 0.564 À0.308 0.640 New 0.790 À0.129 0.447 Fall Traditional 0.563 À0.642 0.641 New 0.700 À0.427 0.534 Teton Traditional 0.432 0.268 0.731 New 0.601 À0.298 0.613 Abbreviations: NSE, Nash-Sutcliffe efficiency; PBIAS, percent bias; RMSE, root mean square error.
In this region, denser forests combined with warmer springs and summers set the scene for changing snow-vegetation interactions and greater evaporative demand to have a noticeable effect on local hydrology. Follow-up studies should investigate the utility of VWA in watersheds that have seen more widespread changes in snow accumulation and melt and/or recent vegetation disturbances. The similarity index of SWI and NDVI integrates many different aspects of snow-vegetation dynamics, and it would be valuable to investigate whether the effects a disturbance such as wildfire has on both SWI and NDVI translates to a meaningful shift in VWA that corresponds to a change in runoff efficiency.
Here, we present a case study looking at how the relationships used in seasonal runoff predictive models have changed over the past three decades in the Henrys Fork of the Snake River, a forested headwater system located in the Greater Yellowstone Ecosystem. We see that predictive relationships used prior to a severe drought in the early 2000s over-predict runoff during and after the drought in two of the three watersheds, and that the remaining watershed shows declining performance in the post-drought period. Our analysis of climate and vegetation change over the same period shows that, although this region has not experienced widespread shifts in precipitation or snow accumulation since the late 1980s, late springs and summers are generally warmer and forests are greener with more tree-covered area as forests recover from multiple disturbances. These results suggest the changing snowpack-streamflow relationship is affected not just by changes in snow accumulation and melt, but also by vegetation dynamics, an aspect of watershed hydrology not historically captured in regression-based predictive models. To this end, we propose a new metric, the VWA, that quantifies water availability-growing season alignment. New predictive models show that this index, along with a metric of atmospheric demand, are significant predictors of runoff and can be used to improve the runoff predictions in the post-drought period. Furthermore, VWA for a given calendar year is a significant negative predictor of the following year's summer runoff. This supports the hypothesis that the more synchronized the timing of water inputs are with the growing season, the lower future runoff efficiency will be. This finding emphasizes that our understanding of runoff dynamics and our ability to predict runoff depends not only on our ability to model snowpack but also on our understanding of and ability to model how vegetation responds to warming and water stress and how these stress responses affect groundwater storage and recharge.
The VWA proposed here is a step forward in integrating vegetation dynamics into seasonal predictive models and further work should be done to assess the applicability of this metric outside of the case study presented here.
10991085, 2024, 6, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/hyp.15211 by Sarah Newcomb -Idaho State University , Wiley Online Library on [29/06/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
of 19