Streamflow generation and deep groundwater recharge may be vulnerable to
loss of snow, making it important to quantify how snowmelt is partitioned
between soil storage, deep drainage, evapotranspiration, and runoff. Based
on previous findings, we hypothesize that snowmelt produces greater
streamflow and deep drainage than rainfall and that this effect is greatest
in dry climates. To test this hypothesis we examine how snowmelt and
rainfall partitioning vary with climate and soil properties using a
physically based variably saturated subsurface flow model, HYDRUS-1D. We
developed model experiments using observed climate from mountain regions and artificial climate inputs that convert all precipitation to rain, and then
evaluated how climate variability affects partitioning in soils with
different hydraulic properties and depths. Results indicate that event-scale runoff is higher for snowmelt than for rainfall due to higher antecedent moisture and input rates in both wet and dry climates. Annual runoff also increases with snowmelt fraction, whereas deep drainage is not correlated with snowmelt fraction. Deep drainage is less affected by changes from snowmelt to rainfall because it is controlled by deep soil moisture changes over longer timescales. Soil texture modifies daily wetting and drying patterns but has limited effect on annual water budget partitioning, whereas increases in soil depth lead to lower runoff and greater deep drainage. Overall these results indicate that runoff may be substantially reduced with seasonal snowpack decline in all climates, whereas the effects of snowpack decline on deep drainage are less consistent. These mechanisms help explain recent observations of streamflow sensitivity to changing snowpack and highlight the importance of developing strategies to plan for changes in water budgets in areas most at risk for shifts from snow to rain.

Changes in precipitation phase from snow to rain can modify hydrological
partitioning by altering the timing and rate of inputs. Daily snowmelt rates
typically do not reach the extreme intensities of rainfall (Yan et al.,
2018), meaning that most areas (i.e., the Cascades) are predicted to receive
more intense water inputs with more winter rainfall, whereas some other
areas (i.e., Southern Rockies) will likely experience a decline in input
intensity with snow loss (Harpold and Kohler, 2017). Warmer areas like the
maritime western US may experience near complete loss of snowpack as snow
fully transitions to rain by the end of the 21st century (Klos et al.,
2014). Unlike rainfall, which is typically episodic, snow can accumulate
over time and melt as a concentrated pulse of soil water input (Loik et al.,
2004). This means that at 7 to 30 d scales snowmelt inputs are of greater
magnitude than rainfall (Harpold and Kohler, 2017). Concentrated snowmelt
can lead to a large proportion of runoff and deep drainage (Earman et al.,
2006; Berghuijs et al., 2014; Li et al., 2017). With climate warming, future
snowmelt rates may be reduced in many areas because earlier melt occurs when
solar radiation is lower (Musselman et al., 2017). Along with warmer
temperatures, increasing atmospheric humidity is leading to more frequent
mid-winter melt events in humid regions yet increased snowpack sublimation
and/or evaporation in dry regions (Harpold and Brooks, 2018). Given the
considerable heterogeneity in climate, soils, topography, and vegetation
across mountain ranges, the water budgets of different locations respond
unevenly to a loss of snow.

Both surface runoff and deep drainage are affected by soil texture, soil
depth, rooting depth (Cho and Olivera, 2009; Seyfried et al., 2005), and
topography. These properties have limited variability over time spans of
hydrologic analysis and can produce temporally stable spatial patterns of
soil moisture, where some parts of the landscape are consistently wetter
than others (Williams et al., 2009; Kaiser and McGlynn, 2018). Aspect
modifies the snowpack energy balance, leading to higher sustained soil
moisture content on north-facing slopes compared to south-facing slopes with
the same input (in the Northern Hemisphere; Williams et al., 2009; Hinckley
et al., 2014; Webb et al., 2015, 2018). Landscape evolution may
lead to deeper profiles and more deeply weathered rock due to wetter
conditions on north-facing slopes, making these slopes more conducive to
deep drainage in some locations (Hinckley et al., 2014; Langston et al.,
2015). Where soils are shallow, winter precipitation may exceed the
soil storage capacity, leading to both runoff generation and deep drainage
(Smith et al., 2011). Deeper soil profiles have greater storage capacity,
which can sustain streamflow, even with snow loss; however, given
consecutive years of low input these profiles will be depleted of storage,
leading to lower flows (Markovich et al., 2016). Deeper soils can also help
sustain transpiration during the late spring and summer, when shallow soils
have dried (Foster et al., 2016; Jepsen et al., 2016). Streamflow can be
insensitive to inputs under dry conditions, but respond rapidly once a
threshold soil moisture storage value is exceeded (McNamara et al., 2005;
Liu et al., 2008; Seyfried et al., 2009). McNamara et al. (2005)
hypothesized that when dry-soil barriers are breached, there is sudden
connection to upslope soils, leading to delivery of water to areas that were
previously disconnected. In their semi-arid study area, such breaching of
dry-soil barriers was only observed for periods of concentrated and
sustained input from high-magnitude spring snowmelt. Together, the complex
interactions of soil properties, antecedent conditions, water inputs, and
evaporative demand make it challenging to determine how changes from snow to
rain affect hydrologic response even in idealized settings.

Figure 1Conceptual illustration of study hypotheses indicating the
importance of concentrated snowmelt input (c, d) vs. intermittent input (a, b) for runoff generation. The wet climate (b, d) generates more runoff (Q) and deep drainage (D) and less
evapotranspiration (ET) compared to the dry climate (a, c). In
both climates, concentrated input can increase both Q and D because it is more likely to allow soil saturation than intermittent input, which allows ET
during periods of drying. The concentrated input from snowmelt leads to
greater increases in Q and D in the dry climate than in the wet climate because snowmelt is the most likely cause of soil saturation in dry climates.

