Interactions between snow cover and evaporation lead to higher sensitivity of streamflow to temperature

Estimates of potential evaporation often neglect the effects of snow cover on evaporation process. Here, we present a definition of potential evaporation that explicitly accounts for landscapes that are partially covered by snow. We show that, in the presence of snowpack, our evaporation estimates differ from conventional methods that assume evaporation from a free water surface. Specifically, we find that conventional methods overestimate potential evaporation as well as aridity, taken as the ratio of atmospheric water demand to supply, in landscapes where snowfall is significant. With dwindling snow-cover, actual aridity increases, which could explain the reduction in streamflow with decreasing snowfall. We suggest that streamflow, and hence water availability, is more sensitive to temperature changes in colder than in warmer regions. Colder, snow-dominated regions are more likely to experience greater changes in water availability with warming, suggests an analysis that explicitly includes snow cover in potential evaporation estimates.

F reshwater from snowmelt accounts for almost one-sixth of the global human consumption and agricultural use of water 1 , whereas the economic and societal impacts of snowcover loss can surmount trillions of dollars 2,3 . Existing trends in the reduction of snowfall 4,5 and snow coverage [6][7][8][9] are expected to continue under the projected global warming 10 .
Potential evaporation or evapotranspiration (PET) refers to the maximum rate of water vaporization from a saturated (or unstressed) surface to an unsaturated atmosphere mainly driven by the atmospheric water vapor pressure deficit (D) and the available energy (A E ; primarily net radiation) absorbed by the surface 11 . Operationally, it is estimated as the evaporation occurring from a free water surface [11][12][13][14] or the evapotranspiration of a well-watered short grass 12,14,15 with fixed surface albedo (0.08 for water and 0.23 for a reference crop) under existing atmospheric conditions. PET is widely used to assess continental drying through aridity indices 11,16 and drought indices 17 , whereas it is also used to estimate the atmospheric demand of water in hydrologic models [18][19][20] and offline assessments of climate change impacts on water resources 21 . However, the lack of a definition for PET over snow-covered surfaces has hampered progress towards understanding the fate of water resources for colder, snow-dominated regions, which are commonly excluded from global studies 11,21 . Under such conditions, the assumptions of a free water surface or the reference crop as the choice of evaporative surface may lead to unrealistic estimates of PET.
Here we define PET from a landscape partially covered by snow as the maximum water vaporization (evaporation and/or sublimation) from an unstressed water surface covered by an unstressed snow surface, of which the coverage is changing seasonally. The term "unstressed" refers to a zero surface resistance to molecular diffusion of water vapor through the roughness sublayer to the turbulence-dominated surface layer. To account for the dynamically changing surface albedo caused by seasonal snow accumulation/ablation, we used a physically based land surface model (LSM) driven by observed atmospheric forcing to estimate PET. We report the results from this alternative method for computing PET and compare these results to estimates of PET following a conventional definition that assumes the evaporation of a free water surface 13,14 . We use both methods to compute daily values of PET for 671 catchments from the Catchment Attributes and Meteorology for Large-Sample Studies (CAMELS) 22 located within the conterminous United States ( Supplementary Fig. 1a), over a period of 30 hydrologic years . We explore some immediate implications of the proposed changes by comparing 30-year climatology of both estimates as a function of average catchment temperature ( T) and also compute the aridity index (ϕ), a variable commonly utilized to assess global patterns of water availability accounting for both atmospheric water demand (PET) and supply (Precipitation). Following that, we offer an explanation of the existing paradigm stating that a reduction in snowfall fraction leads to a decrease in the mean streamflow at the catchment scale 23 , based on a framework 24 that assumes water balance to be controlled by ϕ. Finally, we demonstrate how temperature increases might lead to much higher reductions in streamflow in colder regions as compared to warmer regions.

