Depth‐ and Time‐Resolved Distributions of Snowmelt‐Driven Hillslope Subsurface Flow and Transport and Their Contributions to Surface Waters

Major components of hydrologic and elemental cycles reside underground, where their complex dynamics and linkages to surface waters are obscure. We delineated seasonal subsurface flow and transport dynamics along a hillslope in the Rocky Mountains (USA), where precipitation occurs primarily as winter snow and drainage discharges into the East River, a tributary of the Gunnison River. Hydraulic and geochemical measurements down to 10 m below ground surface supported application of transmissivity feedback of snowmelt to describe subsurface flow and transport through three zones: soil, weathering shale, and saturated fractured shale. Groundwater flow is predicted to depths of at least 176 m, although a shallower limit exists if hillslope‐scale hydraulic conductivities are higher than our local measurements. Snowmelt during the high snowpack water year 2017 sustained flow along the weathering zone and downslope within the soil, while negligible downslope flow occurred along the soil during the low snowpack water year 2018. We introduce subsurface concentration‐discharge (C‐Q) relations for explaining hillslope contributions to C‐Q observed in rivers and demonstrate their calculations based on transmissivity fluxes and measured pore water specific conductance and dissolved organic carbon. The specific conductance data show that major ions in the hillslope pore waters, primarily from the weathering and fractured shale, are about six times more concentrated than in the river, indicating hillslope solute loads are disproportionately high, while flow from this site and similar regions are relatively smaller. This methodology is applicable in different representative environments within snow‐dominated watersheds for linking their subsurface exports to surface waters.