The goal of this study is to improve our understanding of how changes in
precipitation phase from snow to rain affect hydrological partitioning in a
one-dimensional (1-D) representation of the critical zone. Partitioning of
precipitation input, P, can be into runoff, Q, defined as lateral export of water from the domain, evaporation, E, transpiration, T, deep drainage below the root zone, D, and storage within the soil root zone, ΔS.
Throughout this study, the term runoff refers to non-infiltrated input that
exits the domain laterally due to infiltration or saturation excess
mechanisms. Over a given time increment, partitioning can be tracked using
the water balance (Eq. 1).

(1)P=Q+E+T+D+ΔS

We address the following questions: (1) are snowmelt and rain partitioned differently
between Q, ET, and D? And (2) how is snowmelt and rain partitioning affected by climate, soil type, and soil depth? We use a physically based 1-D modeling approach to address these questions and systematically test hypotheses about hydrologic response to snow loss. The 1-D modeling approach allows for isolated comparison of climatic and edaphic factors on input partitioning; it is a simplified approach that neglects lateral redistribution of water.

We hypothesize that reducing the fraction of precipitation falling as snow
leads to lower runoff and deep drainage because it reduces the concentration
of input in time (Fig. 1). Concentrated input during melt of a seasonal
snowpack often saturates soils, causing saturation excess runoff and deep
drainage below the root zone (Hunsaker et al., 2012; Kampf et al., 2015;
Webb et al., 2015; Barnhart et al., 2016). Diffuse input over time reduces
the likelihood of saturation because it allows more water redistribution and
evapotranspiration between inputs. We also hypothesize that snowmelt is
critical for runoff generation and deep drainage in dry climates and deep
soils, where snowmelt is the dominant cause of soil saturation (McNamara et
al., 2005; Tague and Peng, 2013), whereas the partitioning of rain and
snowmelt may be more similar in wet climates and shallow soils, which are
more frequently saturated by either rain or snowmelt inputs (Loik et al.,
2004) (Fig. 1).

To evaluate soil moisture response to rainfall and snowmelt over a wide
range of climate and soil conditions, we used HYDRUS-1D (Šimůnek et
al., 1998), a physically based finite-element numerical model for simulating
one-dimensional water movement in variably saturated, multi-layer, porous
media.

2.1 Study design, site selection, and data sources

We utilized daily input data from five United States Department of
Agriculture Natural Resources Conservation Service (NRCS) snow telemetry (SNOTEL) sites in each of three regions that span a wide range of climate and snow conditions: the Cascades, Sierra Nevada, and Uinta mountains, for a total of 15 sites. Daily rather than hourly data were chosen to limit the effects of missing and incorrect values on the analyses. The three regions were chosen to represent dominant climate types in the western US, and within each region, sites were selected to span a snow persistence (SP)
gradient, where SP is the mean annual fraction of time that an area is snow
covered between 1 January and 3 July (Moore et al., 2015) (Fig. 2a, Table 1).

Figure 2(a) SNOTEL sites utilized for climate scenarios in this study with insets displaying snow zones classified by mean annual snow persistence (Moore et al., 2015). (b) Modeling domain layout with yellow points showing 5, 20, and 50 cm depths where volumetric water content time series were used for model calibration. The deepest yellow point is the depth where time series were extracted to calculate deep saturation. Symbols in the graph above the discretized soil profile represent the range of climate scenarios used plotted by mean annual precipitation (P) and mean annual temperature (T), and the three soil profiles below represent the soil parameter sets labeled with italicized capital letters: (a) loam,
(b) sandy clay loam, and (c) sandy loam. Different layers in each soil profile are represented as shades of gray; shading does not indicate any property of the soil layer.

With each climate scenario we simulated vertical profiles of volumetric
water content (VWC), which were depth-integrated to compute soil moisture
storage (S). For these simulations deep drainage (D) is any flux of water
downward below the root zone. Runoff (Q) is any water that does not
infiltrate into the soil, either because of saturated conditions or because
input rates exceed infiltration capacity. Evaporation (E) is direct
evaporation from the soil, and transpiration (T) is mediated by plant roots. For this study, these values are combined into evapotranspiration (ET) to represent the bulk loss of water to the atmosphere.

Daily precipitation (P), snow water equivalent (SWE), and volumetric water
content (VWC) at 5, 20, and 50 cm were obtained for each SNOTEL site using
the NRCS National Weather and Climate Center (NWCC, 2016) Report Generator
(Table 1). The records were quality controlled to ensure reasonable
precipitation, SWE, and VWC values as in Harpold and Molotch (2015).
Unrealistic values were removed (i.e., negative SWE and VWC below zero or above
unity); all daily VWC outside of 3 standard deviations from the mean
were removed, and a manual screening was performed on VWC data to identify
shifts and other artifacts not captured by the first two automated
procedures. Daily potential evapotranspiration (PET) was extracted from
daily gridMET (Abatzoglou, 2013) for the 4 km pixel containing each SNOTEL
site. This product uses the ASCE Penman–Monteith method to compute PET.

