Why does snowmelt-driven streamflow response to warming vary? A data-driven review and predictive framework

Climate change is altering the seasonal accumulation and ablation of snow across mid-latitude mountainous regions in the Northern Hemisphere with profound implications for the water resources available to downstream communities and environments. Despite decades of empirical and model-based research on snowmelt-driven streamflow, our ability to predict whether streamflow will increase or decrease in a changing climate remains limited by two factors. First, predictions are fundamentally hampered by high spatial and temporal variability in the processes that control net snow accumulation and ablation across mountainous environments. Second, we lack a consistent and testable framework to coordinate research to determine which dominant mechanisms influencing seasonal snow dynamics are most and least important for streamflow generation in different basins. Our data-driven review marks a step towards the development of such a framework. We first conduct a systematic literature review that synthesizes knowledge about seasonal snowmelt-driven streamflow and how it is altered by climate change, highlighting unsettled questions about how annual streamflow volume is shaped by changing snow dynamics. Drawing from literature, we then propose a framework comprised of three testable, inter-related mechanisms—snow season mass and energy exchanges, the intensity of snow season liquid water inputs, and the synchrony of energy and water availability. Using data for 537 catchments in the United States, we demonstrate the utility of each mechanism and suggest that streamflow prediction will be more challenging in regions with multiple interacting mechanisms. This framework is intended to inform the research community and improve management predictions as it is tested and refined.