Introduction
Subsurface flow of water in mountainous watersheds influences rates of river discharge and water quality (Ameli et al., 2017;Carroll et al., 2018;Welch & Allen, 2012;Winnick et al., 2017), elemental cycling (Boy et al., 2008), and rock weathering (Anderson & Dietrich, 2001;Bluth & Kump, 1994;Brantley, Lebedeva, et al., 2017). The importance of hillslope and catchment hydrologic and geochemical responses to recharge events has motivated investigations in a wide variety of settings, some which have supported decades of research. Numerous studies across different catchments within the H. J. Andrews Experimental Forest in the western Cascade Range of Oregon have elucidated influences of land management and climate on hydrologic processes and water quality (Harr, 1977;Jones, 2000;McGuire et al., 2007). Research on the steep, humid, forested catchments at Maimai, New Zealand began in the late 1970s and has revealed the importance of preferential flow of rainfall recharge along soil horizons and interfaces with bedrock (McDonnell, 1990;McGlynn et al., 2002;Mosley, 1979;Sklash et al., 1986), bedrock topography , as well as significant lateral variability in flow (Woods & Rowe, 1996). At Coos Bay, Oregon, saturated and unsaturated flow and weathering have been quantified in a steep headwaters catchment Anderson & Dietrich, 2001;Montgomery et al., 2002;Torres et al., 1998). This introduction is not intended to be a comprehensive review of hillslope hydrology but points to research around the world that underscores the importance of better understanding subsurface flow (Brantley, McDowell, et al., 2017;Gu et al., 2018;Hughes et al., 2015;Sidle, 2015;Tetzlaff et al., 2015).
Quantification of subsurface and surface components of flow and transport requires extensive measurements of subsurface properties, pore waters, and overland flow (when present), yet collection of sufficient information needed to mechanistically close water and chemical cycles remains a challenge. While hydraulic conductivity (K) measurements in shallow soils are routine and often numerous within individual study sites, deeper measurements are commonly limited to tests on sparsely space piezometers and wells. Although continuous recording of water table elevations is routine (Leung et al., 2011;McDaniel et al., 2001;Montgomery et al., 2002;Uchida et al., 2004;van Meerveld et al., 2015), direct and continuous measurements of matric potentials needed to determine directions of vadose zone fluxes are less common, especially beyond shallow soil depths (Anderson & Burt, 1978;O'Geen et al., 2018). Troughs and trenches provide continuous measurements of runoff, shallow groundwater flow, and water chemistry, yet suitable collection lengths need to be determined, and limitations on installation depths preclude capture of deep flow paths (McGuire et al., 2007;Woods & Rowe, 1996).
Economic and logistical considerations allow only sparse pore water sampling at the hillslope scale. Soil and vadose zone pore water sampling with soil water samplers or suction lysimeters Boyer et al., 1997;Sklash et al., 1986;Sprenger et al., 2015) is common within shallow depths, but is not only discontinuous, but ineffective during periods when matric potentials become only moderately low. Moreover, pore waters obtained from suction samplers are collected over times that depend strongly on the matric potential (supporting information). Finally, pore water and solute sampling at depths greater than a few tens of meters is uncommon. Thus, sparse information on the permeability distribution, the hydraulic potential field (especially in the vadose zone), and the spatiotemporal distribution of pore water solutes limits the resolution with which subsurface processes can be linked to their integrated outputs measured in surface waters.
In contrast to the limited information on temporal and spatial distributions of subsurface solute fluxes leaving hillslopes and catchments, there exists a wealth of data on C-Q relations for solutes in streams and rivers (Godsey et al., 2009;Herndon et al., 2015;Moatar et al., 2017), and these are providing sensitive indicators of climate change and other anthropogenic influences on the hydrologic cycle (Dawson et al., 2008;Godsey et al., 2009;Todd et al., 2012). Still, decades after the introduction of the chemical and isotope mass balance approaches (Dinçer et al., 1970;Pinder & Jones, 1969), many approaches to hydrograph analyses rely largely or solely on measurements in streams and rivers to infer subsurface sources (Chanat et al., 2002;Miller et al., 2015;Stewart et al., 2007). Thus, there remains a need to develop practical strategies that provide mechanistic predictions of how different subsurface compartments export water and solute.
While the importance of understanding distributions of fluxes in hillslope environments has long been appreciated, large uncertainties are still commonly associated with flow and transport because of limited and sparse measurements. Phenomena such as groundwater ridging (Gillham, 1984;Sklash & Farvolden, 1979;Waswa & Lorentz, 2015), perched water flow (McDaniel et al., 2001;Salve et al., 2012;van Meerveld et al., 2015;Wilson et al., 1990), bypass flow by piping (Bryan & Jones, 1997;Jones, 1971;McDonnell, 1990;Uchida et al., 2005;Wilson et al., 2016), and transmissivity feedback Kendall et al., 1999;Laine-Kaulio et al., 2014;Laudon et al., 2004;Rodhe, 1989) have been invoked to conceptually explain field observations. Given that mechanistic predictions of subsurface flow and transport are based on Darcy's law, information on distributions of K and hydraulic potentials are essential. Although impacts of different hypothetical K distributions on hydraulic responses to recharge can be efficiently examined numerically (Hazenberg et al., 2016;Sprenger et al., 2016;Weiler & McDonnell, 2004), it is expensive to characterize distributions of permeabilities and hydraulic potentials in hillslopes over the lateral extents and depth ranges needed to assure that the flow field is adequately known. Thus, it is desirable to develop approaches for examining subsurface hydraulic responses that rely on shallower measurements yet lead to constraints on deeper flow.
Measurements to be described later revealed lack of a perched water table and K profiles that decreased with depth. Because these are key prerequisites of the transmissivity feedback model, our analysis builds on this 10.1029/2019WR025093 Water Resources Research approach which will be briefly described here. In many hillslope environments, the strata within which water table fluctuations occur have K that increase toward the soil surface. In such settings, downslope subsurface flow is amplified when the water table rises, not only because of the increased thickness of the saturated zone, but also because of the enhancement of transmissivity by inclusion of shallower zones with higher K. Transmissivity feedback describes the process where recharge from rainfall or snowmelt not only raises the water table but also drives a significant fraction of flow rapidly along the previously unsaturated zone overlying the pre-event water table, thereby mobilizing pre-event vadose zone water Rodhe, 1989).
Even when the spatial distributions of K are reasonably well constrained, an impermeable lower boundary needs to be assigned at some depth in order to define the subsurface flow domain. The bottom boundary is often defined to have a constant slope that is either horizontal (Ameli et al., 2018), parallel to the hillslope surface (Hazenberg et al., 2016;Salvucci & Entekhabi, 1995), or with a dip that differs from that of the topography (Troch et al., 2003). A commonly used alternative to assigning a depth for the impermeable lower boundary is to assume that K decreases exponentially with depth (Ameli et al., 2016;Beven, 1997;Brooks et al., 2004;Jiang et al., 2009). Although K profiles measured in some soils and shallow weathering rock do show rough exponential decay with depth (Brooks et al., 2004;Elsenbeer et al., 1992;Harr, 1977), other studies have shown more complex, nonmonotonic depth dependence (Meyer et al., 2014;Ross & McElwee, 2007). The lower boundary of the Maimai watershed has been previously described as effectively or nearly impermeable (Mosley, 1979; to poorly permeable, with about 100 mm of annual seepage losses to deep groundwater (McGlynn et al., 2002). However, a recent study that included K measurements in 1.51 to 8.80 m deep wells drilled into the bedrock yielded a range of K from 6.3 × 10 −9 m s −1 to 2.7 × 10 −7 m s −1 , with a mean K of 5.6 × 10 −8 m s −1 . Indeed, numerous studies have documented moderately high K, >10 −7 m s −1 , in rocks at depths of tens to hundreds of meters (Chen et al., 2015;Cook et al., 1996;Meyer et al., 2008Meyer et al., , 2014Welch & Allen, 2014), and even to several km (Kuang & Jiao, 2014;Ranjram et al., 2015).
Based on spatially and temporally resolved measurements of hydraulic and geochemical characteristics, we present a simple approach for characterizing subsurface distributions of flow and nonreactive solute transport. Subsurface flow and transport were analyzed within the lower section of a hillslope transect based on hydraulic and geochemical measurements within soil, weathering bedrock, and fractured bedrock to depths of about 10 m. By balancing time-dependent transmissivity fluxes against estimated annual net recharge (P minus evapotranspiration ET) entering the transect and its upslope contributing area, we calculated depthand time-dependent downslope fluxes through the soil, the underlying weathering rock zone, and deeper fractured parent rock. The advantages of this mass balanced transmissivity feedback analysis include (1) obtaining calculate subsurface discharge responses that are consistent with measurements of K and hydraulic potentials, (2) estimation of the depth to an operationally defined impermeable lower boundary, and (3) calculation of subsurface C-Q relations. These subsurface C-Q relations predict the time flow and Darcy flow rate dependence of pore water solute concentrations exported from the hillslope, based on subsurface hydrologic and geochemical measurements. This approach also provides guidance on additional measurement needs, especially with respect to deeper bedrock, and can be improved with additional data.
This study was conducted primarily within water year (WY) 2017 and WY2018, which contained substantially above-and below-average snowfall, respectively. The lower snowpack during WY2018 was followed by earlier snowmelt, characteristic of responses to warming trends observed over the past several decades (Clow, 2010). During both WYs, the water table along our monitored transect oscillated continuously between minima and maxima prior to snowmelt and immediately following snowmelt. As shown later, predictions of solute fluxes toward the river can be obtained based on combining time-dependent, downslope transmissivity fluxes through the soil, weathering zone, and fractured bedrock, with the fairly distinct pore water geochemistry in each of these layers. The sequence of high and low snowpack years provided an opportunity to observe contrasting responses that can be used to help improve estimates of K and fluxes and can help explore subsurface responses to warming trends.

Site Description
The field site is situated in the Rocky Mountains, within the upper East River watershed, 5.3 km southsouthwest of Gothic, Colorado, on a northeast facing hillslope that drains into a floodplain along the East 10.1029/2019WR025093

Water Resources Research
River (Dwivedi et al., 2018;Hubbard et al., 2018). The East River drains into the Gunnison River, a major tributary of the Colorado River. Measurements obtained for this study come primarily from three instrumented stations (PLM1, PLM2, and PLM3) spanning an elevation change of 27.1 m over a horizontal distance of 137 m along the lower portion of the hillslope (Figure 1). Beyond the lowest hillslope station (PLM3), at 2.1 m lower elevation and a horizontal distance of 43 m, an additional monitoring station (PLM4) was installed on the floodplain, about 50 m away from the East River. An additional location 35 m upslope from PLM1, designated PLM0, was only characterized for basic physical and chemical properties. As shown in Figure 1, the hillslope extends beyond the instrumented transect 690 m in the southwestward direction to a peak (elevation 2,936 m). The shallow soil profiles are vegetated with grasses, forbs, and shrubs representative of the regional lower montane life zone and underlain by Cretaceous Mancos Shale. Mancos Shale is broadly distributed throughout the southwestern

Water Resources Research
United States and releases elevated concentrations of major ions, trace elements, and potential contaminants during weathering (Morrison et al., 2012;Tuttle et al., 2014aTuttle et al., , 2014b. The area has a mean annual temperature of 3°C, mean annual precipitation of about 680 mm (~70% as snow). The Community Land Model (CLM), (Oleson et al., 2013) was recently used to delineate spatiotemporal variations in ET with 900 m resolution across the East River Watershed (Tran et al., 2019). Application of this analysis to the hillslope field site for WY2017* and WY2018 yielded ET of 310 and 288 mm, respectively. WY2017* refers to the interval from 1 December 2016 to 30 November 2017, rather than the standard WY, because continuous logging of water table levels was initiated on 1 December 2016. Seasonal characteristics of temperature, P, and ET ( Figure 2) support snowpack accumulation from November through mid-April, followed by warming-driven snowmelt and a rapid rise in ET. Note that the CLM-based calculations show that the onset of rising spring ET shifted earlier by about 10 days in WY2018 relative to WY2017 (Figure 2c). The North American summer monsoon brings significant additional P, although it occurs during months with high ET. Over the interval from snowmelt through mid-September in WY2016, 2017, and 2018, the calculated ET minus P amounted to 68, 180, and 128 mm, respectively, indicating that ET was large enough to prevent net recharge to groundwater from monsoon precipitation. This inference will be supported later with soil moisture profiles collected in late summer of 2016 and 2018.

Drilling, Instrumentation, and Monitoring
A track mounted drill rig was used to drill five boreholes in September 2016. Sites PLM1-4 were drilled down to a depth of 10.0 m below ground surface (bgs) with a 0.14 m diameter ODEX bit, with ejected drill cuttings collected in 0.6 to 0.75 m depth increments. An additional borehole (PLM6) was first drilled to obtain core samples, then for installation of a 0.051 m ID monitoring well slotted over the depth interval of 6.10 to 9.15 m bgs. The hillslope boreholes (PLM1,2,3) were instrumented to monitor vadose zone and groundwater interactions in a manner similar to that of a previous study (Tokunaga et al., 2016;Wan et al., 2018). Most instrument packages included a pore water sampler (pressure/vacuum soil water sampler, Soilmoisture Equipment 1920F1L06-B02M2), matric potential sensor (Decagon MPS-6), water content sensor (Decagon 5TE), and gas sampler (1.3 cm diameter, 2.6 cm long porous stainless steel filter tips, 50 μm pores, McMaster-Carr). Tubing (1/4 inch OD, stainless steel for gas samplers and nylon for water samplers) and electrical leads (MPS-6 and 5TE) were bundled together in lengths sufficient to extend above ground. For installation of instruments, a pump was lowered to the bottom of the hole to drain inflowing water and to allow video inspection of the borehole wall for identifying fractured zones to target for instrumentation. The camera and pump were removed, and specific depth intervals shown in Figure 3 were equipped with instrument clusters, backfilled with 150 to 425 μm silica sand (typically 15 to 30 cm below and above instruments), and isolated with bentonite (3/8" chips). The uppermost segments of boreholes were backfilled with native soil. The shallowest instruments were emplaced in individual 0.09 m diameter holes, hand-augered to the desired depth. It should be noted that flow of shallow sediments under positive pore pressures of the near-river floodplain environment caused the original PLM4 borehole to collapse before instruments could be installed. Therefore, the PLM4 instruments were only installed to 1.28 m bgs, in individual handaugered boreholes.

Soil Sampling and Characterization
Prior to drilling, soils were collected in 0.10 m depth increments by hand augering and coring, with samples sealed in Ziploc freezer bags and shipped back to the laboratory. These samples were used to determine

10.1029/2019WR025093
Water Resources Research depth profiles of particle-size distributions, water contents, and water potential profiles. The water potentials were determined with a Decagon WP4C, on subsamples taken from interior regions of intact cores and aggregates. An additional set of soil samples were collected in September 2018 at the hillslope stations for additional water potential measurements.

Hydraulic Conductivities
Hydraulic conductivities of soils along the hillslope were measured with a Guelph permeameter (Reynolds & Elrick, 1987;Zhang et al., 1998). The Guelph permeameter was also used to obtain measurements of K at a depth of 1.40 m within weathering shale in the vicinity of PLM3 and PLM6. Hydraulic conductivities in the fractured shale where measured at PLM1 (6.3 to 7.3 m depth bgs screened) and PLM6 (4.8 to 9.8 m depth bgs screened), based on bailing and monitoring water level recovery rates (Bouwer & Rice, 1976). Assignment of suitable effective K values is recognized to be both essential and challenging for predicting subsurface flow (Binley et al., 1989;Brooks et al., 2004). Therefore, impacts of varying K values relative to initial best estimates from field measurements are also explored.

Hydraulic Potentials
Water table depths at PLM1 and PLM6 were monitored with continuously logging pressure transducers (AquaTROLL 200) that were periodically checked with manual tape measurements of depths to water. Depths to the water table at PLM1, 2, and 3 were also obtained from equilibrium pressure measurements in pore water samplers using the "tensisampler" approach (Tokunaga, 1992). Although equilibrated tensisampler readings have high accuracy (±0.02 m H 2 O) for low magnitude matric head, response time limitations make them unreliable for matric head values more negative than about −1 m H 2 O. On the other hand, MPS-6 measurements are not able to resolve matric and pressure head greater than −1.0 m H 2 O. Therefore, by combining MPS-6 and tensisampler measurements, a very wide range of matric/pressure potentials can be monitored without loss of resolution at near-zero matric/pressure potentials. Procedures involved in estimating the depth of the water table based on trends in moisture and matric potential sensors (Decagon 5TE and MPS-6, respectively) are provided in the supporting information. At individual stations, vertical hydraulic head profiles take their local soil surface as the gravitation potential datum. Coordinates and elevations at the soil surface of each PLM station, were determined by GPS, thereby allowing determination of hydraulic head gradients across the hillslope transect.

Transmissivity Feedback Calculations of Subsurface Flow
As shown later, higher K is associated with shallower depths, and the water table rises and declines within the instrumented section of the hillslope without evidence of surface runoff and perched water tables. These observations justified application of transmissivity calculations, where transmissivity is the product of K times the saturated thickness of a given zone, to characterize subsurface flow from the hillslope transect toward the river. This approach calculates fluxes through zones below the dynamic water table based on the product of water table height-dependent transmissivities and the hydraulic head gradient. The latter was equated with the difference in water table elevation at PLM1 relative to PLM6, divided by their horizontal distance of 135.3 m. Based on both the K values and pore water chemistries, the subsurface was simplified into a 0.5 m thick surface soil zone, a 0.5 m thick subsoil zone, 2.6 m thick weathering shale zone, and an underlying fractured shale of thickness to be estimated through a water mass balance calculation. Within each zone, the average of measured K is multiplied by its water table elevation-dependent saturated thickness to obtain its time-dependent transmissivity. The influence of capillarity in enhancing the thickness of the effectively saturated soil was approximated by applying hysteretic tension saturation matric head values of −0.05 (estimated) and −0.12 m (measured) for wetting and draining, respectively, based on soil moisture retention relations ( Figure S2). Capillary effects were not added to the weathering zone because of very low magnitude matric potentials expected for draining major fractures.

Defining the Depth of the Impermeable Boundary
Because the depth-dependent K distribution in the fractured shale and depth to the effectively impermeable bottom boundary are unknown; treatment of deep flow needs to be simple and constrained by available information. We assigned the depth to the impermeable boundary through a water mass balance constraint applied to the transmissivity feedback calculations of flow during an annual cycle. The mass balance constraint involved estimating the net flow of water through the transect per unit width (transverse to the 10.1029/2019WR025093

Water Resources Research
flow direction) based on annual P minus annual ET, (P − ET) yr , applied over both the instrumented hillslope transect and the upslope contributing area. These two areas were assigned 137 and 345 m 2 , respectively, with the latter area obtained from an isosceles triangle with unit width base and its apex at the topographic peak ( Figure 1). Thus, the system consists of a narrowly diverging upper hillslope supplying the narrow parallelepiped-shaped lower hillslope. Note that this approach could be generalized to apply in other locations having simple diverging and converging domains (Salvucci & Entekhabi, 1995;Troch et al., 2003). The sum of these two surface areas for our system, 482 m 2 , was also used in a later comparison between predicted hillslope discharges and measured East River discharges. The net recharge through the transect's soil surface during WY2017* was then balanced with the sum of transmissive fluxes along the soil, weathering shale, and fractured bedrock layers over this same period by adjusting the thickness of the fractured bedrock. Because this adjustment requires extrapolating to depths considerably deeper than those of our instruments and K measurements, our deepest measured K values were assumed to be representative of even deeper fractured bedrock. Balancing recharge with the depth-integrated transmissivity flux yields an operationally defined depth for the impermeable boundary. Note that calibration of our system for the high precipitation WY2017* allowed comparison with transmissivity flux predictions in the very low WY2018. As shown later, such comparison is useful for guiding improvement in values of K used in the model.

Pore Water Sampling and Analyses
Pore water samples were collected at intervals ranging from weekly to approximately monthly, excluding winter months when collections were not possible, usually because of ice in aboveground tubing. Prior to sample collection, the air head space pressure was measured with a Tensimeter (Soil Measurement Systems) for determination of the matric/pressure head. Shallow pore water samples were collected via hand vacuum pump into a sampling flask, while deeper samples required pressure-lift using a bicycle pump. As noted earlier, pressure/vacuum soil water samplers do not yield instantaneous pore waters, but instead collect pore waters over varying intervals of time, limited by the hydraulic conductance of the ceramic and its surrounding medium (Tokunaga, 1992). Therefore, characteristic dates of the samples are earlier than their collection dates and were estimated based on response times of samplers under unsaturated versus saturated conditions (supporting information). Pore water specific conductance (SC) was measured with a VWR 61161 meter, and nonpurgeable organic carbon (referred to here as dissolved organic carbon, DOC) was measured with a Shimadzu TOC-VCPH analyzer.

Water isotopic analyses
The stable isotope compositions of groundwater deuterium δD and δ 18 O have been widely used to determine sources of the groundwater . Both snow (light isotopic compositions) and groundwater (heavier isotopic compositions) samples were collected for water analyses. A series of four snow pits were sampled in early March 2017 along a transect extending from above PLM1 (elevation of 2,900 m, Figure 1) down to the floodplain near PLM4, before the beginning of significant snowmelt. Groundwater samples for isotopic analyses were collected from the deeper pore water samplers at PLM1-3 on an approximately monthly basis from December 2016 through August 2018.
Water isotopic analyses were done at the Center for Isotope Geochemistry at Lawrence Berkeley National Laboratory using a Los Gatos Research Liquid Water Isotope Analyzer (LWIA). 1.2 ml of each sample was filtered and loaded into vials on the auto-sampler associated with the LWIA. Each sample was injected into the analyzer eight times, and the first three data points were discarded to avoid any effects resulting from carryover between samples. The remaining data points were averaged and corrected based on the results of internal standards with known isotopic compositions analyzed every third sample. Hydrogen and oxygen isotope ratios are reported using the conventional δ notation relative to Vienna Standard Mean Ocean Water.

Calculating Concentrations of Solute Discharges from the Hillslope
By combining transmissivity-based pore water fluxes described in section 3.5 with their associated solute concentrations, solute transport from the hillslope into the floodplain can be calculated. Specifically, where the F i (t) are the fractional transmissive fluxes in each zone relative to the total downslope transmissive flow, subscripts ws and fs refer to the weathering shale and fractured shale zones, respectively, and zonal 10.1029/2019WR025093

Water Resources Research
C i are different in the three approaches. Three different approaches were examined to represent solute contributions in the discharge waters. Each of these approaches is analogous to three-component mass balance models used to generate watershed C-Q relations (Chanat et al., 2002;Evans & Davies, 1998) but are based on different ways of treating measured subsurface concentrations. The first approach simply used the average of all data (SC or DOC) available within each zone (soil, weathering shale, fractured shale) to represent the zone-specific constant concentration. In the second approach, smooth curves were fit through the time trends in concentrations within each zone, and the time-averaged value of estimated concentration trends were taken to be representative of the constant zonal concentration. The third approach to estimating the composite concentrations in discharging hillslope pore waters attempts to be more realistic, through weighting time-dependent concentrations by their relative fluxes. Because of the sparse availability of pore waters from the soil, pore waters were not subdivided into surface soil and subsoil zones in any of the three approaches. Although overland flow was negligible at our site, equation (1) could be extended to include surface runoff when supporting measurements are available.

Soil Textures
The thickness of the soil was about 1.0 m along most of the hillslope transect (PLM0,1,2) and ranged from 1.2 to 1.8 m at PLM3. The underlying weathering shale could not be sampled by hand augering. Auger sampling of the floodplain site (PLM4) was limited to a depth of 0.80 m, where cohesionless, water-saturated sandy gravels were encountered. The hillslope soil textures were generally loam to silt loam, while the floodplain soils were silty clay and silty clay loam. Particle-size analyses show overall trends of slightly decreasing sand fractions and increasing silt and clay fractions in moving from PLM0 toward PLM4 (Table S1).

Hydraulic Conductivities
A total of 69 K measurements were made within the soil profiles at the hillslope and floodplain stations. A range of K in weathering shale at a depth of 1.40 m bgs near PLM6 was obtained based on Guelph permeameter fluxes at two water levels, using the nominal minimum and maximum diameters of the irregular borehole. The depth distribution of these K measurements in soils, weathering shale, and the two deeper piezometer and well K measurements are presented in Figure 4. For estimating subsurface flow along the hillslope, the hillslope soils were divided into shallow (0 to 0.50 m bgs) and subsoil (0.50 to 1.00 m bgs) zones, having average K of 9.7 × 10 −6 and 7.9 × 10 −6 m s −1 , respectively. The weathering zone shale and fractured

Water Resources Research
shale were assigned their average K values of 1.1 × 10 −5 and 1.6 × 10 −7 m s −1 , respectively. Impacts of different K values are examined in the Discussion section 5.3. Relative to the hillslope soils, the finer textured floodplain soils had distinctly lower K (Figure 4b). Also note from Figure 4b that average K values decrease slightly along the hillslope, but variability indicated by the standard deviation of the lognormal K is large.

Local, Vertical Hydraulic Potential Profiles and Time Trends Along the Hillslope
At PLM1, 2, and 3, pressure measurements on tensisamplers in combination with MPS-6 data yielded local matric/pressure head profiles essential for understanding flow within the hillslope subsurface. These data were also useful for constraining the depth above which vadose zone water is lost to ET and below which water is flowing toward the water table, that is, the zero-flux boundary, identified by the depth of maximum hydraulic head (Salve & Tokunaga, 2000;Wellings, 1984). Example hydraulic head profiles from the hillslope over the months following snowmelt obtained with this combined approach are shown in Figure 5. Given the large magnitude of the snowmelt accompanied with high ET during subsequent summer months, the zero-flux boundary generally follows the declining water table. It should be noted that even during WY2017, surface runoff was not observed along the hillslope, although localized surface flow was noted below the hillslope-floodplain interface. The effectiveness of ET in removing shallow soil water during summer months was reflected by lack of sustained downward hydraulic head gradients ( Figure 5) and very negative water potentials, despite precipitation occurring during the summer monsoon season. Thus, following snowmelt, the pore waters and dissolved solutes in the hillslope subsurface bifurcate into two directions at the zero-flux boundary. Fluxes below this boundary recharge baseflow, while fluxes from above the boundary are directed toward the atmosphere. After snowmelt, the zero-flux boundary extends progressively deeper ( Figure 5) until the following year's snowmelt event. Secondary hydraulic potential maxima and zero-flux boundaries can be generated by summer rainfall as illustrated in the 9-19-17 profile.

Soil Water Contents and Potentials in Late Summer
Water content and water potential profiles measured on soils collected on 14 and 15 September 2016, and on 22 September 2018 (only for PLM1, 2, and 3 in 2018) generally showed wetter conditions within 0.5 m of the surface, resulting from summer monsoon rainfall ( Figure 6). During the preceding 30 days, this amounted to about 70 and 46 mm of precipitation, for WY2016 and WY2018, respectively. At PLM0, the coarser textured soils (Table S1) and local abundance of Artemisia tridentata may explain the lower water contents and water potentials (Campbell & Harris, 1977;Ryel et al., 2003). The underlying soil along the hillslope was relatively drier, reflecting efficient evapotranspiration and lack of deep infiltration during summer. As noted earlier, the CLM-based calculations (Tran et al., 2019) for the summer months (WY2016, 2017, 2018) predicted that ET was large enough to prevent net infiltration of monsoon precipitation at this site. Although not part of the hillslope, the floodplain station PLM4 is included for contrast, with its higher water contents and potentials resulting from a shallow depth to the water table (0.76 m bgs during sampling) and higher clay content.

Water Table Depths and Time Trends Along the Hillslope
The depths to the water table along the hillslope over WY2017 and WY2018 are shown in Figure 7a, combining results from piezometers, tensisamplers, and moisture content sensors (5TE). The pressure transducerbased water elevations were periodically checked against manual water level measurements and corrected to within ±0.05 m. These data show slowing declining water levels at the beginning of each WY, followed by abrupt rises in April in response to warming temperatures ( Figure 2) and snowmelt. Throughout these annual cycles, the slope of the water table remains remarkably close to that of the soil surface so that the   (Figure 7b). The uncertainty in individual water levels (PLM1 and PLM6) introduced an uncertainty of less than 1 mm m −1 in the overall gradient. Recall that the characteristic K of the hillslope soils and weathering shale are relatively high ( Figure 4) and therefore can accommodate fairly rapid flow from snowmelt infiltration. Indeed, overland flow along the hillslope was not observed, although surface flow was noted along parts of the floodplain during snowmelt.

Pore Water Chemistry
Pore waters in the unsaturated zone and groundwater were collected through both WY2017 and WY2018, with the water table annually fluctuating between the upper fractured shale zone and reaching into the soil. Weather conditions (freezing of water in sampling tubes) and extended periods of very negative matric potentials in soils and the upper weathering zone prevented year-round collection of shallower pore water profiles. However, fairly comprehensive pore water sample sets were obtained within the three hillslope monitoring stations at and shortly after the end of snowmelt (May and June). A few of these more complete profiles, and less complete late summer profiles of salinity (as SC) and dissolved organic carbon (DOC), are presented in Figure 8. Highest salinities generally occur within the 1.0 to 3.6 m depth interval, consistent with weathering release of major ions from Mancos Shale (Tuttle et al., 2014a(Tuttle et al., , 2014b, and leaching of ions from the soil profile during snowmelt. In contrast, DOC levels are commonly high within the soil (when matric potentials are sufficiently high to allow pore water sampling), reflecting the dominant influence of plant/soil/microbial interactions in organic matter production and release.
Time trends of pore water SC and DOC ( Figure 9) from all of the hillslope stations and depths show considerable variability, but distinctions between the soil, weathering, and fractured shale zones persist. It should

10.1029/2019WR025093
Water Resources Research be noted that times associated with pore water samplers were based on a response correction dependent on proximity to the water table (supporting information). Soil SC (Figure 9a) range from being moderately elevated to very low, with the latter condition reflecting dilution and leaching from snowmelt. Highest levels of SC generally occur within the weathering zone, which experiences snowmelt dilution to a lesser extent than the soil. Lowest SC values occur within the fractured shale, with negligible temporal variation.
For later calculations of solute transport, continuous curves representing the time dependence of representative SC in the soil and weathering zones were generated by fitting through points having at least three samples on a given date and by assuming that SC values plateau during summer through winter at the maximum measured average value. The fit curve was generated by piecewise scaling the error function to declining and rising intervals as described in supporting information. Based on the abundant measurements of fractured shale pore waters, its SC was assumed to remain constant and set at the average of all values obtained within this zone. Here, it is worth noting that the constant zonal SC approximation commonly employed in hydrograph separation (Stewart et al., 2007) is supported only in our measurements of fractured shale pore waters, while snowmelt-driven dilution and displacement impart clear time dependence to SC in the soil and weathering zones.
For DOC, large annual fluctuations in concentrations occurred only within the soil, where sharp decreases accompanied snowmelt infiltration ( Figure 9b). Thus, approximate time dependent soil DOC behavior was generated in a manner similar to that described previously for SC, while the averaged DOC value was used as its constant concentration in the fractured shale. The DOC concentration trends in the weathering zone exhibited a decline during the WY2017 snowmelt, followed with further gradual decline through WY2018. Therefore, the error function fit was applied only to the initial WY2017 decline, followed with a linear fit through the remaining data. Representative zonal SC and DOC obtained by averaging the available data and the curve fits are indicated along the concentration axes. The higher values obtained by averaging the estimated time dependent is a consequence of expected persistence of higher concentrations from summer through winter.

Water Isotopes
The weighted average isotopic composition of snow from all four snow pits was −20.0‰ for δ 18 O and −148‰ for δD, with only a 1‰ range in δ 18 O and 3‰ in δD between the average values for the pits. The isotopic composition of the summer monsoon rain is much less well constrained but is generally much heavier (very approximately averaging −10‰ for δ 18 O and −70‰ for δD).  Since the pioneering study by Dinçer et al. (1970), δ 18 O values of fresh snow and average δ 18 O values of snowpacks have often been used as end members for calculations of snowmelt recharge into groundwater (Taylor et al., 2002). However, isotopic fractionation of the snowpack and snowmelt over time is significant, causing enrichments of the heavy isotopes in the recharging meltwater (Taylor et al., 2002). The significant and consistent decreases in δ 18 O over time in the groundwater within the fractured shale reflect the influence of the large WY2017 snowmelt recharge event ( Figure 10). The general decreases in δ 18 O with depth are consistent with earlier stage snowmelt percolating deeper into the fractured shale bedrock. The exception to this trend found in the 6.89 m depth pore waters at PLM1 may result from complex flow paths occurring in this fractured shale. Given the high contrast between these groundwater δ 18 O values and those of rainwaters, the water isotopes indicate that rainfall contributions to the hillslope's groundwater are very small. As noted above, other hydrologic considerations also show that monsoon precipitation did not contribute to groundwater recharge, even during WY2017.

Transmissivity Feedback Calculations of Subsurface Flow
The time-dependent water table measured within the monitored transect (Figure 7a, "piezometer average" curve) was used to determine the time-dependent transmissivities of the saturated portions of soil, weathering shale, and fractured shale. The average of measured K in each zone (Figure 4) was multiplied by its timedependent saturated thickness to obtain its time-dependent transmissivity. Note from Figure 7a that the soil transmissivity with respect to downslope flow is only active when the water table rises above 1.0 m bgs to saturate varying thicknesses within the soil, and that the transmissivity of the weathering shale is variably active throughout the year. In contrast, the transmissivity of the underlying fractured shale remains constant throughout our study interval because the water table does not drop below its interface with the weathering shale. The water table-dependent downslope fluxes within effectively saturated soil, weathering zone, and fractured shale is qualitatively depicted in the upper panels of Figure 11.
Based on our measurements of water level recovery in wells and piezometers, we used a K of 1.6 × 10 −7 m s −1 to represent the fractured shale ( Figure 4a) and assign the depth to the impermeable boundary through a water mass balance constraint. The continuously measured gradient (Figure 7b), only varied between 0.189 and 0.207 m m −1 , remaining close to the overall slope between PLM1 and both PLM3/PLM6 (0.200 m m −1 ).
With these approximations, daily subsurface flow was calculated from 12-1-2016 to 12-1-2017 (WY2017*), with total P and ET of 872 and 310 mm, respectively, with the thickness of the fractured shale layer

10.1029/2019WR025093
Water Resources Research adjusted to provide a net flow through all three layers equal to 562 mm (P − ET). Adjustment to accommodate this net recharge resulted in a calculated thickness of 172.3 m for the fractured shale and depth to the effectively impermeable boundary at 176 m bgs. Based on these values, the time-dependent contribution to baseflow fluxes from the fractured shale, weathering shale, and soil layers were calculated. It is emphasized that a constant K within the fractured shale was assumed in order to obtain this estimated thickness and that the transmissivity of this deeper zone is the more general parameter to vary. Thus, trends of net decreases or increases in fractured shale K with depth would lead to deeper or shallower depths to the assigned impermeable layer, respectively.
As shown in Figure 11b, the calculated fluxes within the fractured shale were nearly constant throughout the year, with small fluctuations that reflected the time-varying hydraulic gradient (water level differences between PLM1 and PLM6, Figure 7b). In contrast, contributions of flow though the weathering shale and soil layers are highly seasonal, reflecting annual snowmelt and ET responses. Also shown in Figure 11b are the calculated total subsurface fluxes for WY2017* and WY2018 normalized to the combined transect Nevertheless, snowmelt recharge through soil is certainly generally important in supplying water into the weathering shale and fractured shale layers, even in years with low precipitation.

Solute Discharges From the Hillslope
The measured concentrations of major ions (represented by SC) and DOC ( Figure 9) were combined with their corresponding transmissivity-based pore water fluxes through the soil, weathering shale, and fracture shale zones depicted in Figure 11b to estimate the time-dependent exports of these solutes from the hillslope into the floodplain. As mentioned in the Methods, the solute concentrations were applied in three different ways to their corresponding fractional fluxes (Figure 11c) to generate flow-weighted concentrations in the discharging pore waters. The first approach to calculating discharge concentrations simply used the average of all data (SC or DOC) available within each zone (soil, weathering shale, and fractured shale) to represent zone-specific constant concentrations. For the soil and weathering zone, these average concentration values are indicated by the dashed arrows along the concentration axes of Figures 9a and 9b. Inspection of Figure 9 shows that this procedure is satisfactory only for the fractured shale. The second approach estimated timeaverage concentrations of fit curves (solid arrows along the concentration axes of Figure 9), and these are higher in the soil and weathering zones because the less frequently recovered summer through winter pore waters are expected to have higher solute concentrations. Although the curves shown in Figure 9 are qualitative, they capture the important influence of dilution during and shortly after snowmelt recharge events and the typically more concentrated (yet undersampled because of lower matric potentials) pore waters associated with summer through winter. The third approach to estimating the composite concentrations in discharging hillslope pore waters attempts to be more realistic, through weighting time-dependent concentrations by their relative fluxes. This last approach amounted to weighting the continuous curves in Figure 9 by their associated fractional fluxes (Figure 11c).
The predicted time-dependent discharge concentrations obtained via these three approached are shown in Figures 12a and 12b for SC and DOC, respectively. For discharges of major ions, all three approaches predict that peak concentrations occur prior to the calculated time of maximum subsurface discharge, reach local minima during and following maximum flow, and annually reach minimum concentrations just prior to initiation of snowmelt. The predicted maxima in discharge SC ahead of maximum pore water flow is a consequence of highest SC being associated with the weathering zone and mobilized during water table rise, while further increases in the hillslope water table lead to dilution through downslope fluxes of less saline soil water. When average SC values from each zone are applied, the depressions in SC associated with peak flow are smaller in magnitude and duration relative to that obtained with the time-dependent SC relation. In addition, secondary maxima in discharge SC are predicted during late summer to fall when the estimated

Water Resources Research
time dependence of SC is included in calculations. It should be noted that the secondary maxima reflect assumed sustained elevated SC in the weathering zone which was based on measured increasing SC trends during summer months (Figure 9a).
Concentrations of DOC in pore waters discharging from the hillslope are predicted to reach maximum values during or shortly ahead of the annual maxima in subsurface flow, followed by steady decline until the next snowmelt cycle (Figure 12b). When average zonal DOC values were used, the rise and fall of discharge DOC concentrations occurs in synchrony with water levels, and this simply reflects the fact that highest average DOC concentrations are associated with the soil, which is very transmissive when saturated. When the time dependence of DOC concentrations within different zones is included, the concentration maxima occur earlier than the flux maxima, reflecting advection of previously elevated DOC levels immediately before influences of dilution by snowmelt. The generally simpler behavior predicted for DOC relative to that of SC follows from the fact that the distinctly higher DOC concentrations in soil undergo dilution during and immediately after snowmelt, while DOC concentrations in the weathering zone and fractured shale are similar to each other. Because of the latter condition, switching between periods with flow dominated by the weathering zone and fractured shale (Figure 9c) has less impact on discharge DOC concentrations than it does on discharge SC.
The relations presented in Figure 12 can be replotted to depict the dependence of predicted discharge concentrations on overall subsurface flow rates, leading to subsurface C-Q relations shown in Figure 13.
When average values of SC and DOC are applied in each zone, predicted discharge concentrations are solely determined by these averages and the water table elevation. Therefore, the WY2017 and WY 2018 C-Q obtained with average zonal values overlap, without hysteresis. In both the averaging and time-dependent treatments, discharge SC reach maximum values at intermediate flow rates because highest SC reside within the weathering zone, while discharge DOC reach maxima at higher flow rates that support downslope export of DOC-rich soil pore waters. Most of the C-Q relations obtained when time-dependent concentrations are flux-weighted exhibit clockwise hysteresis as well as year to year variation. The hillslope subsurface DOC showed flux-weighted C-Q with increasing and wide clockwise hysteresis during WY2017 because of the significant period of downslope flow through the DOC-rich soil zone.

Applicability of Transmissivity Feedback
Relative to the fractured shale, K of the weathering shale and soil zones increase by over an order of magnitude ( Figure 4) and therefore substantially enhance transmissivity when the water table rises in response to snowmelt. The increased K in shallower zones along with the measured hydraulic head trends all clearly support application of transmissivity feedback during snowmelt. Although flow along transient perched saturated zones has sometimes been observed (Dunne & Black, 1970;Leung et al., 2011;McDaniel et al., 2001) or assumed at other sites, the measured hydraulic potential profiles ( Figure 5) indicate that perched water is unimportant on this hillslope, even during the 2017 snowmelt event. Neither groundwater ridging (Sklash & Farvolden, 1979) was observed in this system either nor was it expected because of the strong

Water Resources Research
topographic influence that keeps the water table approximately parallel to the hillslope surface during its annual oscillations (Figure 7). While the overall depth trends in our K profiles justifies use of transmissivity feedback, significant uncertainties remain with respect to suitable values of K and their impacts on predicted discharge concentrations.

Impacts of the K Profiles and its Variation
Given the relatively uniform thickness of the hillslope soils and the practically constant slope of the water table (Figure 7b), the greatest uncertainties associated with predicting subsurface flow are associated with values of K assigned to the soil (K s ), weathering zone (K ws ), and fractured shale (K fs ). The importance of K distributions and their uncertainties in controlling predictions of subsurface flow are well recognized (Freeze & Witherspoon, 1967;Gleeson & Manning, 2008;Toth, 1963;Welch & Allen, 2014). Although a relatively large number of K measurements were obtained within the soil zone, and downslope transmissive flow (as opposed to infiltration into the weathering zone) occurs over relatively short periods of time within the soil, previous studies suggest that the K s and K ws suitable for hillslope-scale flow are greater than that obtained from permeameter measurements (Brooks et al., 2004). We address the uncertainty associated with K fs in the next section in the context of defining the depth to the impermeable boundary. Here, we first examine impacts of changing K ws and then also increase K s .
Recognizing the potentially wide range of possible K values, order of magnitude increases and decreases in K ws relative to the previously applied average value (1.1 × 10 −5 m s −1 ) were initially examined. However, the maximum possible K ws is only about 2.75 times greater than the reference value because the mass balance constraint would require that the impermeable boundary be set at our existing borehole depths of 10.0 m bgs, and this depth lacked evidence of reduced permeability. Thus, the upper limit of K fs was set to 1.65 × 10 −5 m s −1 , only 50% greater than the reference value. With this upper limit K ws , the depth of the impermeable boundary is at 128.7 m bgs. For exploring changes in the opposite direction, we lowered K ws to 3.67 × 10 −6 m s −1 , one third of the reference value, and this required setting the impermeable boundary at 239.0 m bgs in order to maintain mass balance.
Following the above procedures, the influences of increasing and decreasing K ws on predicted subsurface flow are shown in Figure 14a. The dotted, solid, and dashed curves represent cases with K ws equal to 3.67 ×10 −6 , 1.1×10 −5 , and 1.65×10 −5 m s −1 , respectively. Under higher K ws , fluxes through the fractured bedrock must be reduced to maintain mass balance. As shown in Figure 14a, the higher K ws generates higher maximum fluxes in the weathering zone and higher peak flows from the hillslope, offset by lower minimum fluxes during winter baseflow. The latter effect results from baseflow occurring through a shallower fractured bedrock (lower transmissivity) relative to the reference case. Conversely, a lower K ws dampens the maximum fluxes through the weathering zone and suppress peak flows along the hillslope, while increasing winter baseflow rates through a thicker fractured bedrock. The lack of K ws induced changes in transmissive fluxes through the soil in Figure 14a is simply a consequence of not changing K soil .
Unlike the restricted potential increases in K ws , large increases in K s are allowable without violating water mass balance because periods of water table residence within the soil zone are relatively short. Therefore, the example of combining effects of increasing K ws and K s shown in Figure 14b involves the same 50% increase in K ws , and a 25-fold increase in K s relative to their respective reference values. The water mass balance constraint with these K values requires setting the impermeable boundary at 29.8 m bgs. While a 25fold increase in K s appears large, a comparison of hillslope plot scale K s and Guelph permeameter-based K s measurements yielded ratios ranging from 15.1 to 48.7 (Brooks et al., 2004). These changes to K ws and K s greatly increase the predicted fluxes along the weathering zone, soil, and total fluxes during snowmelt, while decreasing the overall fluxes through the fractured bedrock and total hillslope fluxes during winter baseflow.
Recall that the water mass balance constraint was applied for WY2017* to calibrate the effective depth of the hillslope, using our measurements-based K values. Thus, the fluxes for WY2018 are predictive. Assuming that water storage within the system is nearly the same at the beginning of each WY, the predicted WY2018 subsurface fluxes obtained with different K values can be compared against the estimated P − ET of 270 mm for WY2018. Under the "standard" scenario, the WY2018 discharge was equivalent to 523 mm (Figure 11b), nearly double the estimated P − ET for the WY. With the example enhanced K in the weathering zone and soil, the WY2018 subsurface discharge amounts to 311 mm, much closer to the estimated P − ET. This comparison suggests that our Guelph permeameter-measured K in fact substantially underestimates downslope flow at the field scale. Figure 14. Impacts of different K values on the predicted water fluxes, discharge SC, and discharge DOC. Time trends for (a) influences of different K ws on subsurface fluxes within different zones, (b) influences of increasing both K s and K ws on subsurface fluxes within different zones, (c) subsurface flow-weighted SC of discharging pore waters, and (d) subsurface flow-weighted DOC concentrations in discharging pore waters. The K for fractured shale bedrock was kept constant, while the thickness of the fractured shale was adjusted to achieve water mass balance. Vertical dashed lines indicate dates associated with calculated annual peak subsurface discharge rates.

Composite Solute Concentrations in Hillslope Discharge Waters and Influences of K Variations
The most distinct features of the solute concentration time trends are the dilution of major ions (SC) and DOC from both the soil and weathering shale along the hillslope during snowmelt leaching (Figure 9), similar to trends measured in other catchments, whether recharged by snowmelt (Boyer et al., 1997(Boyer et al., , 2000Hornberger et al., 1994) or rainfall (Anderson & Dietrich, 2001;Herndon et al., 2015). In contrast to these other catchments, the magnitudes of major ion exports from this East River lower montane hillslope are substantially higher, as expected from the salt-rich Mancos Shale (Tuttle et al., 2014b). It should be noted that the hysteresis loops for the SC C-Q (Figure 13a) do not conform to any of the diagnostic shapes categorized in Figure 2 and Table 1 of Evans and Davies (1998) because of salinity dilution and reordering of relative concentrations within zones during snowmelt recharge (Figure 9a). The WY2017 DOC C-Q loop conforms roughly to Type C2 of Evans and Davies (1998) because the relative concentration order (decreasing DOC concentration with depth) is preserved, even during and after snowmelt (Figure 9b). In contrast, the relatively small variations in DOC discharge concentrations (nearly chemostatic) calculated for the low flow WY2018 results from the water table remained too low to export DOC along the soil zone and the DOC concentrations in the weathering zone becoming very similar (Figure 9b).
The impacts of varying K ws on predicted changes in exports of major ions (SC) and DOC in the discharging pore waters are shown in Figures 14c and 14d, respectively. Increased K ws (dashed curves) generally results in higher concentrations of these solutes because fractional contributions of more dilute pore waters in the fractured bedrock must be decreased to obtain water mass balance. All maxima in SC are amplified when K ws is increased because saturated flow is always present to varying extents within the solute-rich weathering zone. In contrast, the DOC concentration maxima were only amplified during snowmelt in WY2017, when water table rise extended into the DOC-rich soil zone. When K ws was lowered, predicted export concentrations of SC and DOC both were reduced because of decreased flow through the solute-rich weathering zone and compensating increased flow through the more dilute fractured bedrock. Lastly, when the 25-fold increases in K s were applied, further enhancement of peak discharge concentrations are predicted for SC and DOC (red curved in Figures 14c and 14d, respectively), again because of the higher concentrations of these solutes in soils relative to bedrock. The notable exception is seen for SC during peak water fluxes, reflecting strong dilution from snowmelt.

Defining the Depth to the Impermeable Boundary
Recall that the impermeable boundary at the bottom boundary of this hillslope was assigned based on the subsurface flow mass balance constraint described in sections 3.6 and 4.8. A key simplification applied in this analysis is that the average K of the fractured bedrock obtained on measurements that extend only to

Water Resources Research
10 m bgs is representative of the deeper bedrock until an impermeable boundary is encountered. Given that the calculated depth to the impermeable boundary under our standard conditions is 176 m bgs, far deeper than any of our current measurements, further drilling to allow deeper K measurements, pore water sampling, and pressure monitoring is clearly warranted. During recent drilling down to 70 m bgs between PLM1 and PLM2, fractures were commonly present along the full length of the core samples. The fact that an impermeable boundary is required at 29.8 m bgs when the likely more realistic higher K ws and K s are used, indicates that K ws must be decreasing at depths greater than 10 m. Thus, our impermeable boundary is an operationally defined and can be improved upon when deeper K measurements are available.

Comparison with East River Discharges and Concentrations
Continuous measurements of East River discharge and SC recorded at the Pumphouse gauging station, located 0.3 km away from the hillslope transect (Figure 1), are very helpful for understanding the lower montane hillslope's hydrogeologic characteristics within the context of the broader watershed. Measurements from that station are part of a detailed analysis on hydrologic responses of multiple catchments within the upper East River watershed (Carroll et al., 2018;Winnick et al., 2017). The Pumphouse station integrates flow and transport from 85 km 2 of the East River watershed, over terrain that rises to elevations that are nearly 1.4 km above that of the lower montane hillslope transect. Relative to the lower montane hillslope, snowpack at higher elevations is greater, especially within the upper subalpine regions of the watershed (Carroll et al., 2018). The East River's hydrograph and SC time trends from the Pumphouse station show large decreases in SC resulting from snowmelt in both years, and much larger discharges in WY2017 (7.30 ×10 7 m 3 ) relative to WY2018 (3.45 ×10 7 m 3 ) (Figure 14a), while river C-Q relations for SC show dilution with increased flow and decreasing clockwise hysteresis reflecting the importance of runoff and dilute surface pore waters during snowmelt (Figure 14b). Earlier C-Q relations obtained at the Pumphouse station (2014 to 2015) for major ions and DOC exhibited clockwise decreasing and increasing hysteresis, respectively (Winnick et al., 2017).
Comparisons between the SC of hillslope pore waters (Figures 9a, 12a, 13a, and 14c) and East River waters (Figures 15a and 15b) show that the hillslope discharges major ions at substantially higher concentrations than most other areas upstream of the Pumphouse. The flow-weighted average SC was calculated by summing the product of daily flow times SC over the year and dividing by the annual flow. The SC associated with hillslope discharges are about five to six times greater than that of the larger watershed (Table 1).
Recall that the bedrock underlying the hillslope transect consists of Mancos Shale, which is known to release solutes at elevated concentrations. While Mancos Shale also occurs over many portions of the upper East River watershed, this larger area also contains an abundance of igneous rocks and sandstones (Carroll et al., 2018;Gaskill et al., 1991) that generally release solutes at lower rates (Bluth & Kump, 1994). Overland flow at the higher elevation regions of the watershed is also important in lowering SC in the river.
Comparisons of flow-weighted SC from high snowpack WY2017 and low snowpack WY2018 may help predict more general discharge responses to changing climate. With lower snowpack and snowmelt, discharges

Water Resources Research
to the river are more strongly influenced by deeper flow. Hence, vertical solute profile associated with different settings can strongly influence overall discharge chemistry. In the context of transmissivity feedback, melting of low snowpack is less effective at direct transport of near-surface solutes. At the lower montane hillslope, lower inventories of salts (SC) and DOC were transported downslope during WY2018. Inspection of Table 1 shows that the calculated hillslope average discharge SC decreased from 1,172 μS cm −1 in WY2017 to 1,139 μS cm −1 in WY2018. In contrast, the average SC increased from 203 μS cm −1 in WY2017 to 232 μS cm −1 in WY2018 over the larger Upper East River watershed, indicating greater influences of leaching of solutes from deeper bedrock following low snowmelt in the higher elevation catchments.
Surface area-normalized estimates of flow contributions to river discharge from the hillslope relative to the larger watershed provide another assessment of the hydraulic behavior of the lower montane region. We compared the net flow of water per unit area for WY2017 and WY 2018 by normalizing discharges from the hillslope and the East River above Pumphouse to their respective areas ( Table 1). The lower values for the hillslope relative to the upper East River watershed are consistent with its lower precipitation (especially as snow). Although surface runoff was not important along the lower montane hillslope, it is in the upper watershed with its higher discharge.

Generalization
The approach we have developed is expected to be broadly useful at other hillslope sites where recharge occurs primarily through snowmelt, provided that key information is available and a few commonly occurring conditions are met. The approach requires estimates of net recharge (P − ET) and the contributing upslope surface area estimated from topography. Within the lower section of the hillslope, it requires measurements of K profiles, hydraulic potential profiles (or at least water table depths), and solute concentration profiles. Applicability of transmissivity feedback also requires that the K profile generally increases toward the soil surface, a condition that is often met. Within the lower section of hillslopes, the approach also requires two elevationally separated monitoring stations with vertically resolved measurements, in order to track hydraulic gradients. In some settings including our study site, the water table remained nearly parallel to the soil surface, such that paired stations may not be needed for determining the hydraulic gradient. Although ET estimates are often uncertain in mountainous terrain (Henn et al., 2018), the recent approach by Tran et al. (Tran et al., 2019) provided predictions that are remarkably consistent with soil mass balance measurements. Lateral variability of flow and transport transverse to the mean flow direction can limit predictions based on single transect data (Woods & Rowe, 1996). With these uncertainties in mind, the strategy developed here to calculate subsurface flow and transport may be applied to multiple representative hillslope environments within a watershed to elucidate their individual subsurface contributions to river recharge.

Conclusions
Hydrologic and geochemical measurements showed that the subsurface response to snowmelt along a lower montane hillslope exhibited distinctly different characteristics during the high-and low snowpack WY2017 and WY2018. Above-average snowmelt in WY2017 resulted in nearly fully saturated soils along the hillslope and efficient mobilization of DOC in soils zone and salts within the weathering Mancos Shale. The weathering shale zone maintained highest levels of soluble salts, and the hillslope was generally higher in salinity than most of the upper East River watershed, reflecting the dominant influence of Mancos Shale at this site. In WY2018, early melting of the low snowpack mobilized smaller inventories of salts and negligible DOC. The basic subsurface flow behavior was described using the transmissivity feedback approach, and in combination with water mass balance, we estimated the effective maximum depth of flow at 176 m bgs. If higher the hillslope-scale effective K of the soil and weathering shale are substantially larger than indicated by borehole measurements, the effective maximum flow depth is shallower. Multiple field soil moisture measurements and CLM calculations indicated that ET prevents summer monsoon precipitation from contribution to recharge at this site. The δ 18 O values of the hillslope groundwater show that it is largely derived from snowmelt and that early snowmelt percolates deeper into the fractured bedrock. Inclusion of estimated time-dependent solute concentrations within each zone yielded predictions of hysteretic subsurface C-Q relations that can help explain hillslope exports to the river. The mass balance transmissivity feedback-based approach permits spreadsheet calculations that are further simplified for snow-dominated mountainous environments where the water table response is dominated by the annual spring snowmelt event. This approach can be used to guide the design of more complex models and can be refined as additional field data become available. Finally, these measurement-based calculations can be used to easily explore hillslope responses to different climate (snowmelt) scenarios, as demonstrated in our analyses of subsurface flow and transport for the very different WY2017 and WY2018. drilling, and Craig Goodknight (Navarro, Inc.) for field logging of drill cuttings, and Jose Suarez (HRL Compliance Solutions) for drilling. We gratefully acknowledge Billy Barr (Rocky Mountain Biological Laboratory) for the Gothic weather and snow data. We appreciate the valuable comments provided by the reviewers and the Associate Editor, which helped improve presentation of this work. Data used in this paper are deposited in the U.S. DOE Environmental Systems Science Data Infrastructure for a Virtual Ecosystem (ESS-DIVE).