We chose three SNOTEL sites (432 Currant Creek, 698 Pole Creek R.S., 979 Van Wyck) to represent soil profile characteristics. While 365 of the 747 SNOTEL sites in the western US have soil moisture sensors, only a fraction of these sites have detailed soil profile data. The sites with soil profile data have information obtained from soil samples taken in the soil pits and processed in the NRCS Soil Survey Laboratory in Lincoln, NE, for texture, water retention properties, and hydraulic conductivity. We obtained detailed soil profile data in the form of pedon primary characterization files from the NRCS, and selected profiles (Fig. 2b, Table 2) that represent the range of soil textures and hydraulic conductivity values present at SNOTEL locations. Each had detailed NRCS pedon primary characterizations to depths greater than 100 cm and >15 years of daily soil moisture records
at 5, 20, and 50 cm depths.

2.2 Simulations

In HYDRUS-1D, we simulated water flow and root water uptake for a vertical
domain 10 m deep (Fig. 2b). The domain was discretized into 500 nodes with
higher node density near the surface (∼0.15 cm for the top 5 to ∼5 cm for the bottom of the profile). For the surface boundary, we used a time-variable atmospheric boundary condition, which allows specification of input (snowmelt and rain) and potential evapotranspiration (PET) time series. Runoff can also be generated at the surface boundary. For the lower boundary, we allowed free drainage from the bottom of the soil profile at 10 m. Surface soil water input was calculated by totaling snowmelt and rainfall input at the daily time step from SNOTEL precipitation and SWE values. Melt was computed for any day when SWE decreased; if SWE decreased, and the precipitation was greater than 0, total soil water input was assumed to be melt plus precipitation. The atmospheric boundary condition requires PET, leaf area index (LAI), and a radiation extinction coefficient used in the estimation and separation of potential evaporation and transpiration. We assigned a constant LAI of 3, as this value generally fits the mixed conifer forests (Jensen et al., 2011) where SNOTEL sites are installed. We assumed a radiative extinction coefficient of 0.39,
which is the default value. Root water uptake in the model was estimated
using Feddes parameters for a conifer forest (Lv, 2014: h1=0 cm, h2=0 cm, h3h=-5100 cm, h3l=-12800 cm, h4=-21500 cm, TPlow=0.5 cm d−1, TPhigh=0.1 cm d−1), with roots uniformly distributed from the soil surface to the interface with a lower hydraulic conductivity layer, as we lacked any more detailed information on root distribution or soil depth at these sites.

We created soil layers with depths and textures taken from the NRCS soil
pedon measurements. From this information we applied the neural network
capability of HYDRUS-1D, which draws from the USDA ROSETTA pedotransfer
function model (Schaap et al., 2001), to determine soil hydraulic
parameters. Using the NRCS pedon primary characterizations we input percent
sand, silt and clay, bulk density, wilting point, and field capacity. The
neural network model then estimates soil hydraulic parameters based on these
inputs. Below the depth of the soil pedon measurements, we configured the
simulations to have a deep “bedrock” or regolith layer with lower
saturated hydraulic conductivity (Ks) but the same water retention parameters as the layer above. Any water entering this lower layer is considered deep drainage. The hydraulic conductivity of this lower layer was set at one-tenth that of the layer above. This value was determined through iterative testing of Ks values (see Supplement). We extended the “bedrock” or regolith layer to 10 m depth to allow for deep drainage to occur without boundary effects that could be caused by a shallower regolith. The initial VWC for all layers in each simulation was 0.2, and simulations were run with a year of surface boundary condition inputs to establish initial conditions. We tested the simulation configuration by comparing to observed VWC at 5, 20, and 50 cm depths for the three selected soil profile sites (Fig. S1, Table S1 in the Supplement). Rather than force fitting, our goal was to produce simulations with similar timing of wetting and drying to observations. This approach is consistent with other studies using HYDRUS-1D, which also started with basic soils data and application of the ROSETTA pedotransfer function (Scott et al., 2000) and then calibrated to observed water content measurements by adjusting permeability of the “bedrock” layer (Flint et al., 2008).

We applied climate scenarios from each of the 15 SNOTEL sites selected
(Table 1) to each of the soil profiles to examine how climate and soil type
affect partitioning. We then conducted additional experiments to modify
inputs using just the loam profile. First, to examine whether snowmelt and
rainfall are partitioned differently, we changed all precipitation to rain
and compared the simulation output to those with the original climate data.
Second, to examine the effects of input concentration, the temporal
clustering of input through time, we artificially produced intermittent
input (four 5-day periods of low magnitude) and concentrated input (one
20-day period of high magnitude) of the same annual total for one wet (559) and one dry (375) site using the loam profile (1056) for all years of data. Third, to examine how soil depth affects partitioning, we altered the
depth of rooting zones to 0.5, 1.5, and 2 times their original depth. For 0.5 depth scenarios, we replaced soil layers deeper than 0.5 times the original depth with the bedrock/regolith layer. For 1.5× and 2× scenarios, the layer above bedrock/regolith was extended downward and the rooting zone extended to the new soil depth.

2.3 Analysis

Using the simulation results, we examined how rain and snowmelt were
partitioned into soil storage (S), deep drainage (D), evapotranspiration (ET), and runoff (Q). Daily soil storage is reported as the total soil water within the rooting zone only, and D is any water passing below the rooting zone (106–127 cm depending on the soil profile). We assessed partition components both in units of length (cm) and as ratios to total input (unitless, e.g., Q∕P) at both event and annual timescales.