Introduction
Snowmelt-driven streamflow is a critical-and increasingly vulnerable-freshwater resource for downstream environments supporting agriculture, municipalities, and native ecosystems (Viviroli et al 2007, Mankin et al 2015, Li et al 2017, Immerzeel et al 2020, Viviroli et al 2020. In the United States (US) alone, snowmelt runoff from high elevation, mountainous (upland) catchments serves water to over 60 million people and supports billions of dollars in economic productivity (Barnett et al 2005, Sturm et al 2017 with additional impacts on ecosystem health (Allan and Castillo 2007), the extent and severity of wildfires (Westerling et al 2006, Holden et al 2012, flood risk (Hamlet and Lettenmaier 2007), and reservoir management (Ehsani et al 2017). Snowmeltderived water resources are perhaps most important in semi-arid regions where the predictable cycle of seasonal snow accumulation and melt, together with vast networks of physical and legal infrastructure, enable development in otherwise water-limited environments (Church 1933, Meybeck et al 2001, Qin et al 2020. Large interannual variability in snowmelt water supply is driven primarily by variability in winter precipitation (P), but also by 30% or higher variability in runoff efficiency across basins (Mote 2003, Harpold et al 2012, Brooks et al 2021. In spite of this variability, decades of observations have resulted in reasonably skilled water supply models dependent on readily observable metrics such as 1 April Snow Water Equivalent (SWE) and temperaturedriven melt models (Pagano et al 2004). The last few decades, however, have seen a large decrease in the skill of these models (Pagano et al 2004) leading to diverse suggestions about limitations in understanding and predictability of water supply (Milly et al 2008, Montanari and Koutsoyiannis 2014, Milly et al 2015. This review focuses on the mechanisms that give rise to spatial and temporal variability in runoff efficiency as snow cover changes in response to a warming climate.
Seasonal snowmelt-driven streamflow is one of the fastest changing aspects of the global hydrologic cycle in response to climate change (Musselman et al 2017). Warmer winter and spring temperatures are decreasing the fraction of precipitation falling as snow (fs, Knowles et al 2006, Klos et al 2014, delaying the initiation of consistent snow cover, increasing soil frost (Burakowski et al 2008, Wobus et al 2017, increasing water vapor exchanges between snowpack and the atmosphere (Harpold and Brooks 2018, Sexstone et al 2018), advancing the timing and slowing the rate of snowmelt (Musselman et al 2017), and decreasing the persistence of snow cover (Stewart 2009). In contrast to a general consensus on seasonal snow cover decline under a warming climate, predictions about streamflow response are much more diverse. For example, recent research suggests that changing snow dynamics may or may not lead to earlier streamflow timing (Stewart et al 2005, Moore et al 2007, Fritze et al 2011, and may or may not alter annual runoff efficiency (Berghuijs et al 2014, McCabe et al 2018. These uncertainties in how changes in seasonal snow dynamics will influence streamflow on the scale at which management decisions are made (e.g. 100-10 000 km 2 ) remain and have been the focus of research seeking to improve empirical and model-based forecasting tools (Huang et al 2017, Siirila-Woodburn et al 2021.
The conflicting reports on snowmelt-driven streamflow response to a changing climate result in large part from regional differences in the mechanisms controlling seasonal snow accumulation and ablation. To better anticipate and adapt to changing upland snow resources requires a consistent and testable framework to characterize variability in the dominant processes that drive differences in energy-water coupling during the snow season. Such a framework would help focus research priorities as well as help the scientific community describe variability in how, where, and why streamflow will be impacted by changing snow dynamics.
Our data-driven review presents a consistent framework to answer a question fundamental to hydrologic sciences and water management: What are the mechanisms that give rise to variable streamflow response to changes in the amount, timing, accumulation, and ablation of winter snow cover? We first synthesize current understanding of the seasonal cycle of snow accumulation and ablation before summarizing prior work-and outstanding questionson how changes in this cycle under climate change are altering the timing, intra-annual distribution, and annual volume of streamflow. We then use our literature review to derive a simple framework centered around three potentially competing mechanisms that capture how abiotic and biotic factors differentially impact the interactions between energy and water during the snow season. These mechanisms include: (a) snow season water vapor losses, (b) the intensity of liquid water inputs (LWI), which we define as the amount of liquid water that reaches the ground surface from either rain or melting snow, and (c) the synchrony of water and energy availability. We include a number of demonstrative experiments that highlight which regions across the continental United States (CONUS) may be most sensitive to changes in a single snow metric (fs) before suggesting research opportunities.

Accumulation and ablation of seasonal snow cover
The annual cycle of snow accumulation and ablation is driven primarily by incoming solar radiation, secondarily by longwave radiation exchanges between the snowpack and atmosphere, and to a lesser degree by turbulent exchanges of sensible and latent heat and ground heat flux (Marks and Dozier 1992). Feedbacks between radiative and turbulent energy exchanges may either exacerbate or buffer the effects of climate change, resulting in spatially variable responses to the widespread warming observed across western North America (Harpold et al 2012, Bormann et al 2018, Harpold and Brooks 2018. To capture spatially and temporally variable responses in snow conditions, studies have proposed a number of snow metrics, which we review below.

Snow dynamics, terminology, and measurement
Point-scale metrics such as event snowfall, accumulated snow depth, and peak SWE have been broadly adopted to quantify the amount of snow on the ground (Bohr and Aguado 2001, Broxton et al 2016. Lapse rate models and assumptions regarding the relationship between temperature (T) and  (Lundquist et al 2008, Maggioni et al 2016, Skofronick-Jackson et al 2018 of snowfall in complex terrain is helping to distinguish between snow fall and redistribution or ablation of the snowpack. Contemporaneous advances in airborne and spacebased remote sensing have also enhanced the measurement of snow cover and SWE at large spatial scales (Tedesco et al 2015, JianCheng et al 2016, facilitating improved gridded snow products (e.g. Broxton et al 2016). In spite of these advances, accurate estimation of the amount, phase, and intensity of P in complex mountain environments remains a challenge (Rasmussen et al 2012) with evidence that high-resolution atmospheric models may be particularly useful tools for enabling further improvement (Lundquist et al 2019). In this vein, combinations of process-based snow modeling (Painter et al 2016) and observations are beginning to uncover terrain-mediated complexities in snow cover (Dong et al 2005). Together, these advances-particularly with respect to satellite-based sensors like the Moderate Resolution Imaging Spectroradiometer-have improved our ability to characterize the temporal and spatial patterns of snow dynamics across regional Historical challenges in obtaining direct measurements and developing scalable metrics for snowfall, snow fraction, or SWE led to the focus on choosing a fixed date (e.g. 1 April) to estimate net snow water input (Changnon et al 1991, Cayan 1996). However, Hamlet et al (2005 used trend analyses on 1 March, 1 April, and 1 May SWE to show how regions experience differential sensitivity to changes in P and T (e.g. higher sensitivity of SWE to P in cold regions and higher sensitivity of SWE to T in warm regions). Alternative, although less common, approaches have relied on the day of the water year (DoWY) of peak SWE to account for regional variance in snow accumulation and ablation (Bohr andAguado 2001, Girotto et al 2014). We summarize several of these metrics and provide a standardized definition to assist in their evaluation in table 1.

Snowpack energy balance
A warming climate interacts with seasonal snow cover through energy exchanges that occur between the snow surface and the atmosphere. We ground our discussion of the potential effects of these changes in the energy balance of a snowpack, written as: where ∆q is the change in energy [J m −2 ], t is time [s] and F is net energy flux [W m −2 ] and: Net radiative fluxes (net radiation) [R n ; W m −2 ] dominates snowpack energy balance and is composed of incoming solar (positive, towards the snow surface) and bi-directional longwave radiation (see equation (6)) (Marks and Dozier 1992). H is net sensible heat flux [W m −2 ] and maybe be positive (downwards towards the snow surface) or negative (upwards away from the snow surface) depending on snow and air T (Marks and Dozier 1992). LE is net latent heat flux [W m −2 ] and is typically negative (outgoing away from the snow surface) in continental snowpacks but may be positive in maritime environments. Both H and LE fluxes strongly relate to boundary layer turbulence and wind speed that drive air exchanges between snowpack and overlying atmosphere (Massman et al 1997, Lee et al 2004 where q cc [J m −2 ] is commonly known as the cold content and is the total energy required to raise the T s to 0 • C: And q melt [J m −2 ] is the energy associated with phase change: where c i is the heat capacity of ice [2102 J kg −1 K −1 ], T s is the average T of the snowpack, T m is the melting point of ice (0 • C), ρ w is the density of water (∼1000 kg m −3 ), and h s is the snow depth [m], h swe is the SWE [m], and γ is the latent heat of fusion Determining the response of snowpack to climate change is complicated by the fact that turbulent exchanges associated with T and LE are typically much small than radiative fluxes (Marks and Dozier 1992), with R n varying seasonally as a function of solar angle, day length, T, cloudiness, and spatially due to near surface topography, terrain, and vegetation structure. These complexities require a closer examination of snowpack radiative energy balance: where R s is incoming (positive) solar radiation During sunny days, R s -which is driven by solar angle (a function of day of year, latitude), aspect, topographic reflectance, and both topographic and vegetative shading-is positive and typically the largest energy flux in equation (6). The albedo (α) of fresh snow is high although it can be modified by snow grain size (typically related to time since last snowfall) and impurities in the snowpack (Skiles et al 2012, Deems et al 2013. Net longwave radiation (R l−in −R l−out ) is a function of the T of the snowpack, T of the overlying atmosphere, and differences in emissivity of snow and air (Marks and Dozier 1992). For example, the atmosphere has a lower emissivity than snow resulting in a cooling of the snowpack below ambient air T, especially at night. In contrast, clouds have similar emissivity as snow which may prevent snowpacks from cooling, especially at night (Ambach 1974, Plüss andOhmura 1997). The amount, extent, persistence, and freshness of snow strongly influence R n (e.g. Meira-Neto et al 2020) by altering the fraction of R s that is reflected (Schneider andDickinson 1974, Ingram et al 1989). Snow has much higher albedo than most terrestrial surfaces, and thus as R n increases the climate system reduces snow cover in a positive feedback process known as the snow-albedo feedback (Hall 2004, Déry and Brown 2007, Fletcher et al 2012, Qu and Hall 2014, Thackeray and Fletcher 2016. If early season snowfall is sufficient to cover the lower albedo terrestrial surfaces, the snow-albedo feedback will reduce R n and tend to preserve snow cover until solar angles increase in the spring (e.g. Koster et al 2010). A reduction in snow cover during spring when Rs is higher can enhance local warming (Hall 2004, Déry and Brown 2007, Thackeray and Fletcher 2016 with the future potential for a 1% reduction in surface albedo per degree of warming (Hall et al 2008, Fletcher et al 2012, Qu and Hall 2014. In contrast, spring snowfall events may increase albedo, reducing R n , cooling the local environment, and delaying melt. Superimposed on these larger scale radiative feedbacks are the effects of landscape heterogeneity, including aspect and vegetation (Harding and Pomeroy 1996, Broxton et al 2015, Tennant et al 2017, Safa et al 2021 that remain difficult to measure in complex, mountainous terrain (Reba et al 2009) and under variable snow cover conditions (Schlögl et al 2018). These feedbacks may result in high local variability in the partitioning of available energy to sublimation fluxes, cooling the snowpack, versus greater energy fluxes causing melt, advancing snowmelt, and causing snow-albedo feedbacks (Sexstone et al 2018).

Water vapor fluxes between atmosphere and seasonal snowpacks
Complex interactions between the snow surface and the atmosphere drive variability in the amount of water vapor lost during the snow season. Strong T lapse rates associated with either orographic or convective uplift can result in snowfall during most months in high mountains; however, consistent seasonal snow cover does not begin to accumulate until solar angles are low and R n favors net cooling of the land surface following a fresh snowfall (Bales et al 2006, Markovich et al 2019. During the snow season, air T is low and transpiration is typically assumed to be limited (Day et al 1989, Huxman et al 2003, Goulden and Bales 2014, Bowling et al 2018; however, other water vapor losses are possible. For example, exposed snowpacks above treeline, in meadows, forest gaps, fields, and in forest canopies are subject to considerable vapor loss over winter from sublimation (Rinehart et al 2008, Veatch et al 2009, Gustafson et al 2010, Harpold et al 2012, Biederman et al 2015, Sexstone et al 2018. In continental alpine systems, 20%-30% of winter snowfall may sublimate before melt (Hood et al 1999, Sexstone et al 2018. Sublimation effects are most dominant during snowpack accumulation and increase with low atmospheric pressure, low humidity, increased solar radiation and high wind speed (Lundberg and Halldin 2001, Earman et al 2006, Stigter et al 2018. Sublimation of snow intercepted by forest canopies are estimated at roughly 30% of local snowfall depending on leaf area (Storck et al 2002, Essery et al 2003, but similar to sublimation from open exposed snowpack on the ground, their sensitivity to climate change is poorly characterized (Lundquist et al 2021). Potential water vapor losses from the snowpack to the atmosphere may increase due to greater energy availability (both R n and T) (Meira-Neto et al 2020) and lower LE associated with evaporation of liquid water within the snowpack relative to sublimation of ice crystals in winter (Jambon-Puillet et al 2018).

Snowmelt and catchment liquid water input (LWI)
During the snow season, the intensity and timing of LWI-the amount of liquid from either rain or melting snow that reaches the ground surface in a given control volume (e.g. catchment)-is typically a function of snowmelt rates (Trujillo and Molotch 2014, Harpold and Kohler 2017, Musselman et al 2017, Harpold and Brooks 2018 and is thus relatively predictable (Harpold and Kohler 2017). However, the season-long interaction between P, snow accumulation and redistribution, and ablation processes causes highly heterogeneous snowmelt (Tennant et al 2017). Snowpacks become isothermal at 0 • C and begin to melt as solar angle increases, days lengthen, and surface albedo decreases (Cline 1997, Skiles et al 2012), which can be influenced by warming T, cloud cover, increased humidity (Clow 2010, Harpold and Brooks 2018), as well as snowpack depth and aspect (Kormos et al 2014, Christensen et al 2021. For snowmelt water infiltration into the soil zone sufficient melt must occur to overcome matric forces within the snowpack (including preferential flowpaths), which is often referred to as the snowpack becoming 'ripe' Woo 1984, Leroux andPomeroy 2017). As such, large snow-covered areas in the catchment which receive sufficient energy to overcome cold content and become ripe experience a relatively predictable seasonal increase in LWI driven by snowmelt (Dunne and Black 1971). When driven by snowmelt, LWI tend to occur over a longer duration than when driven by rain and at an intensity well below the infiltration capacity of most mountain soils, particularly those with well-developed organic layers (Liu et al 2008). In the absence of rare extremely warm rainfall and condensation (rain on snow) events, for example, snowmelt rates in the western CONUS rarely exceed 10 cm per day and average closer to 2.5 cm per day (Harpold and Kohler 2017). An important exception is frozen soils, where lower infiltration rates can be exceeded by LWI (Shanley andChalmers 1999, Bayard et al 2005). As a result, LWI driven primarily by snowmelt during the snow season are often associated with more consistent hydrological effects on shallow subsurface flow response, which is the dominant pathway for water redistribution and streamflow generation in upland catchments (Barthold and Woods 2015).
Because of moisture thresholds imposed by catchment-scale properties (e.g. infiltration capacity or soil water holding capacity), both the timing and intensity of LWI control how it is partitioned between subsurface and surface pathways (Harr 1981, Barnhart et al 2016, Berghuijs et al 2016. In general, snowmelt lags P inputs and there is some evidence suggesting that sequential melt results in more substantial soil moisture response-especially at deeper depths-than ephemeral snowmelt and rainfall (

Streamflow response to changing snow dynamics
Interactions between snowmelt-driven LWI and the subsurface, atmosphere, and vegetation exert a complex control over streamflow generation in mountainous catchments. For example, subsurface and surface lateral flow arising from seasonal increases in snowmelt-driven LWI leads to seasonal increases in streamflow generation via a number of mechanisms including infiltration excess overland flow (Horton 1933), saturation excess overland flow (Dunne 1978), preferential subsurface flow, as well as fill and spill flow in certain cases McDonnell 2006, McDonnell et al 2021). To help illustrate these connections, we adopt the below form of the snow season water budget for upland catchments following Godsey et al (2014), which we modify to include error: where  (Miralles et al 2020). We adopt ET in this paper due to its widespread use, however, we clarify that it refers to the loss of water to the atmosphere via evaporation, transpiration, and sublimation including canopy interception effects and blowing snow sublimation (Mcmahon et al 2013). The theoretical maximum rate of vaporization from a saturated surface is often referred to as potential evaporation or evapotranspiration (PET), which is a function of both the available energy (primarily R n , see also equations (2) and (6)) and the atmospheric water vapor pressure deficit (Meira-Neto et al 2020). Without limitations on available water, ET would be expected to equal PET. During the snow season, the equation (7) assumes that groundwater inflows and outflows from neighboring catchments are minimal although this assumption may be problematic in upland catchments over longer periods (Fan 2019). Interactions between variables in equation (7) lead to a distinct hydrograph in snowmelt-dominated systems typified by relative predictability in the timing and distribution of annual Q volume (Pagano et al 2004). Given the importance of stationarity for water supply prediction (Milly et al 2008), a diversity of metrics have been developed to better characterize the timing, distribution, and volume of Q. We detail several of these Q metrics below, including their derivation and significance for a discussion of the snowmelt-driven hydrograph and evidence for potential changes in Q.   Subsurface hydrological processes are typically invoked to explain catchment to regional scale streamflow timing sensitivity to changing snowpack. Grant (2009), Safeeq et al (2013), Maurer and Bowling (2014), and Harpold and Molotch (2015) all suggested that regional-scale subsurface hydrology provides a mechanistic explanation for the variable sensitivity of Q timing to climate change. Harpold and Molotch (2015), for example, emphasized that the timing of peak soil moisture can either exacerbate or moderate the sensitivity of Q to changes in snowmelt timing. Work in the Pacific Northwest highlights the role of bedrock properties and larger subsurface storage in moderating spring flows Grant 2009, Safeeq et al 2013). Additional synthesis efforts have highlighted T and/or P, elevation, and atmospheric circulation variations to explain differences in Q timing sensitivity with mid to lower elevation basins exhibiting the greatest potential for earlier Q in western North America (Stewart 2009).

Sensitivity of annual streamflow volume to changing snow dynamics
summary of spring Q timing. We note that spring timing is measured in a multitude of ways (e.g. runoff timing, melt-out, peak Q, DOQ25) in the studies reviewed. We elected to use spring Q or runoff timing as an umbrella term to reflect the different ways in which changes in Q timing is measured. * Denotes a study that presented variable results, but where some results outside the scope of this review and thus categorized as earlier. In the case of Jeton et al (1996), we excluded their cooler scenario results, Arnell (1999) conducted a global analysis and we included only results for western North America. We also note that some results presented evidence for stronger trends in certain regions, as is the case for earlier spring Q or runoff timing at mid-elevation basins. In these cases, we followed the authors description of their results in categorizing them as earlier versus variable to the best of our ability.

Snow effects on intra-annual streamflow distribution
Empirical work in the CONUS clearly connected advances in streamflow timing to an increase in the seasonal fraction of streamflow occurring during the snow season (Dettinger and Cayan 1995, Stewart et al 2005). Dettinger and Cayan (1995) and Stewart et al (2005) found statistically significant increases in winter Q and decreases in warm season summer Q, which were echoed in smaller-scale analysis by Nayak et al (2010) in Idaho. Using hydrological modeling, Godsey et al (2014) later showed that simulated future changes in fs lead to a 10% decrease in the volume of warm season Q with evidence of considerable intersite sensitivity to changes in fs. A hydrogeologic analysis by Safeeq et al (2014) suggested that future changes in fs might render areas with higher summer Q (greater subsurface storage) particularly vulnerable to climate change. On the whole, many of these lines of evidence about the intra-annual distribution of Q (Dettinger and Cayan 1995, Regonda et al 2005, Stewart et al 2005 emphasized that they did not translate into statistically significant changes in the amount of interannual Q. Consistent with this conclusion, some later research reported changes in the seasonal fraction of Q without attendant changes in the magnitude of annual Q (Nayak et al 2010). The legacy of these findings is one line of evidence suggesting that changes the intra-annual distribution of Q without necessarily impacting the volume of annual Q. . Some research has proposed that earlier streamflow timing and changes in the distribution of Q (e.g. increase in fractional Q during the snow season) increase the amount of runoff-and subsequently the annual volume of Qassuming that the timing of energy inputs remains the same (Tague and Peng 2013). Jeton et al (1996), for example, used a process-based model to suggest that increased asynchrony between water and energy inputs may increase Q from high elevations basins and decrease Q in middle elevation basins. However, other research (Risbey andEntekhabi 1996, Dettinger et al 2004) found that advances in Q may also increase water-limitations on ET, thereby offsetting impacts on Q. Reflecting this uncertainty, other have recorded scattered trends in both the mean and median annual flow (Stewart et  However, physical explanations for the hypothesis that declines in snowfall will drive declines in streamflow remain elusive (Berghuijs et al 2014). Some research has posited that increases in spring P may also buffer Q against declines in fs (Pederson et al 2011) or that rainfall and mixed P inputs during the winter may countervail reductions in Q from declining fs (Hammond and Kampf 2020) and model results suggest that higher snowmelt rates may have a larger effect on runoff ratios (Barnhart et al 2016). There is also recent evidence that efforts reliant on Budyko-based estimates of streamflow response to snow may be sensitive to poor assumptions about PET (Meira-Neto et al 2020). Complicating matters further, McCabe et al (2018) found that declines in fs have not altered runoff ratios in an empirical analysis of streamflow in the Pacific Northwest. These mixed findings on annual Q volume and runoff efficiency to changing snow conditions limit our ability to anticipate and respond to climate change.

Towards a framework linking snow processes and streamflow generation
Despite abundant research on changes in snowmeltdriven Q (section 3), we lack a robust, consistent, and readily testable framework to explain varying Q response to climate change. Below, we present a conceptual framework that distills interactions between snow and the atmosphere, vegetation, and the subsurface into three inter-related mechanisms (figure 2) that can be tested using different snow (e.g. fs as in our demonstrations) and Q metrics (e.g. annual volume and runoff efficiency):

(a) Mechanism 1-Snow Season Water Vapor
Fluxes. Snow dynamics influence the available energy via equations (1)-(6), which influence the timing and amount of water vapor fluxes to the atmosphere during the snow season (figure 2).

(b) Mechanism 2-Intensity of Liquid Water
Inputs. Because snow persists after it falls, snow dynamics can also modify the intensity of LWI, which in combination with site-specific thresholds (e.g. soil water holding capacity or infiltration capacity) can alter how water is partitioned to other water budget variables (figure 2). (c) Mechanism 3-Energy-Water Synchrony.
Snow enables the release of LWI after P has fallen. As such, snow dynamics facilitate greater temporal synchrony between LWI and PET during periods of higher radiation (figure 2).
We designed several data-driven demonstrations that rely on publicly available data from the Catchment Attributes for Large Sample studies (CAMELs) database (Addor et al 2017) to illustrate variable sensitivity to each mechanism across a range of study sites in the CONUS. Please see text S1 (available online at stacks.iop.org/ERL/17/053004/mmedia) for a full description of the data and figure S1 for a map of study sites. Importantly, the intention of these experiments is not to quantitatively relate our mechanisms to Q or examine the effects of other snow dynamics, but rather to establish consistent mechanisms that can be further developed and tested in future empirical and process-based modeling work. As such, we focus explicitly on demonstrating the maximum potential sensitivity of each mechanism to a decline in fs to zero (e.g. an all rain future) because it is well-connected to both snow and Q metrics (see section 4.1 below), reliably quantified without incorporating additional remotely-sensed data, and used to investigate Q response in several widely-cited studies (Berghuijs et al 2014, McCabe et al 2018). However, our framework leaves much room to be improved with additional snow and Q metrics to coordinate research on how and why Q response to changing snow dynamics varies. Additionally, our framework and demonstration are intended to establish a consistent set of mechanisms, thus we elect to treat each mechanism as distinct, but discuss the implications of interacting mechanisms in the conclusion.

Isolating snow metrics important for streamflow across the CONUS
Through correlation statistics, we assess relationships between snow metrics presented in table 1 and Q metrics in table 2 in 537 catchments to justify our  focus on fs. Figure 3(A) illustrates that snow metrics are highly correlated to each other, with few metrics exhibiting Spearman correlation values below ∼0.4. The mean annual fs captures a variety of snow dynamics similarly to snow persistence. Although each metric ultimately conveys different information about the characteristics of snow, both metrics exhibit the strongest and broadest correlation with other snow dynamics (figure 3(A)) and Q characteristics ( figure 3(B)) of interest. As expected from past studies, relationships between the fs and Q volume metrics were weaker than timing metrics and fs negatively correlated to flood metrics (Davenport et al 2020). All of the metrics used in this evaluation with the

) as illustrated by linear regression between snow season
Rn and fs as illustrated using Daymet data from the CAMELs database (please see text S1). (A)-(C) Grouped site-year regression based on average daily snow season incoming shortwave (Rs) for low, medium, and high P environments (n ∼ 5000 per low, medium, and high P). Grey bounds represent the 95% confidence interval for binned regressions. (D)-(F) Map of the maximum potential daily increase in snow season ET is calculated between average historical fs and fs = 0, then normalized by snow season P. exception of streamflow timing (e.g. DOQ 25,50,75 ), runoff ratio, baseflow index, and flood metrics rely heavily upon modeled data.

Mechanism 1: changes in snow season water vapor fluxes
Snow dynamics influence abiotic interactions between the land surface and the atmosphere by exerting a first-order control over the amount of available energy that can drive water vapor fluxes (Mechanism 1, figure 2). To investigate differences in the regional sensitivity to Mechanism 1, we performed several linear regressions between fs and R n grouped by snow season P amount (figures 4(A)-(C)). Grey bounds in figures 4(A)-(C) reflect the 95% confidence interval for regressions. We refer to the reader to figure S2 for scatterplots of the underlying data. Grouping by P accounts for the correlation between snow season fs and snow season length. Within snow season P groups, annual data were binned into nine groups of roughly equal number of catchments based on mean daily snow season solar radiation. We then used linear regression to identify the expected value of R n when fs = 0 (figures 4(D)-(F)) and estimated the maximum potential change in R n as the difference between modeled R n using the mean historical snow season fs and the modeled R n when snow season fs = 0. Assuming that the maximum potential change R n was exclusively available to latent heat flux per equation (2), we then converted this value to a water flux (i.e. ET using latent heat of vaporization into mm d −1 ) and normalized the resulting value by mean annual snow season P. We note that although the 95% confidence interval for our regressions is narrow in many cases (grey bounds in figures 4(A)-(C)), our demonstration relies on a regression model (please see figure S2 for more indepth presentation of the underlying data).
Linear regressions support the moderating role of higher snow season fs on R n in sunny, moderate to high P environments (orange and red lines in figures 4(B), (C) and (E), (F), dark green circles). We summarize these linear regressions in tables S1-S3. Consistent with our linear regressions, we observe strong coherence between circle color (indicating higher potential changes in the ratio of ET to P) and circle size (indicating historically higher fs) in Figure 5. Potential for changes in LWI intensity (Mechanism 2) as illustrated using Daymet data from the CAMELs database (please see text S1). Here we calculate the potential shift in LWI intensity (mm d −1 ) over 1 d, 3 d, and 14 d intervals as the difference between LWI and P. medium and high P environments (figures 4(E) and (F)). Results suggest that areas with historically larger snowfall during the snow season have the highest potential sensitivity to increases in snow season water vapor losses. These dynamics are most important in the interior western and central CONUS with some evidence of sensitivity at higher latitudes in the northeastern CONUS.

Mechanism 2: changes in LWI intensity
In certain environments, the unique characteristic of seasonal snow to persist after falling can facilitate the release of LWI later in time and at a lower or higher intensity than incoming P (Mechanism 2, figure 2). Potential for changes in LWI together with physical hydrological thresholds (e.g. soil water holding capacity) exerts a first-order control on infiltration and runoff generation. We investigated the maximum potential change in peak LWI intensity (mm d −1 ) under a complete transition from snowfall to rainfall (fs = 0) in figure 5. We selected the annual peak intensity of LWI and P inputs for each catchment for a running 1 d, 3 d, and 14 d mean value. Windows were selected to capture different potential effects of LWI intensity on streamflow generation. Assuming stationarity in the intensity of P, we then estimated the maximum potential change in the intensity of LWI as the difference between the intensity of snowmeltdriven LWI and P inputs when fs = 0 (no snow storage). Increases in rainfall intensity due to higher saturated vapor pressure with rising T (i.e. 7% increase per 1 • C of warming) (Trenberth 2011) or melt during rain-on-snow events (Li et al 2019) could further intensify LWI beyond what we consider here.
Our results indicate that maximum potential LWI intensity at 1, 3, and 14 d durations will increase from historical snowmelt values as snowfall turns to rain, especially in catchments with higher historical fs (size of the grey circles in figures 5(A)-(C)). This effect is particularly apparent for the 1 and 3 d LWI intensities (figures 5(A) and (B)), with substantially lower shifts in the 14 d LWI intensities after shifts to rainfall ( figure 5(C)). Intuitively, in places with historically low fs there is already relatively little difference in LWI and P intensity, which is reflected in the clustering of small grey circles at x = 0 in figures 5(A)-(C). The map in figures 5(D)-(F) highlights broad regional trends in sensitivity at the different temporal scales. The intensity of LWI in catchments at higher latitude and in the interior western CONUS appears most sensitive an increase as fs declines. However, the relationship between annual fs and maximum Figure 6. Potential for changes in the synchrony of water-energy inputs (Mechanism 1) as illustrated using Daymet data in the CAMELS database (please see text S1). (A), (B) Mean annual historical DoWY when 25% and 50% of LWI occur versus maximum shift in days between 25% and 50% of LWI timing and P timing. (C), (D) Map of the potential maximum shift in the timing of 25% and 50% of LWI when fs = 0 (e.g. all LWI driven by rainfall).
shift in LWI intensity does not fully explain the regional patterning in figures 5(D)-(F). For example, catchments along the central CONUS-Canadian border and some catchments along the western slope of the Appalachian Mountains also show large differences between LWI and P intensity. This suggests that other factors, such as intense spring rain, might explain or modulate the sensitivity to changes in LWI intensity.

Mechanism 3: water-energy synchrony
Snowpacks provide temporary storage of higher winter P that is released as LWI later in the season, which is a first-order control on Q via the partitioning of stored water to ET versus Q generation or groundwater (Mechanism 3, figure 2). We simulate the maximum potential difference in the timing of LWI and P inputs under a complete transition from snow to rain (fs = 0). For each catchment, we calculated the mean DoWY in which the catchment achieved 25% and 50% of LWI and P inputs across all years. Using these data, we then approximated the maximum shift in the timing of LWI inputs as the difference between 25% or 50% LWI and P inputs, respectively (i.e. assuming that the timing of LWI would equal the timing of P inputs if fs = 0 and there is no snowmelt to modulate the timing of LWI inputs). Consistent with Tague and Peng (2013) and our own analysis of the PET timing variability, we assumed that changes in the timing of LWI was the largest driver of potential asynchrony between LWI and energy inputs.
We observed that catchments with a higher annual fs also experienced greater maximum shifts between 25% and 50% LWI and P inputs as indicated by increase in circle size along the y-axis in figures 6(A) and (B). This relationship between fs and the maximum potential shift in the timing of LWI is evident across catchments regardless of the DoWY they reach 25% or 50% of their annual LWI (e.g. x-axis in figures 6(A) and (B)), although it does appear to scale the magnitude of the shift between LWI and P inputs. Intuitively, catchments with less annual snowfall experience relative harmony between LWI and P inputs. The translation of these results to geographic locations in figures 6(C) and (D) highlights distinct regional trends in sensitivity, with largest potential shifts in the timing of 25% and 50% LWI in the interior western CONUS (dark purple Figure 7. Summary graphic highlighting the potential utility of the framework proposed as part of this data-driven review. Here, we use both experimental results and literature review with experience of the authorship team to highlight where each mechanism-and the set of processes it represents-may be most important. Hydrologic regions correspond to United States Geological Survey (USGS) HUC (Hydrologic Unit Code) two boundaries. We show the study sites and HUC 2 boundaries in figure S1. circles). In both the 25% (early streamflow generation) and 50% LWI (peak streamflow generation) cases, the maximum risk for potential shifts in LWI timing were well connected to mean annual fs, as demonstrated by the coherence between circle size and shading in figures 6(C) and (D).

Summary and conclusions
Our data-driven review focuses on how regional variability in climate differentially influences the partitioning of winter precipitation along a gradient of fractional snow cover within the CONUS. We identify three inter-dependent mechanisms based on how snowpack mass and energy balance interacts with local climate, infiltration, and catchment water storage. Our framework leads to testable hypotheses useful for evaluating regional variability in streamflow response under a warmer climate.
• Mechanism 1 addresses how a warmer and drier climate will impact snow season water vapor fluxes. Although often ignored, snow season ET can be a large component of the annual water budget, which is likely to increase in a warming climate. These losses are likely to be greatest in the Great Basin, Missouri, Upper and Lower Colorado, Souris-Red-Rainey, and Arkansas White-Red basins (figure 7), with the potential to exacerbate summer drought stress and reduce annual Q consistent with Milly and Dunne (2020) who used a physically-based model in the Upper Colorado River Basin to estimate a 9.3% decrease in Q per degree Celsius of warming because of increased ET due primarily to snow albedo feedbacks. • Mechanism 2 addresses how increased LWI intensity driven by snow season rainfall (e.g. fs = 0) will interact with physical hydrological controls on infiltration and routing. The maximum potential risk for more intense LWI are greatest in the Pacific Northwest, California, Great Basin, Upper Colorado, and Missouri Basins in the western CONUS as well as the Great Lakes and New England Basins in the eastern CONUS (figure 7). These changes (figure 5) would be expected to increase the amount of Q occurring during the snow season consistent with Davenport et al (2020), who showed that declines in fs led to proportionally larger increases in streamflow and advance the timing of spring Q. As such, changes could potentially exacerbate summer drought stress per Harpold and Molotch (2015) with variable impacts on annual Q volume. • Mechanism 3 addresses the role of subsurface storage in buffering earlier water inputs (energywater synchrony) in a warmer, rainier climate. The maximum potential risk for decreased temporal synchrony between water and energy inputs is greatest in the Pacific Northwest, California, Great Basin, Upper Colorado, Rio Grande, and Missouri Basins in western CONUS ( figure 7). However, the extent to which temporal asynchrony may or may not impact seasonal and annual Q volume or drought stress remains difficult to parse consistent with Jeton et al (1996) who showed that higher and lower elevation basins experience bi-directional changes in annual Q volume in a warmer climate.
Our review provides a consistent framework for assessing the impacts of climate change on snowmeltderived streamflow, highlighting differential risks across regions of the western U.S. There are, however, limitations in our initial demonstration worth considering. First, using one snow metric (fs) as a proxy for climate change may not capture all the potential nuances in each of our mechanisms, especially at smaller spatial scales. Related to this, there is also a clear opportunity for future research to evaluate how the mechanisms we identify influence streamflow generation using the metrics aggregated for this review. Finally, although mechanisms may interact, an in-depth investigation of these interactions was beyond the scope of our demonstration. However, our initial results can help to characterize end-members for future research by establishing a set of testable, inter-related mechanisms that reflect the dominant processes connecting snow to streamflow across CONUS. Snowmelt-driven streamflow in areas with either less persistent snow cover, small historical fs, and with low intensity P during the snow season (e.g. the HUCs in the northeastern CONUS in figure 7) may be easier to predict than in areas with persistent snow cover, large historical fs, and higher intensity P during the snow season (e.g. HUCs in the interior western CONUS in figure 7), where snowmelt-driven streamflow is likely to have unique feedbacks not completely captured in our framework.
Ultimately, hydrologic models need to capture potential interactions of these three mechanisms to accurately predict future changes in snowmelt-driven streamflow. Specifically, improvements in modeling capabilities should be focused in three areas: (a) representing complexities in snowpack-atmosphere energy fluxes and how this variability influences snow ablation; (b) representing variability in subsurface storage as soil and groundwater and how these stores are partitioned to either atmosphere fluxes or streamflow; and (c) continuing to develop high quality forcing datasets, including P, T, R n , humidity, and wind speed, in complex terrain. Progress on these fronts requires advances in hydrologic models, process understanding, conceptual frameworks, and observations. As a critical link in this chain, our mechanistic framework offers value for evaluating and communicating changes in critical mountain water supplies in an increasingly complex and uncertain hydrologic future.