Results
Estimation of PET. Under prolonged sub-freezing temperatures, water will not exist in its liquid form but will be frozen or exist as accumulated snow over the land surface, leading to significant changes in the surface albedo of the evaporative surface. For most catchments, snow is the dominant form of frozen water and the presence of snow will greatly increase the amount of reflected solar radiation and therefore reduce the net radiation 25 . Although a variety of analytical equations exist for estimation of PET for liquid water surfaces 14 , there are no simple analytical schemes that account for the vaporization from a surface subject to varying snow cover. For this reason, we used the physically based Noah-MP LSM 26 ("Methods") to compute PET over two scenarios of unstressed surfaces as follows: (1) a liquid water surface (PET water ) and (2) liquid water surface with seasonal snow cover (PET snow ). We ran the model for 30 years from 1983 to 2013 driven by observed hourly near-surface atmospheric data. We then computed monthly PET climatology (the average of the monthly model outputs over the 31 years) and PET climatology (the average of all data over the 30 years). Figure 1a shows monthly PET values (in mm equivalents) for a catchment with significant amount of snowfall (f s = 0.55), located in the state of Montana, USA. The monthly climatology of PET water (mm) is greater than that of PET snow in snow seasons but converges in summer. Over all CAMELS catchments, the difference in the 31-year PET climatology becomes more apparent for colder catchments with T ranging from −5°C to 10°C (Fig. 1b). PET water becomes much larger than PET snow by up to 50% in colder catchments with T around 0°C (Fig. 1c) and snow-dominated catchments with f s above 0.7 (see "Methods") ( Supplementary Fig. 2), suggesting that the conventional approach overestimates PET, progressively increasing towards colder, snow-dominated catchments.
The sensitivity of PET to increasing temperatures. To answer the question of how PET will change over snow-dominated catchments under a warming climate, we performed additional model experiments to compute PET water and PET snow with +1.0, +2.0, and +3.0°C increases in surface air temperature. To discern the physical mechanisms that may contribute to the warming-induced changes in PET, we used an attribution analysis 21,27 based on the well-established Penman 28 equation using monthly North American Land Data Assimilation System project phase-2 (NLDAS-2) atmospheric data and model outputs. Under the +3.0°C warming in South Crow Creek, Montana (USA), warming-induced increase in the water vapor pressure deficit (ΔD) appears to be the most dominant contributor to the increases in PET estimated by the conventional method (Δ PET water ) followed by changes in the slope of water vapor pressure vs temperature (ΔS) and available energy (ΔA E ), which reduce the positive contributions by ΔD (Fig. 2a). The negative effect of ΔA E on Δ PET water are caused by the energy loss through increases in the outgoing longwave radiation due to the surface warming (increases in the downward longwave radiation under a warming climate is not considered in the warming experiment). When the snow albedo effects are considered, ΔA E becomes a positive contributor to the increase in PET (Δ PET snow , Fig. 2b) due to warming-induced decreases in snow coverage and thus surface albedo. Between the two cases ( Fig. 2a, b), the largest change in these contributors is ΔA E followed by changes in aerodynamic factor (Δf u ). Warming-induced changes in f u (see Eq. (2) in "Methods") are due mainly to the changes in the surface drag coefficient, C M , for both cases. For Δ PET water , Δf u is negligible, whereas for Δ PET snow , Δf u becomes a positive contributor, because the warming may induce increases in surface roughness length, z 0, which are weighted by fractional snow cover (SCF), due to reduced snow cover (z 0 for water is 0.01 m, whereas 0.002 for snow surface). Over all the CAMELS catchments, warminginduced changes in PET are generally greater for partially snowcovered surfaces than for water due mainly to the warminginduced increases in available energy as a result of reduced surface albedo or surface darkening (Fig. 2c). The warming effects progressively increase as the extent of warming increases from +1.0°C to +3.0°C. The catchments with T around 5°C appears to be the most sensitive catchments, where the duration of snow coverage, changing more dynamically with seasons, is more sensitive to the warming.
Revisiting the snowfall-streamflow paradigm. A decrease in the fraction of precipitation falling as snow has recently been suggested to lead to a decrease in the mean streamflow at the catchment scale 23 when long-term catchment fluxes are analyzed through the widely adopted Budyko 24 framework. The Budyko 24 framework assumes that, at sufficiently long timescales, the mean annual precipitation is partitioned into evapotranspiration (E) and streamflow (Q), and this partitioning is controlled by the aridity index (ϕ), which represents the competition between atmospheric demand and supply of water. ϕ is commonly estimated as the ratio between the mean annual values of PET and mean annual precipitation (P). The basis for studying catchmentscale mean annual water balance within the Budyko framework is the definition of PET and ϕ, as their estimation determines the location of catchments along the x-axis in the Budyko space (ϕ vs. E/P). The snowfall-streamflow paradigm (SSP) 23 suggests that catchments with a greater snowfall fraction show a lower value of E/P than that predicted by the Budyko framework, suggesting that catchments with a greater snowfall fraction can produce more runoff, Q/P. Therefore, a warming-induced shift snow to rain would lead to a decrease in runoff. The PET estimates used in ref. 23 do not account for the presence of snow, as PET was estimated based on pan evaporation measurements 29 , from the months of May through October, whereas the estimates for the remaining months were extrapolated through a non-physically based procedure 30 . A physically based explanation for the SSP has been recently suggested for a study 31 within the context of the Colorado River Basin, where the reductions of streamflow were associated to an increase in evapotranspiration due to reduction in snow cover. However, such conclusions were drawn through a more expensive Monte Carlo simulation with a radiation-aware hydrologic model. Here we attempt to generalize the previous findings of the decreasing streamflow under a warming climate 23,31 through the lenses of our snow-corrected PET estimation.
We analyzed the positions of the catchments in the Budyko space to assess whether there are significant differences between the normalized streamflow (Q/P = 1 − E/P) of snow-dominated catchments and others with significantly less or no snowfall. For this analysis we have reduced the sample size to 214 catchments ( Supplementary Fig. 1b), after a careful quality-control assessment ("Methods"). We first assessed whether f s is associated with anomalies in Q/P when using PET water and PET snow to compute ϕ water and ϕ snow , respectively ("Methods"). Figure 3a, b show the positions of the CAMELS catchments based on ϕ water and ϕ snow ("Methods"), in which catchments are differentiated with respect to snowfall fraction, where a value of f s = 0.5 was chosen to highlight the differences in the catchment location within the ϕ − E/P space. Although an outlier behavior is suggested for snowdominated catchments when ϕ water is used (Fig. 3a), the same pattern cannot be seen in the case of ϕ snow (Fig. 3b). We then assessed how Q/P changes with varying degrees of f s (Δ Q/P , "Methods"). A significant positive relationship (r = 0.3, p = 10 −6 ) between Δ Q/P and f s exists when ϕ water is used, in which a 20% increase in Q/P occurs for each incremental increase in f s (Fig. 3c), consistent with the SSP as shown in ref. 23 . However, when ϕ snow is used, the trends are not significant (r = 0.02, p = 0.6) and negligible for Δ Q/P , suggesting the normalized runoff of snow-free and snow-dominated catchments to be approximately the same for catchments having similar values of aridity index. This also suggests the Budyko 24 framework to be valid for both groups, allowing the SSP to be properly hypothesized through the Budyko 24 framework as follows: In a warming climate, the increasing temperatures leading to a shift from snow to rain (i.e., decreasing f s ) will increase the atmospheric demand for water (i.e., the increasing PET and ϕ), leading to a reduction in streamflow at the mean annual timescale.
To further explore the interpretation of streamflow changes with the extent of warming for catchments with different snowfall extents, we estimated the changes in the normalized streamflow Δð Q P Þ caused by changes in PET (and therefore ϕ) computed for the different warming scenarios ("Methods") by assuming both PET water and PET snow leading to Δð Q P Þ water and Δð Q P Þ snow . Under the same +3.0°C warming scenario, colder catchments appear to lose more water and thus result in a larger drop in the normalized streamflow (Δð Q P Þ) when using PET snow than using PET water , with the difference being inversely proportional to catchment-averaged temperature (Fig. 4). Our results point to a much higher sensitivity of colder, snow-dominated regions than for warmer snow-free regions. For example, when comparing Δð Q P Þ snow of catchments within 0 ± 2°C average temperature with catchments within T ≈ 20 ± 2°C, this difference is of almost four times (Δð Q P Þ snow % 0:11, T ≈ 0°C vs. Δð Q P Þ snow 0:03, T ≈ 20°C), with similar results for other warming scenario ( Supplementary Fig. 3). This is explained by the higher positive changes in mean annual PET and therefore increased aridity experienced by colder regions due to reduction in snow-driven atmospheric demand as previously explained.