To analyze hydrologic partitioning at event timescale we defined rainfall
events as days with precipitation while SWE equaled zero and snowmelt events
for days with declining SWE and no simultaneous precipitation. To focus on
differences between rainfall and snowmelt, only events with entirely
rainfall or entirely snowmelt input were considered in this analysis; mixed
events were excluded, though mixed input accounts for an average of 47 %
of annual input across all sites and years (Table S6). Events could last as
long as the conditions were continuously satisfied, and only those followed
by at least 5 days of no input were used in analysis. Total depths of
each variable were computed for each defined event time period. Input rain
and snowmelt were summed over the event time period, and response variables
(Q, ET, D) also included the day after the event ended to account for lag in event response. Antecedent S for each event was determined by taking the root zone storage from the day prior to the first event input.

At the annual scale, soil water input and partitioning components (rain,
snowmelt, Q, ET, D) were totaled for each year and the change in water year storage (ΔS) determined by subtracting the values of S at the end of the year from the value at the beginning of the year. In addition to ΔS, mean saturation (Sat) at each observed depth was calculated as the average annual VWC divided by soil porosity. We use mean saturation (Sat) as an alternative to change in water year storage (ΔS) because mean saturation is much easier to quantify at a field site than root zone storage, and this extends the application of our study to other areas with daily VWC data. Sat also provides a measure of soil water conditions throughout the year as opposed to ΔS, which represents only changes between the start and end of the water year.

To characterize climate conditions at the mean annual scale, each site was
classified as dry (precipitation deficit, PET >P) or wet
(precipitation surplus, PET <P). This separation by aridity index is
based on our hypothesis that the influence of concentrated snowmelt is
greater in dry climates than in wet climates (Hammond et al., 2018a). We also
report the maximum SWE and snowmelt fraction as the annual total snowmelt
divided by annual total input. Following the methods for computing the
precipitation concentration index (PCI), which represents the continuity or
discrete nature of input through time (Martin-Vide, 2004; Raziei et al.,
2008; Li et al., 2011), we computed the input concentration index (ICI)
using snowmelt and rain input. When calculated with daily data on an annual
basis, PCI commonly ranges between 0 and 100, where higher values correspond
to precipitation that is irregularly spaced in time and low values
correspond to precipitation evenly distributed throughout the year (Cortesi
et al., 2012). Our use of the terms input concentration and the input
concentration index refer to the temporal clustering of input in time and
do not refer to the intensity of melt. Pearson correlation tests were
conducted between explanatory variables (P, PET, P∕PET, peak SWE, average melt rate, and ICI) and dependent variables (Q, ET, D, mean saturation at 100 cm: Sat100).

Using both the event and annual results, we examined (1) whether
partitioning of rainfall input differed from that of snowmelt input, and
(2) how partitioning was affected by climate, soil texture, and soil depth. For question 1, we tested for differences in event partitioning between input
type (rain or snowmelt) and differences in annual partitioning between
historical and all rain scenarios using ANOVA. For question 2, we tested for
differences in annual partitioning between climate (wet, dry) and soil depth
groupings, also using ANOVA. Additionally, for question 2 we tested the
pairwise difference in linear regression slopes using the regression with
interaction test in JMP (SAS-based statistical software) to determine
whether the rate of change between explanatory and response variables
differed by climate or soil depth grouping. By comparing the slopes of
regressions run on standardized data, it is possible to assess the influence
of independent variables on dependent variables in different groupings. In
this study, we use this test to assess the influence of snowmelt fraction of
input and input concentration index on runoff and deep drainage response for
all, wet, and dry groupings as well as soil texture groupings.

Simulations for each of the 15 climate scenarios exhibit substantial
variability at the annual scale in precipitation (P), runoff (Q), and deep drainage (D) (Fig. 3). All regions have a wide range of annual P, but overall the highest P was in the Cascades region and lowest in the Uinta. The wide range of climate conditions simulated allows for an evaluation of climate effects on Q, ET, D, and Sat100 (Table S3). Annual precipitation (P) is positively correlated with runoff (Q, r=0.97), deep drainage (D, r=0.92), and Sat100 (r=0.73) (Table S3). The relationship is linear for Q but nonlinear for D and Sat100. Sat100 plateaus at ∼250 cm P with further P partitioned to Q instead of D. Evapotranspiration (ET) has the weakest correlations with P (r=0.08) of all partitioned components. Q∕P increases with P up to around 250 cm of P, and D∕P increases with P up to around 100 cm (Fig. 3). ET∕P decreases with precipitation, whereas ΔS∕P is unrelated to P. At values of P greater than around 300 cm, all variables have relatively consistent values even as P increases.

3.1 Snowmelt vs. rainfall and climatic influences on partitioning

Our first research question asks whether snowmelt and rainfall are
partitioned differently. At the event scale, input rates are significantly
greater on average for snowmelt than for rainfall in each of the three
regions and for the full dataset (ANOVA p<0.0001, mean snowmelt 1.1 cm d−1, mean rainfall 0.9 cm d−1, Fig. 4), though rainfall events have a higher maximum input rate (maximum snowmelt 8.0 cm d−1, maximum rainfall 14.7 cm d−1). Snowmelt events tend to occur on wetter soils, as estimated by antecedent soil moisture storage for the rooting zone (ANOVA p<0.0001, mean S for snowmelt 56.6 cm, mean S for rainfall 48.2 cm). Average runoff ratios (Q∕P) are higher for snowmelt than for rainfall (ANOVA p<0.0001, mean Q∕P snowmelt 0.20, mean Q∕P rainfall 0.03), whereas ET∕P is lower for snowmelt as compared to rainfall (mean snowmelt 0.24, mean rainfall 0.40). Deep drainage responses are affected by longer timescales than single events, so we did not include these in the event analysis. This event analysis only considered binary snowmelt or rainfall events.

At the annual scale, input at all sites is a mixture of rain and snowmelt.
To examine the importance of snow to partitioning, we used snowmelt
fraction, defined as the fraction of snowmelt to total precipitation, and
ICI. Snowmelt fraction and snow persistence are
generally positively correlated with ICI at dry sites in the Uinta and
Sierra, but this correlation declines with wetter sites in the Cascades
(Fig. S7). Q∕P increases with snowmelt fraction (r=0.41), most noticeably
where snowmelt fraction is >0.5, and increases with ICI (r=0.80) (Fig. 5). The ranges of Q∕P are higher in wet than in dry climates, though dry climates show greater rates of change with increasing snowmelt fraction and input concentration (Table S4). D∕P is somewhat correlated with snowmelt fraction (r=0.20) and ICI (r=0.43). D∕P ranges are higher in wet than in dry climates, and many years in dry climates do not generate D. ET∕P is not related to snowmelt fraction and generally declines with
ICI (r=-0.75). Ranges are lower for wet climates, where greater input is
partitioned to Q and D.

1 Values of ET/P>1 indicate root uptake from soil storage for years with low input (Fig. S5). 2 For a dry site (375) and a wet site (559). Intermittent simulations have annual total input separated into four 5-day periods, whereas concentrated input simulations have all input in a 20-day period of high magnitude.

Another effect of snow loss can be a decrease in input concentration.
Experimental scenarios with constant P separated into intermittent and
concentrated inputs for a wet site (375) and a dry site (559) show that
increasing input concentration leads to significantly greater Q∕P in the dry site (p<0.05, intermittent mean 0.54, concentrated mean 0.68, Table 3, Fig. 6) but no significant difference in the wet site. In contrast, D∕P is significantly greater (p<0.0001) for the concentrated input scenarios for both dry and wet sites, as no deep drainage is produced with intermittent input. ET∕P is significantly lower in concentrated input scenarios,
with a greater difference in dry climates (p=0.004, mean intermittent 0.80
vs. concentrated 0.66) than in wet climates (p=0.013, mean intermittent 0.34 vs. concentrated 0.28).

3.2 Soil property influences on partitioning

Soil stores water that may later be partitioned into Q, ET, and D. Using Sat100 as an indicator of soil moisture storage, Fig. 7 displays the
relationships between Q∕P, D∕P, and ET∕P vs. Sat100 as separated by climate type, soil texture, and root zone depth. Sat100 has strong relationships with Q∕P, D∕P, and ET∕P for all, wet, and dry sites (Fig. 7, Table S5). Q∕P is generally low (Fig. 7a, <0.3) until Sat100 is greater than >0.5. D∕P in the simulations also increases with Sat100, and many simulation years have limited D when Sat100<0.5. ET∕P generally decreases with saturation for Sat100 values >0.5.

Figure 7(a) Annual surface runoff (Q), deep drainage (D), and
evapotranspiration (ET) as a fraction of annual precipitation (P) vs. annual mean saturation at 100 cm depth (Sat100) and classified by climate type on the loam profile, dry sites P∕PET <= 1, wet P∕PET > 1. (b) The same variables displayed in (a) but classified by soil texture on three different soil profiles. (c) The same variables in A but classified by root zone depth on four different profiles of differing root zone depth. All simulations use historical input.

When these same relationships are separated by soil texture rather than
wet/dry climate (Fig. 7b, Table S5), the response patterns are similar
between soil types except for the sandy loam profile, which displays higher Q∕P and D∕P than the loam and sandy clay loam profiles at similar Sat100 levels. Differences between responses by soil texture are more evident at sub-annual timescales (Fig. 8a). For the example time period shown in Fig. 8a, the 100 cm depth in loam and sandy clay loam profiles wet up each spring during snowmelt 5 d prior to the sandy loam profile, and they generated deep drainage earlier and on more occasions than sandy loam due to higher water retention. The latter soils ultimately reached the highest annual D∕P values at higher Sat100 values, leading to more runoff generation via saturation excess, whereas the drier conditions in sandy loam led to more infiltration excess runoff. While this example time period and site display noticeable differences in cumulative response between soil textures, when the data for all sites and years are combined few significant differences in annual partitioning between soil textures emerge (Figs. 6 and 7).

Figure 8(a) Daily time series of cumulative runoff (Q), saturation at 100 cm depth (Sat100), and cumulative deep drainage (D) for SNOTEL site 698 input on SNOTEL site 515 (sandy loam), 1049 (sandy clay loam), and 1056 (loam) profiles. (b) Daily series for the same variables plotted for four depth scenarios 0.5×, 1×, 1.5×, and 2× original rooting zone depth.

To assess the influence of soil profile depths on partitioning, we altered
the loam soil profile to be 0.5×, 1.5×, and 2× times its original depth (Fig. 6, Table 3). For historical input, Q∕P and D∕P are greatest for the 0.5× depth scenario, and Q∕P declines significantly with deeper soils for both dry and wet sites (p<0.0001), with the greatest declines between 0.5× and 1× (original) depth. D∕P declines significantly between 0.5× and 1× depth, and then increases slightly for all sites with subsequent increases in depth to 1.5×
and 2× (Fig. 6, Table 3). Q∕P and D∕P differences by depth are significant between 0.5× and 1× depth, but not for all subsequent depth comparisons for all, wet, and dry site classifications (Table 3). In pairwise comparisons between depth scenarios Q∕P is only significantly different between 0.5× and 1× depth categories
(p<0.0001). Changes in ET∕P with soil depth are not significant according to ANOVA tests.