Discussion
The modeling assumptions adopted here represent an attempt towards a conceptual and physically meaningful estimation of PET for snow conditions, and we believe that additional investigations might be useful for its proper advancement. Also, the warming experiments performed with Noah-MP in this study do not include potential changes in other climatic variables that may contribute to changes in PET 10 , for instance, increases in the downward longwave radiation and specific humidity under a warming climate as consistently projected by many earth system models. Nonetheless, the prescribed warming experiments performed in this study provide a first-order approach towards the understanding of how PET and ϕ may respond to future climates. Our results also carry some important uncertainties arising from the limitations embedded in our water-balance analysis. The Budyko framework 24 is based on the simple climate-based reasoning on the long-term water balance, lacking detailed description of physical processes controlling the partitioning of precipitation into runoff and evapotranspiration. Moreover, the lack of a robust way for its application when predicting changes occurring in time vs. its ability to explain spatial patterns of water-balance partitioning 32 , as well potential issues arising from the water-balance closure assumption 33 should be carefully considered if a similar analysis is intended for prediction purposes. We conclude that the conventional approaches, e.g., the open water PET, overestimate PET and thus aridity index in snow-dominated catchments. We have shown that colder regions might experience a larger decrease in streamflow than warmer regions due mainly to increases in net radiation caused by decreases in surface albedo or surface "darkening." This result is consistent with a recent study, suggesting "warming-driven loss of reflective snow energizes evaporation" through a Monte Carlo simulation with a radiation-aware hydrologic model 31 and provides an explanation for the suggested catchment-scale reduction of streamflow occurring from a shift from snowfall to rainfall at the mean annual timescale 23 . Although current global-scale assessments of future changes in water availability performed through the lens of the Budyko framework do not include cold regions in their analysis 11,21 , we believe this study provide a promising step towards a generalized understanding of how water availability of colder regions is affected by increasing global temperatures. Finally, our results confirm recent concerns [1][2][3]23,31 and reinforce the call for greater attention towards water supply policies in snow affected environments.