Figure 8b displays daily time series of surface runoff, deep saturation,
deep drainage, and cumulative deep drainage during an example period for the
four different soil root zone depth scenarios. The shallowest rooting zone
of 0.5× original depth produces the greatest surface runoff as well as
cumulative deep drainage throughout the example period. Each depth reaches
and remains at saturation for different amounts of time, with the shallowest
profile reaching saturation earliest and remaining saturated longest, but
also decreasing more rapidly to the lowest ending saturation. The deepest
profile takes the longest to increase Sat100, not reaching as high a peak, yet remaining higher at the end of the period. Deep drainage begins earliest for the shallowest depth scenario, though reaching a lower daily flux than the original depth. Deep drainage from the 1×, 1.5× and 2× original depth scenarios lags behind the 0.5× scenario following the same succession as their Sat100 patterns. These patterns in daily Sat100 and deep drainage
result in the highest cumulative deep drainage for the shallowest scenario.

The initial hypotheses for this study were that runoff and deep drainage
would be greater from snowmelt than rainfall and that snowmelt is more
important to generating runoff and deep drainage in deep soils and dry
climates than in shallow soils and wet climates. Our results indicate that
snowmelt is an efficient runoff generator, though not necessarily an
efficient generator of deep drainage. Deep drainage is less affected by
input type because it is controlled by deep soil moisture patterns over
longer timescales. Soil texture modifies daily wetting and drying patterns
but has limited overall effects on annual partitioning, whereas increases in
soil depth decrease runoff and increase deep drainage. Overall these results
indicate that runoff may be substantially reduced with seasonal snowpack
decline in all climates, whereas the effects of snowpack decline on deep
drainage are less consistent. We expand on these key findings in the
paragraphs below and suggest that areas in dry watersheds with storage
similar to peak SWE may be most likely to experience reductions in deep
drainage with continued slow loss.

Multiple lines of evidence confirm snowmelt as a more efficient runoff
generator on average than rainfall. At event scale runoff efficiency was
elevated for snowmelt because of the 22 % greater input rate and 17 %
wetter soils than rainfall. This is consistent with previous studies showing
that snowpack development and subsequent melt tend to occur when soils are
at elevated moisture contents due to lower ET (Liu et al., 2008; Williams et
al., 2009; Bales et al., 2011). The effects of snowmelt vs. rainfall are
weaker at annual timescales (Fig. 5, Table S3) because these longer time
periods include a combination of snow, mixed, and rainfall inputs in
contrast to the event analysis in which we analyzed only events that were
exclusively snowmelt- or rainfall-dominated. Forcing all input into the
extreme case of all rain produces 67 % declines in runoff efficiency (dry:
0.13 vs. 0.04; wet: 0.46 vs. 0.29) (Table 3, Fig. 6), likely because the
input becomes less concentrated in time for the all-rain scenario, allowing
more ET. We also hypothesized that the effects of changing snowpacks would be
greatest in dry climates, where soil saturation is less frequent. However,
simulations suggest that both wet and dry climates are as likely to show
reduced surface runoff with declining snow water inputs.

The effects of snow loss on D were not as consistent across our simulations
as the effects on Q. Prior research has demonstrated strong seasonality in
groundwater recharge, attributable to thresholds in input intensity
(Jasechko and Taylor, 2015) and seasonal differences in evapotranspiration
(Jasechko et al., 2014, 2017). We had hypothesized based on additional research (Hunsaker et al., 2012; Langston et al., 2015; Barnhart et al., 2016; Li et al., 2017; Hammond et al., 2018a) that input concentration along with evapotranspiration seasonality would be the primary reason for
elevated Q and D from snowmelt relative to rainfall. In this study, changes from snow to rain both increased and decreased D∕P (Figs. 6 and S2), and D∕P was not correlated with either snowmelt fraction or ICI in
wet climates. In general, Q∕P was greater than D∕P, and the D∕P response to changing input was weaker because S mediates the connection between input and D. In the 1-D model Q is affected by infiltration rate and near-surface storage and can more rapidly respond to input changes. In the simulations shown here once subsurface storage is zero, D will plateau, and Q will increase with further input due to the saturation excess mechanism. Although these processes were simulated in 1-D, they are consistent with observations of saturation excess overland flow documented in the elevation bands of many SNOTEL sites (Newman et al., 2004; Eiriksson et al., 2013; Kampf et al., 2015). In wet climates, D∕P is less affected by input type because conditions are more likely to be wet, regardless of whether input is snow or rain. D∕P is more affected by
changes from snow to rain in dry climates, likely because of the role that
concentrated snowmelt can play in allowing water to reach the base of the
soil column.

Soil texture and depth generally did not change partitioning at the annual
timescale as much as the varying climate scenarios (Fig. 6), with the
exception of changes in the shallowest soils (1× depth to 0.5× depth results in 12 % Q increase and 180 % D increase). D∕P generally increased with increasing soil depth, demonstrating the importance of lower boundary conditions to shallow vs. deep partitioning. Altering soil profile depth and the associated root zone depth produced the largest effects on Q∕P and D∕P from 0.5× to 1× depth. The responsiveness of fluxes to changes in soil depth from 0.5 to 1× may relate to storage capacity relative to input. The soil depths ranged from 106 to 127 cm, which with a porosity of 0.4 gives a storage capacity of 42–51 cm, large enough to store the mean annual precipitation in some dry watersheds. When this storage is reduced by half to 21–25 cm, it is smaller than the mean annual precipitation at the wetter sites, increasing the likelihood of soil saturation that leads to D and Q. Consequently, the change in profile depth from 0.5 to 1 m represents a shift from annual input greatly exceeding profile storage to storage approximately accommodating annual input. At the sites used in this study, mean annual P ranged from 0.8 to 11.3 times the storage of the 1× soil profile, and peak SWE ranges from 0.1 to 5.9 times the storage. Prior field-based studies have also documented SWE that is similar in magnitude to the maximum amount of water storage in the upper meter of soil (Bales et al., 2011) and have shown that reducing soil depth increases surface runoff and deep drainage (Smith et al., 2011).

Focusing on the influence of soil texture, simulations indicate that shorter
durations of deep drainage for the coarser sandy loam compared to the finer
texture soils are offset by higher rates of flux during deep drainage in the
coarser profile (Fig. 8a). Similarly, lower likelihood of surface saturation in the sandy loam soil compared to other soils is offset by greater likelihood of infiltration excess runoff. Therefore, in this 1-D approach, soil depth exerts a stronger control on annual total input partitioning to Q and D, whereas soil texture has a limited effect on annual partitioning but can affect the timing of partitioning and water availability during different times of year. In natural landscapes, texture differences can result in spatially variable soil moisture (Williams et al., 2009; Kaiser and McGlynn, 2018). Combined variations in soil texture and depth within a watershed may result in significant differences in soil moisture storage across the basin (Bales et al., 2011), resulting in substantial differences in response throughout a watershed. The distribution of soil water storage capacity across the watershed likely exerts a strong control on locations where surface runoff, streamflow generation, and deep drainage are most efficiently generated, especially in dry watersheds where soil moisture is generally low except during snowmelt (Atkinson et al., 2002; Seyfried et al., 2009). Additionally, unsaturated soil water storage may be the dominant control on streamflow activation during dry periods, while total input depth is the dominant control on streamflow generation during wetter periods (Farrick and Branfireun, 2014). Combining the role of soil storage capacity in space and time, areas in dry watersheds with storage similar to peak SWE may be most likely to experience reductions in deep drainage with continued slow loss.

4.2 Limiting assumptions

Given the complex nature of soil water movement in heterogeneous mountain
topography, this study makes several assumptions and simplifications. The
simulations do not include the intricacies of vegetation water use, assuming
a static leaf area index (LAI) with root uptake controlled only by PET and
soil moisture, and we assume free drainage from the bottom boundary of the
modeled domain. Changing static LAI has a substantial effect on soil
moisture dynamics (Chen et al., 2014), though model performance to match
simulated and observed soil moisture does not necessarily improve with the
assimilation of dynamic LAI values (Pauwels et al., 2007). Incorporating
site-specific constant LAI from field measurements or remotely sensed data
may have improved model performance, especially during spring green-up and
fall senescence, and is recommended for future site-specific studies. The
water balance in hydrologic models can be highly sensitive to the method
chosen to represent root uptake and plant water use (Gerten et al., 2004),
and hydrologic models generally poorly capture or replicate the interactions
between soil, vegetation, and atmospheric properties that combine to control
plant water use (Gómez-Plaza et al., 2001; Gerten et al., 2004; Zeng et
al., 2005). In addition, we did not allow for frozen soils in our
simulations, but this can be a strong influence on soil input partitioning
in places where snow depth was <50 cm and incapable of insulating the soil (Slater et al., 2017).

Additionally, simulations are generally wetter than measured water contents;
therefore, the representation of partitioning shown here displays relative
response between climates and soil profiles rather than absolute
quantification of these partitioned components. The profile depths we
simulated represent the minimum likely soil depth, as the collection of the
pedon reports was limited by the depth of refusal for sample collection.
Shallow soil profiles can also lead to a wet bias in simulations, and this
can artificially elevate saturation excess flow, leading to our observations
of greater Q∕P than D∕P in most site years. Our modeled domain included an extended “bedrock” or regolith layer to 10 m depth to allow for deep
drainage without lower boundary effects. The choice of lower boundary
condition affects the simulation of soil moisture and water balance
partitioning, with free drainage generally resulting in lower soil moisture,
evapotranspiration, and runoff than with a no-flux boundary condition
controlled by an impervious layer or fluctuating water table (Chen et al.,
2018). We created a domain much deeper than the soil zone to minimize this
boundary condition effect; effects of lower boundary conditions are
generally seen in deeper layers of the soil profile and during transition
periods between soil water input events when capillary rise can influence
transpiration and deep drainage (Leterme et al., 2012; Brantley et al.,
2017). Though a no-flux boundary condition may be appropriate for sites
where relatively shallow water tables exert a strong influence on soil
moisture dynamics, the inclusion of a no-flux lower boundary for the sites
in this study would have made simulations wetter, furthering the difference
between observed and modeled VWC.

Sub-daily dynamics in snowmelt and ET are not captured by our use of a
daily time step. We chose to model soil water response to rainfall and
snowmelt at the daily time step due to better data quality, but processes
affecting partitioning of these inputs take place at sub-daily scales.
Comparisons of results from simulations using daily vs. hourly input
demonstrate similar timing of response but greater cumulative surface
runoff from hourly simulations and greater cumulative deep drainage from
daily simulations (Table S2). The short hourly time period allows for higher
intensity input, which causes infiltration excess overland flow, whereas
daily input is of lower intensity, allowing for greater deep percolation.
Additionally, SNOTEL sites do not have measured values of PET, so we relied
on a modeled 4 km gridded product (Abatzoglou, 2013), which may better
represent some sites than others. It was beyond the scope of this study to
perform a sensitivity analysis of PET data source.