Methods
Estimation of PET and f s . To estimate PET over unstressed surfaces of free water and partially snow covered, we employed a widely used LSM, Noah with multi-physics (Noah-MP) 26 . Noah-MP is a widely used LSM with demonstrated ability in simulating land surface fluxes 26,34,35 . Noah-MP considers detailed snowpack physics (e.g., liquid water retention and refreezing, densification, heat released or absorbed due to phase changes, and surface albedo changes due to snow aging), and surface energy and mass exchanges with the atmosphere (e.g., sublimation/deposition, patchy snow effects on surface albedo, and surface turbulent heat and moisture exchanges, etc.). Noah-MP explicitly represents sublimation/deposition by converting the latent heat flux over subfreezing surfaces by the latent of sublimation, which is the sum of fusion (334 kJ kg −1 ) and vaporization (2265 kJ kg −1 ). The surface albedo and roughness length of a water surface partially covered by snow is the areal average weighted by fractional snow cover, SCF, which is diagnosed from the model predicted snow depth using a snow depletion curve. Noah-MP was driven by the hourly North America Land Data Assimilation System Phase-2 near surface atmospheric forcing data over the contiguous US (CONUS) from 1980 to 2015 at 1/8th degree resolution 36 . We assumed two distinct scenarios, both of which assumed "unstressed" surfaces with a zero surface resistance to Fig. 3 The Budyko framework analysis for the selected 214 CAMELS catchments. Catchment location using a ϕ water and c ϕ snow [magenta dots represents catchments with f s > 0.5]. The relationship between the deviation of the observed Q/P (circles) from a fitted Budyko curve (red lines in a and c) and snowfall fraction, f s using ϕ water (b) and ϕ snow (d). T) for the selected 214 CAMELS under the "water" scenario, whereas red dashed lines show the changes according to the "snow" scenario. Shaded regions of both curves represent the SD taken from the groups of catchments within the same temperature ranges.
water diffusion in the roughness sublayer. Scenario 1 (water) was based on a liquid free water surface (i.e., assuming surface albedo is free water albedo = 0.08), whereas Scenario 2 (snow) allows snowpack to be built on "water" surface. For Scenario 1, we set up the snowfall fraction of the precipitation equal to 0, whereas no snow cover at the water surface was allowed by setting up SCF of the surface equal to 0. The water body temperature is solved at four layers with thicknesses of 0.1, 0.3, 0.6, and 1.0 m (totally 2 m) with a zero-heat-flux lower boundary condition at 8 m. The same as for soil. Noah-MP does not compute the frozen water depth but determines if the water surface is frozen from the resulting water surface temperature (T g ) (frozen if T g ≤ 0°C), where T g is computed by iterating the surface energy balance equation 26 . Noah-MP assumes that snowpack can build up over water surfaces when the surface temperature is below the freezing point (T frz ), while neglecting energy consumed/released during phase change. In addition, to ensure a conceptual agreement between Noah-MP evaporation and that of conventionally employed analytical solutions 11,13,28 , we excluded the stability correction to the surface drag coefficient, though the atmospheric warming may substantially increase the surface layer stability over snow and lake surfaces 27 suppressing the surface turbulent exchanges of water and energy with the atmosphere. We estimated f s based on a new snow-rain partitioning scheme using a sigmoidal function of wet-bulb temperature, which was demonstrated to improve snow simulations in the CONUS 35 .
Attribution analysis of the changes in PET. We adopted an "offline" attribution analysis 21,27 to elucidate the mechanisms contributing to the changes in PET with increasing temperature for both PET scenarios. The temperature increments adopted here were uniformly applied to the historical timeseries. The attribution analysis consists of using the monthly NLDAS atmospheric forcing data 36 and outputs from Noah-MP (net radiation and ground heat flux) to compute PET values through the commonly used Penman 28 equation, in which PET (mm) can be written as: where, s is the slope of saturated water vapor pressure vs. air temperature (Pa K −1 ), γ is the psychrometric constant (Pa K −1 ), A E is available energy (W m −2 ), a residual of net radiation (Rn, in w m −2 ) minus ground heat flux (G, in w m −2 ), f u is the wind function, D is the 2 m water vapor pressure deficit (Pa), and L V is the latent heat of vaporization (J kg −1 ). The wind function encapsulates the aerodynamic effects on evaporation and is defined as: In which ρ is the air density (kg m −3 ), P is the 2 m air pressure (Pa), C M is the unitless surface drag coefficient (or surface exchange coefficient), and u 2 is the 2 m windspeed (m s −1 ). C M controls the turbulence transport of momentum through the surface layer and is parameterized as a function of the surface roughness length, zero-displacement height (lifting up the evaporative surface), and surface layer stability through the Monin-Obukhov similarity theory 26 .
We then approximate the changes in PET (ΔPET) to the first order as: where the terms on the right-hand side represent the individual contribution to the total change ΔPET. These terms can be written as: ∂PET The offline ΔPET values obtained using Eqs. (3) through (9) showed a good agreement with ΔPET values from Noah-MP ( Supplementary Fig. 4), which allows the use of a simpler analytical equation to assess how different drivers might be responsible for changes in PET.
Catchment selection. We used daily streamflow from the CAMELS dataset 22 , which includes catchments from the contiguous United States (CONUS). We specifically looked at the hydrologic years 1 October through 30 September) within the span of 31 years (1983-2013). A few quality-assessment steps where undertaken for selection of catchments for the Budyko framework analysis (Figs. 3 and 4) as follows: (i) we eliminated all catchments with missing streamflow; (ii) the CAMELS dataset provides two different estimates of catchment area (km 2 ). As all fluxes used in here are represented in mm, the estimation of catchment area plays an important role in the flux magnitudes. Therefore, we removed catchments based on the criteria of 1% difference between the two estimates. (iii) Catchments having negative evaporative fraction ( E P < 0), or evaporative fractions > 1 ( E P > 1), were eliminated, to accommodate the assumptions of long-term water balance comprising the Budyko framework and (iv) catchments with areas smaller than 200 km 2 were removed to account for the inaccuracies arising from the transposition of NOAH-MP (pixel size ≈ 170 km 2 ) results onto the basin masks. (v) Catchments with significant surface water bodies (fraction of water of a catchments area >5%) were also removed to reduce the effects of storage in the analysis. The resulting procedure returner in 214 catchments ( Supplementary Fig. S1b) used for the analysis of Figs. 3 and 4.
Aridity Index (ϕ) computation. ϕ is conventionally computed as the ratio between mean annual values of PET (water demand) and P (water supply), in mm equivalents per year. When computed as energy equivalents, ϕ can be regarded as a ratio between the mean annual energy associated with the evaporation at the potential rate (PET ε , in MJ m −2 year −1 ) and the energy required to evaporate all the mean annual precipitation (P ε , in MJ m −2 year −1 ). In this way, P ε allows for the inclusion of the energy required for sublimation (or melt and evaporation) of the solid snow particles, which would not be possible when assuming mm equivalents. Therefore, ϕ was estimated as: where x denotes either scenario "water", leading to PET ε water (MJ m −2 year −1 ) or scenario "snow" leading to PET ε snow (MJ m −2 year −1 ), P is the mean annual precipitation (m year −1 ), λ v is the latent head of vaporization of water (λ v = 2.501 MJ kg −1 ), ρ w the density of water (10 3 kg m −3 ), and λ s is the latent heat of sublimation (λ s = 2.835 MJ m −2 ), being equal to the sum of λ v and the latent heat of fusion (λ f = 0.334 MJ m −2 ), Streamflow comparison. We tested whether there are significant differences between normalized streamflow (Q/P of catchments with and without significant amounts of snowfall across the whole range of observed ϕ for the CAMELS catchments according to the following: first, we fitted a Budyko-type curve (seen in Fig. 3) to the observed values of ϕ and E P : Equation (2) was fitted based on a sub-sample with little influence of snowfall (f s < 0.2), which showed no significant difference between ϕ 1 and ϕ 2 ( Supplementary Fig. 5c, d). The choice for a fitted curve rather than the original Budyko 24 equation was motivated by the bias in estimated vs. observed Q/P obtained by using the latter equation (Supplementary Fig. 5a, b). Following that, we calculated the difference between the predicted and observed values of Q P ðΔ Q=P Þ for all catchments and computed a linear regression (Fig. 3) between the deviation from the curve vs. f s . The results shown in Fig. 3 were also computed using the conventional formulation of the aridity index (ϕ ¼ PET mm ð Þ P mm ð Þ , i.e., assuming no differentiation between snow and rain in the denominator) and did not show significant differences (Supplementary Fig. 6).
Sensitivity of streamflow to increasing temperatures. Values of streamflow sensitivity to changes in temperature were computed as the changes in normalized streamflow (Δ Q/P ) due to changes in temperature through the partial differentiation: where the subscript x stands for scenarios of "water" or "snow," and F represents Q/P obtained through Eq. (11). The term ΔPET x represents the change in PET, obtained from the differences between Noah-MP runs with increments of +1.0, +2.0, and +3.0°C, and a baseline run without warming. The term ΔP ε represents the change in the energy required to evaporate/sublimate the mean annual precipitation and was obtained as the difference in P ε [the denominator of Eq. (10)] between the baseline and warming scenarios.