Hydrologic response in hillslopes and catchments is not fully captured in
the 1-D modeling approach. Water partitioned into Q and D in a 1-D model may not represent the same Q and D observed at a stream: Q generated at a point location may reinfiltrate downslope; D may also emerge downslope to supply streamflow rather than remaining in the deep subsurface. Topography affects both soil moisture and snow patterns (Western et al., 2004; Liator et al., 2008; Williams et al., 2009; Brooks et al., 2015), and it leads to lateral surface or subsurface flow, which can be important in redistributing water
downslope along the soil–snow interface (Webb et al., 2018) and within the
shallow subsurface (Kampf et al., 2015, Kim et al., 2016). Lateral
redistribution of water thus leads to spatially variable patterns of input,
storage, runoff generation, and ET at the hillslope to watershed scales (Brooks et al., 2015). While simulating only vertical flow is reasonable for
SNOTEL sites located in relatively flat forest openings, 1-D simulations
will tend to be biased wet because they do not allow any lateral redistribution. A progression of the work shown here would be to simulate 3-D flow (ex. Weiler and McDonnell, 2007; Seyfried et al., 2009) and examine the spatial variability in effects of snow loss. For example, a decline in deep drainage near a ridge line, where flow paths are predominantly vertical, could reduce subsurface flow emergence at downslope locations, and this decreased groundwater emergence may reduce ET in areas where vegetation is reliant on the emergence of deeper flow paths.

The simulations used here only allow for matrix flow, excluding macropore
flow, for a simplified representation of soil water movement. Preferential
flow though the profile can enhance deep drainage relative to surface
runoff, which is another potential reason why soil moisture simulations were
biased wet; 60 %–80 % of deep drainage has been shown to occur as
preferential rather than interstitial flow (Wood et al., 1997; Jaynes et
al., 2001; Sukhija et al., 2003), although the amount of preferential flow
varies substantially between climates and soils. The magnitudes of fluxes in
our simulations are consistent with observation studies, however, lending
more confidence to the simplified modeling approach. Simulated annual D∕P for dry climates (∼0.05) is similar to values reported from observations (Wood et al., 1997). The simulated Q∕P (∼0.0–0.9) vs. snowmelt fraction plots from HYDRUS-1D simulations follow the same general increasing pattern (r=0.41) as observed Q∕P (∼0.0–1.0) vs. SP in Hammond et al. (2018a) (r=0.39).

Future work could examine the potential sensitivity of the results to these
limiting assumptions, In particular, researchers could examine the extent to
which adding spatially and temporally varying vegetation processes and/or
preferential flow pathways would change water balance partitioning. Simulations could expand to two dimensions to examine how downslope affects
partitioning from ridgelines to valley bottoms or to three dimensions to
examine effects of flow convergence and exfiltration in hillslope hollows.
Because of the complexity of subsurface properties, this work would also
benefit from more information about the hydraulic properties of the deep
subsurface below the measured soil pedons. This could be linked with model
analyses examining how both subsurface properties and boundary conditions
affect the simulations.

This study helps to explain the mechanisms that lead to greater runoff from
snowmelt. At event scale snowmelt generates more runoff because it tends to
have a greater input rate and occurs on wetter soils than rainfall. Seasonal
snowmelt elevates runoff in both wet and dry climates. Deep drainage can
also decline with loss of snow, but it has a weaker response because soil
storage buffers the impacts of snow loss. Soil properties can mediate the
effects of snowmelt on rainfall changes, with soil depth having a greater
effect than texture on input partitioning, particularly where soil water
storage is less than mean annual precipitation. Soils that are shallower
than observed soil depths generate the greatest runoff and deep drainage,
indicating that shallow soils may show the largest changes in partitioning
as input transitions from snowmelt to rainfall. Increasing soil depth above
observed depths gradually reduces surface runoff while increasing deep
drainage. Soil texture modifies short-term timing of soil moisture and
runoff generation, but these effects are not large enough to alter the
annual response of different soil types to changes in snow. The 1-D
simulations provide basic hypotheses for hydrologic partitioning under
changing snowmelt that should be further explored with 2-D or 3-D
hydrological models and direct observations. Although more work is necessary
to translate these findings to watershed-scale streamflow response, the
findings highlight the importance of precipitation-phase shifts for runoff
generation and groundwater recharge.

JH, AH, and SK designed the experiments and JH and SW carried them out. JH and SW performed the simulations. JH conducted statistical analyses on model outputs. JH prepared the manuscript with contributions from all the co-authors.

Li, D., Wrzesien, M. L., Durand, M., Adam, J., and Lettenmaier, D. P.: How
much runoff originates as snow in the western United States, and how will
that change in the future?, Geophys. Res. Lett., 12, 6163–6172, https://doi.org/10.1002/2017GL073551, 2017.

Tague, C. and Peng, H.: The sensitivity of forest water use to the timing of
precipitation and snowmelt recharge in the California Sierra: Implications
for a warming climate, J. Geophys. Res.-Biogeo., 118, 875–887, https://doi.org/10.1002/jgrg.20073, 2013.

Streamflow in high-elevation and high-latitude areas may be vulnerable to snow loss, making it important to quantify how snowmelt and rainfall are divided between soil storage, drainage below plant roots, evapotranspiration and runoff. We examine this separation in different climates and soils using a physically based model. Results show runoff may be reduced with snowpack decline in all climates. The mechanisms responsible help explain recent observations of streamflow sensitivity to snow loss.

Streamflow in high-elevation and high-latitude areas may be vulnerable to snow loss, making it...