Spatial controls on groundwater response dynamics in a snowmelt-dominated montane catchment

The role of spatial variability in water inputs on runoff dynamics has generally not received as much research attention as topography and soils; however, the influence of topography and forest cover on snow surface energy exchanges can result in asynchronous snowmelt throughout a catchment, complicating the space–time patterns of runoff generation. This study investigates temporal variation in the relative importance of spatial controls on the occurrence, duration, and timing of shallow groundwater response, utilizing a highly distributed monitoring network in a snowmeltdominated montane catchment in western Canada. The study findings indicate that deep-soil hydraulic conductivity is a first-order control on the spatial distribution of sites that generate shallow groundwater response versus sites that experience only deep percolation. Upslope contributing area and slope gradient are first-order controls on the duration of groundwater response during peak-flow, recession-flow, and low-flow periods. Shallow runoff response areas expand and contract throughout these periods and follow the general spatial patterns of topographic convergence. However, spatial controls on the timing, intensity, and quantity of snowmelt and controls on vertical versus lateral flux partitioning in the soil overwhelm the influence of topographic convergence on runoff patterns during early spring freshet periods. The study findings suggest that various topographic indices and topography-based rainfall runoff models would not likely be good predictors of runoff patterns in snowmelt-dominated montane catchments during early phases of the spring freshet, but would increase in importance as the freshet and post-freshet periods proceed.


Introduction
Understanding runoff generation processes and streamflow source area dynamics is important for predicting streamflow quantity, quality, and timing, and for assessing the potential impacts of land use and climate changes on water resources (Beschta et al., 2000;Stewart, 2009;Stewart et al., 2005).As a consequence, runoff generation has been a prominent research theme in hydrology for decades, with much greater focus on rainfall runoff than snowmelt runoff.For rainfalldominated catchments, conceptual models of runoff dynamics have typically emphasized the influences of topography and soil characteristics on the downslope flow of water, particularly in relation to flow convergence, connectivity of hillslope flow paths, and threshold responses (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000;Tromp-van Meerveld and McDonnell, 2006a, b;Tromp-van Meerveld et al., 2007;Penna et al., 2011).For example, the hydrogeomorphic concept articulated by Sidle et al. (2000) focuses on the activation of different hydrogeomorphic units as a function of catchment wetness.The storage-excess runoff (i.e.fill and spill) concept similarly examines runoff generation in relation to the effects of soil wetness and storage capacity on flow path continuity (Spence and Woo, 2003;Tromp-van Meerveld and McDonnell, 2006a, b).
Most runoff in montane catchments is generated via subsurface flow, particularly through matrix or macropore flow within saturated soils, and via saturation-excess overland flow or return flow, which are dependent on soils being saturated to the surface (Buttle, 1994;Buttle et al., 2004;McGlynn et al., 1999;Sidle et al., 2000;Sklash and Farvolden, 1979).Infiltration-excess overland flow is rare in montane catchments due to generally high infiltration capacities relative to maximum water input intensities.Exceptions include disturbed sites such as logging roads and locations where soil freezing reduces the infiltration capacity due to the presence of ice-filled soil pores (Dunne and Black, 1971;Laudon et al., 2004;Stadler et al., 1996;Stein et al., 1994).At many montane sites, soils are relatively shallow, highly permeable, and are underlain by relatively impermeable bedrock or glacial till (Freer et al., 2002;Kim et al., 2004;McGlynn et al., 1999;Sidle et al., 2000).Under these conditions, saturated zones form above the confining basal layer and topographic indices have been found to be effective for predicting the spatial patterns of soil saturation, hydrologic connectivity, and runoff generation (Thompson and Moore, 1996;Beven and Kirkby, 1979;Freer et al., 2002).However, at sites with deeper soils, transient saturated zones can form within the soil at depths where the rate of downward percolation exceeds the ability of the soil's hydraulic conductivity to allow drainage, resulting in percolation-excess runoff generation (Latron and Gallart, 2008;Redding andDevito, 2008, 2010).
The role of spatial variability in water inputs (rainfall, snowmelt, or both) on runoff dynamics has generally not received as much attention as topography and soils, particularly at the scale of headwater catchments -where much of the research on rainfall runoff processes has been conducted.While mountainous topography can significantly influence the spatial distribution of rainfall (Goodrich et al., 1995;Guan et al., 2005;Linderson, 2003;Shoji and Kitaura, 2006;Hrachowitz and Weiler, 2011), some studies indicate that the spatial variability of rainfall decreases with increasing event magnitude (Linderson, 2003;Taupin, 1997).In contrast, snowmelt inputs can exhibit significant and systematic spatial variability, even in headwater catchments, due to the influences of slope, aspect, elevation, and forest cover on both snow accumulation and snowmelt processes (Balk and Elder, 2000;Berris and Harr, 1987;DeBeer and Pomeroy, 2010;Ellis et al., 2011;Jost et al., 2007;Toews and Gluns, 1986;Winkler et al., 2005).In regions where melt is dominated by radiation, seasonal melt begins earlier and typically occurs at higher rates at sites with high insolation (e.g.south-facing sites in the Northern Hemisphere) (DeBeer and Pomeroy, 2010;Ellis et al., 2011;Jost et al., 2007Jost et al., , 2012;;Toews and Gluns, 1986) and at sites lacking shading from forest cover (Anderton et al., 2002;Daly et al., 2000;Hock, 1999;Marks et al., 2002;Winkler et al., 2005).Where melt is dominated by the turbulent fluxes of sensible and latent heat, open sites with higher wind speeds experience higher melt rates than those with forest cover (Berris and Harr, 1987).The influence of topography and forest cover on snow surface energy exchanges can result in desynchronization of snowmelt throughout a catchment (Boyer et al., 2000;DeBeer and Pomeroy, 2010;Ellis et al., 2011;Jost et al., 2007Jost et al., , 2012) ) complicating the space-time patterns of runoff generation.Further complexity arises because some physiographic variables can exert contrasting influences on snowmelt runoff.For instance, sites with high insolation might experience more evapotranspiration throughout the growing season and more sublimation throughout the winter season, resulting in drier soils and less snow prior to spring melt.However, the same sites might also experience more rapid snowmelt due to greater energy inputs and, therefore, potentially more rapid runoff response once soils are sufficiently wet.On the other hand, the influences of drier antecedent soil wetness and less snow as water input to the soil could overwhelm the influence of higher melt rates in some circumstances.Because of the differences in water input dynamics between rainfall and snowmelt processes, it cannot be assumed that the existing conceptual models of rainfall runoff (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000) appropriately represent runoff dynamics in snowmelt-dominated mountainous catchments.In particular, topographic and geologic controls on flow path convergence and hydrologic connectivity might be of lower importance than meteorological controls on water input patterns in determining runoff dynamics during snowmelt.
Since runoff generation in forested catchments typically depends on the development of phreatic conditions within the soil, understanding groundwater (used here to refer to phreatic water regardless of the depth below the soil surface) dynamics is important for understanding runoff dynamics (Anderson and Burt, 1978;Jencso et al., 2009;Kuras et al., 2008;Monteith et al., 2006a, b;Seibert et al., 2003).Among studies that investigated groundwater-related runoff dynamics in snowmelt-dominated forested catchments, most have been conducted in small catchments (e.g.0.3-50 ha) with limited elevation ranges (e.g.20-200 m of relief) (Buttle et al., 2004;Dunne and Black, 1971;Flerchinger and Cooley, 2000;Laudon et al., 2004;McDaniel et al., 2008;McNamara et al., 2005;Monteith et al., 2006a, b;Seibert et al., 2003), which would have limited the spatial variability in the timing, quantity, and intensity of snowmelt water inputs and associated impacts on runoff generation patterns.Larger or higherrelief catchments with complex terrain and variable land cover experience large gradients in meteorological and snowpack conditions that could generate asynchronous snowmelt, leading to isolated areas of groundwater response.context of asynchronous water inputs that can occur under snowmelt conditions (Boyer et al., 1995(Boyer et al., , 1997(Boyer et al., , 2000;;Deng et al., 1994;Hinckley et al., 2014;Kuras et al., 2008) and several of these focused more on the flushing of dissolved organic carbon than on runoff generation processes (Boyer et al., 1995(Boyer et al., , 1997(Boyer et al., , 2000)).
The current study addresses spatial controls on groundwater dynamics and their implications for runoff generation in the Cotton Creek Experimental Watershed (CCEW), a snowmelt-dominated montane catchment in southeastern British Columbia, Canada, with complex terrain and variable forest cover.Unlike many other montane study sites that have relatively shallow soils with maximum depths of 1 or 2 m (Sidle et al., 2000;McGlynn et al., 1999;Freer et al., 2002), CCEW is mantled by deep glacial tills in excess of 8 m in some areas.The objectives of the study are (1) to determine what physiographic properties are dominant in controlling the spatial distribution of groundwater response occurrence, duration, and timing accounting for different periods of the annual streamflow hydrograph and within the context of asynchronous water inputs during snowmelt conditions; and (2) to determine whether the spatial patterns of response are consistent with the focus of many rainfall runoff conceptual models that topographic convergence is a dominant influence in controlling if, when, and how much a site contributes to runoff generation in relation to changing catchment wetness (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000).These objectives are addressed by fitting statistical models to explain the spatial variation in groundwater dynamics as a function of landscape properties for hydrologically distinct periods of the annual streamflow hydrograph.
The UEC catchment is 72.0 % forested with two main stand types: (1) stands dominated by subalpine fir (Abies lasiocarpa) and Engelmann spruce (Picea engelmannii), and (2) stands dominated by lodgepole pine (Pinus contorta) and western larch (Larix occidentalis).Clear-cuts in early stages of regeneration (mean tree heights up to 5.5 m) and two bedrock outcrops comprise 27.5 and 0.5 % of the catchment,  Soils throughout the UEC catchment are dominated by sands and silts with abundant coarse fragments.They developed primarily in deep (in excess of 2-8 m) morainal tills with some isolated areas of colluvium (BC Geological Survey, 2012).Except at a limited number of isolated ridge-top outcrops, bedrock is not observed throughout most of the catchment, including along most road cuts, some of which exceed a depth of 8 m.Based on visual observations, the majority of vegetation roots reside within the upper 30 cm of soil, with a lower root density between 30 and 50 cm below the surface.Although soil physical properties (particularly soil texture, coarse fragment content, and porosity) vary considerably across the catchment, vertical variations are, for the most part, gradual with little distinct soil layering.Generally, soils vary gradually from low density and high permeability at the surface to higher density and lower permeability at depth; however, some sites show negligible change (both  (Smith, 2011).Large soil macropores or cracks were generally not observed, likely due to the limited amount of vegetation roots below 30 cm depth, the absence or limited abundance of burrowing animals (e.g.small mammals, earthworms) in the catchment, and the low clay content of the soils.The only exceptions were at the heads of ephemeral streams where large macropores were observed that had likely developed via subsurface erosional processes.Observations were made while handdigging over 400 soil pits (0.4 to 1.7 m depths) in the catchment, including at least 100 for measuring hydraulic conductivity profiles, and along several kilometers of new forestry road (2 to 8 m cut bank depths) constructed at the end of the study.

Sampling design
Fifty hillslope monitoring sites were established (33 in October 2005 and 17 in July 2006) at stratified random locations throughout the UEC catchment (Fig. 1).Stratified random sampling was used to minimize the potential for investigator bias in site selection, and to ensure that statistical inferences could be reliably extrapolated to the entire study catchment.The sample size was selected to maximize statistical power while ensuring the infrastructure could be maintained by one person, particularly during installation and snow sampling.Strata were defined based on elevation (50 % of the sites were established at locations above and 50 % below the mean catchment elevation), insolation (50 % of the sites at locations greater than and 50 % less than the mean annual potential solar radiation within the catchment), forest cover (25 % of the sites in clear-cuts or regenerating stands and 75 % in forested areas), and hydrogeomorphic position (20% of the sites in each of the following classes: riparian, concave-wet, concave-dry, convex-wet, and convex-dry).For the hydrogeomorphic classes, riparian was defined as being located within 10 m of the catchment or sub-basin main-stem channels.Concave versus convex was defined as positive and negative values, respectively, of the Laplacian operator computed from a 3 × 3 neighborhood of cells surrounding each cell of interest in a 25 m resolution digital elevation model (DEM) obtained from the BC Ministry of Forests, Land and Natural Resource Operations.Wet versus dry was defined as a topographic wetness index greater than and less than the catchment mean, respectively.Elevation, insolation, and forest cover were selected for catchment stratification because they strongly influence snow depth, timing and intensity of melt, amount of evapotranspiration, and soil wetness.
Hydrogeomorphic position was selected because of its association with subsurface runoff processes via flow path convergence and divergence.The term hillslope hollow is used hereafter to refer to areas of pronounced surface concavity.The DEM analyses were conducted using Rivertools 3.0 (Rivix LLC, 2012).Potential solar radiation was modeled using Solar Analyst in ArcGIS 9.3.1 (ESRI, 2012).

Site installations
At each site, one groundwater well was manually installed (access limitations prevented bringing heavy equipment to the site) to the greatest depth possible, with the maximum depths limited by the abundance of large cobbles and boulders.Wells were selected rather than piezometers in order to capture the timing of initial groundwater response (i.e.water table measured within the well) and subsequent water table dynamics rather than capturing only hydraulic head at a specific depth in the soil.For the initial installation, PVC wells with a 3.8 cm inner diameter were installed in soil pits of approximately 30 cm diameter that were dug by hand using augers, shovels, picks, and pry-bars.They were screened by cutting narrow slits in the sidewalls with 2 to 3 cm spacing and wrapping the pipes with geotextile to prevent the potential influx of sediments.The soil pits were carefully back-filled with native soil, ensuring the original layering and avoiding compacting, which resulted in porosities that were similar to the original soils.At sites with limited or no observed groundwater responses, up to two additional attempts were made to increase the depths of the wells to improve the chances of observing groundwater responses, including installation of stainless steel drive-point wells using a sledgehammer (15 of 50 sites).The steel wells were screened with narrow slits at 2 mm spacing.Driving the steel wells into the soil inhibited wrapping them with geotextile; however, sedimentation never became an issue.Soil disturbance around the drive-point wells was minimal.After the final installations were complete, the minimum, mean, and maximum well depths were 0.50, 1.09, and 1.64 m, respectively.The mean depth of the unresponsive wells was 1.12 m (SD = 0.34 m) compared to 1.07 m (SD = 0.29 m) for the responsive wells.The wells were screened from the well bottom up to an average depth below the soil surface of 8 cm.Water table depth (relative to the soil surface) at each well was recorded every 30 min using Odyssey capacitance water level recorders (0.8 mm resolution) (Dataflow Systems Pty Limited, 2012).
A PVC pipe was installed within the steel wells to insulate the Odyssey sensors from potential interference.At 6 of the 50 hillslope sites, additional automated instruments were installed in October 2006.ECH 2 O sensors and Decagon loggers (Decagon Devices Inc., 2012) were used to record volumetric soil water content and air temperature, as well as water input (i.e.snowmelt/rainfall) depth from snowmelt lysimeters.These 6 sites are referred to as lysimeter sites, whereas all 50 sites that monitor hillslope runoff processes (including the 6 lysimeter sites, Fig. 1) are referred to as hillslope sites.
Precipitation, air temperature, incoming short-wave radiation, relative humidity, wind speed, and snow depth were obtained from two automated climate stations located in regenerating clear-cuts: the UC climate station (1780 m elevation, 750 m south of the UEC catchment) and the Lower Cotton (LC) climate station (1390 m elevation, 1500 m southwest of the UEC catchment).Additional details of the data collection infrastructure are provided in Table 1 and in Smith (2011).R statistical software (R Development Core Team, 2010) and Aquarius (Aquatic Informatics Inc., 2012) were used for processing the automated time series data.

Soil characterization
Soil samples (approximately 500 g dry mass) from the 45-50 cm depth at each of the 50 hillslope sites were analyzed to quantify the porosity, texture, and fractions of coarse fragment and organic matter.After burning the soil samples at 500 • C for 4 hours to remove organic matter, grain size distributions were determined using wet sieving for particles larger than 0.05 mm and a sedigraph for smaller particles.Due to the relatively small size of the soil samples, the fractions of coarse fragments are not representative of particles larger than approximately 1 cm in diameter.At two sites, multiple samples were gathered from a range of depths up to 1.1 m to assess vertical variations in soil properties.Vertical variation in soil properties at each site was also noted from field-based observations (including manual soil texture tests) made during installation of groundwater wells and soil moisture instruments.
Field-saturated hydraulic conductivity (K s ) was measured using a Guelph Permeameter at approximately 0.25 m (49 sites), 0.50 m (47 sites), 0.75 m (39 sites), 1.00 m (10 sites), 1.25 m (3 sites), and 1.50 m (1 site) soil depths.The resulting K s values were averages of vertical and horizontal conductivities since depth and width dimensions of the water-filled portion of the borehole were approximately equal; however, hydraulic conductivity is likely to be relatively isotropic throughout the catchment since the soils are dominated by sand, silt, and gravel with minimal amounts of clay and since the soils are not stratified (Mitchell, 1993).For each measurement, a 7 cm diameter hole was augered to the desired soil depth and two K s tests were conducted -one with 5 cm of hydraulic head and the other with 10 cm of head.K s was calculated for each test using methods described in the Guelph Permeameter operating instructions (Soil Moisture Equipment Corp, 1991) and the arithmetic mean (the geometic mean resulted in a negligible difference) of both tests was used as the final K s value.

Field measurement campaigns
Between October 2005 and September 2008, field measurement campaigns were conducted in early February and early April, every 2 to 4 weeks from April to early June, and once each in early summer, late summer, and mid-fall.SWE and snow depth were measured manually during site visits throughout winter and spring.Soil saturation was measured year-round by manually inserting an AquaPro capacitance probe (AquaPro Sensors, 2012) to the desired depth in an epoxy access tube that had been installed in the soil during the snow-free season.Measurements were made at 10 cm depth intervals to the maximum installation depth (40 to 90 cm depending on the size and abundance of coarse fragments, with a mean of 66 cm).PVC extension tubes were added to the epoxy access tubes for the winter period to facilitate soil saturation measurements below the snowpack.SWE was measured using a federal snow sampler.At each hillslope site, five snow samples were spaced at 4 m intervals on Hydrol.Earth Syst.Sci., 18, 1835-1856, 2014 www.hydrol-earth-syst-sci.net/18/1835/2014/ a contour across the hillslope centered at 5 m upslope from the groundwater well.

DEM development
A 5 m resolution DEM was developed using photo interpretation methods applied to 1 : 15 000 scale aerial photos (lidar was cost prohibitive).The raw photo interpretation points were supplemented with GPS points (Trimble ProXT GPS and Ranger data logger) gathered over a minimum area of 100 m × 100 m centered over each hillslope site.The final DEM was interpolated to a 5 m resolution using triangulation and was smoothed using a Gaussian filter.R statistical software (R Development Core Team, 2010) was used for merging and filtering the raw point data.SAGA GIS (SAGA User Group Association, 2012) was used for grid-averaging the point data and for interpolating and smoothing the DEM.

Site parameters
Table 2 provides a list of site parameters that were measured or calculated to characterize each hillslope site and were used in statistical analyses to investigate the influences of site physiography on groundwater response.Tree height, diameter, and basal area were obtained from forest cover survey data (7.5 m radius plot around each groundwater well).For calculating the portion of the upslope drainage area that is logged, all areas in an early stage of regeneration were considered logged.Clear-sky fraction was estimated from hemispherical photos taken horizontally over each groundwater well.Slope gradient was measured in the field.Rivertools 3.0 was used for calculating curvature from the 5 m DEM.ArcGIS 9.3.1 (ESRI) was used for delineating and calculating the upslope drainage area (using a D8 grid) for each site and for manually determining the elevation rise and fall to the upslope ridge and downslope channel (along the flowpath), respectively.To account for the effects of site aspect and hillslope shading on energy inputs, potential solar radiation (i.e.direct and diffuse, excluding forest cover effects) was modeled for each day of the year for each site using Ar-cGIS 9.3.1.Radiation was then averaged for three seasons of the year to account for seasonal variation: the snow accumulation season (November through March), the snowmelt season (April and May), and the snow-free season (June through October).The soil parameters were obtained using methods discussed in Sect.2.2.3.

Ordinal logistic regression
Groundwater responses were temporally discontinuous at most sites and detectible groundwater responses were never recorded at 13 of the 50 hillslope sites.Due to this data censoring, statistical analysis methods such as ordinary regression could not be applied without removing the unresponsive sites, which would have led to a substantial loss of spatial information.As a result, the analyses were based on an ordered classification (Sect.2.3.3) of the groundwater response data and ordinal logistic regression (OLR) was applied (Sect.2.3.3) to predict the probability of a site meeting or exceeding each ordered class using the selected site parameters (Table 2) as predictors.OLR is an extension of binary logistic regression (LR).LR forms a linear regression between the natural logarithm of the odds ratio (O) for a response variable and one or more predictor variables: where a is a constant, b is a vector of slope coefficients, x is a vector of predictor variables, and is an error term.The odds ratio for the response is the probability of being in one group divided by the probability of being in the other group, where p is the probability of a response being in a given or higher level category.An extensive review of OLR can be found in McCullagh (1980).An example of OLR applied in a watershed context can be found in McCleary et al. (2011).

Groundwater response classification
The time series data were limited to the period from 1 November 2007, to 20 September 2008, as this period was common to both the groundwater and the streamflow data sets.Although the period of record does not cover a full year, the missing data from late September and October are of little concern since the autumn period is relatively dry in the UEC catchment.The period of record was separated into eight hydrologically distinct streamflow periods to investigate potential links between spatial patterns of groundwater response and seasonal changes in catchment hydrologic conditions.Each period was intended to represent a distinct phase of water input (related to variability in meteorological conditions) and resulting runoff response (Fig. 2): (1) a fall transition period when the catchment experienced limited soil water recharge following the previous summer drought, (2) a winter lowflow period when the catchment experienced minimal water input, (3) an early-melt transition period when the catchment began experiencing active snowmelt input that generated a small streamflow response, (4) a rising limb period when snowmelt in the catchment generated a rapid rise in the streamflow response, (5) a peak-flow period when the streamflow response reached its maximum, (6) a falling limb period when the streamflow decreased quickly and the last of the remaining snow covered areas within the catchment were melting rapidly, (7) a post-melt transition period when no snow remained in the catchment and streamflow continued decreasing at a moderate rate, and (8) a summer lowflow period when streamflow responded to occasional intense rainstorm events.Three types of response variables were defined for the groundwater data: (1) occurrence, in which a well was assigned a value of 1 if a groundwater response was observed during the November 2007 to September 2008 period of record and 0 if it was unresponsive (no formation of a saturated layer within the observed soil profile); (2) duration, computed by determining the fractional portion of time that a water table was recorded in a well at any depth, and then reducing this interval measure to ordered classes for each time period; and (3) timing, in which the date/time (in decimal days since 1 January) of first response and maximum response were classified into ordered classes.Duration classes were defined for the eight streamflow periods individually (Fig. 2), then for the melt period after aggregating streamflow periods 3 through 6, and again for the annual period after aggregating all eight streamflow periods.OLR analyses were applied to all three levels of aggregation.winter, the date/time of the start of the first distinct rise in the water table level during the spring melt was used as the timing of first response.For other sites, the first response was defined as the date/time that a water table was first recorded in the well.Table 3 provides a summary of the classes for each period/timing.Observations of both transient and persistent (i.e.long duration) groundwater were treated as one population for the analyses because of sample size limitations.OLR requires the number of cases within each response class to exceed the number of predictor terms in the model, which restricted the number of classes that could be defined to two or three.As much as possible, natural breaks in the distribution of the response data were used to define class thresholds, but it was necessary to adjust the thresholds slightly for each streamflow period to meet the necessary sample size for each class based on the distribution of durations in the response data.Applying an ordered classification to the duration data also led to a loss of information since duration does not account for variation in groundwater response intensity (e.g.maximum groundwater level, rate of rise or fall); however, it was considered more important to maximize the spatial coverage than to capture more details of the groundwater dynamics.
For calculating duration of response data, we considered using a consistent reference soil depth (rather than calculating the portion of time that a water table was recorded in the well at any depth); however, doing so would have been problematic because of the large ranges in well depths and maximum water table levels measured in the wells.Both factors were governed by the soil conditions, with the latter influenced also by the local runoff processes.For any particular reference soil depth, the sample of positive responses would have been quite small for two reasons: (1) exclusion of sites from the sample due to well depths being above the reference soil depth, and (2) assigning zero values to sites due to their maximum water table levels being below the reference soil depth.Given that statistical approaches utilizing ordered classes (which was necessary due to data censoring) require much larger sample sizes than approaches utilizing continuous response variables (e.g.ordinary regression), the statistical power for investigating the data set would have been overly limiting.In response to this concern, we decided to incorporate well depth into the statistical models to test for its influence, as well as numerous soil parameters, including hydraulic conductivity profile data.Given that the average well depth for unresponsive sites was greater than that for responsive sites (Sect.2.2.2), it is apparent that other factors (e.g.soils, topography, vegetation) determined the occurrence, duration, and timing of groundwater responses, not well depth.

Applying OLR
We considered splitting the response data into independent calibration and test groups, but doing so was not feasible because of the small sample size.Logistic regression requires much larger sample sizes than ordinary regression.As a result, a formal test of the predictive power of the OLR models was not possible.However, the lack of an independent test of predictive accuracy does not diminish the value of conducting a diagnostic analysis focused on testing hypotheses about spatial controls on groundwater response.
In total, 36 site parameters were considered as candidate predictors in the OLR models (  parameters were included because we chose to analyze stationary (or relatively stationary) physiographic variables.Snowpack condition would be captured implicitly through its influence on groundwater response.Moreover, including parameters describing snowpack condition would weaken the apparent effects of the other physiographic parameters in the models due to covariance.It would be more appropriate to treat snowpack as a response variable while specifically investigating physiographic influences on snowpack processes, as studied by Jost et al. (2007) for the CCEW.
For each candidate predictor variable, a histogram was examined to assess whether the data appeared to be normally distributed, and an appropriate transformation was applied if the distribution exhibited a lack of normality (Table 2).Although OLR does not assume any particular sampling distribution, it is known to perform better when the predictor variables are normally distributed (Tabachnick and Fidell, 2007), and this finding was true for the current study.All predictor variables were standardized to minimize computational issues related to multicollinearity and to enhance interpretation of the model coefficients.To reduce the number of potential predictor variables for each model, the classified response data were fit to each potential predictor variable separately and the strongest predictor variable from each parameter group (e.g.forest cover group, soil constituent group; Table 2) was selected to enter the model first.Individual variables and variable interactions were then added, removed, or replaced to achieve a final model for each streamflow period.Any effects that were not physically meaningful or possible were removed from the models.Since the groundwater wells were installed to varying depths, well depth was included as a potential predictor variable to assess whether or not the well installation depths biased the observed responses.The effect of each predictor variable was calculated as the exponential of the following product: the coefficient for the variable multiplied by the range in the data between the 25th and 75th percentile values for the variable.Interaction terms were ignored for calculating main effects.For calculating interaction effects, the corresponding interaction variables were held at their respective minimum or maximum values.
The Wald test statistic (Engle, 1983), Akaike's information criterion (AIC) (Akaike, 1987), and the Bayesian information criterion (BIC) (Schwarz, 1978) were all used for variable selection, with an emphasis on BIC, as it led to the most parsimonious models.OLR assumes that the coefficients that describe the relationship between the predictor variables and the response category are the same for each response category in a model, which is called the proportional odds assumption or the parallel regression assumption.If this assumption were not true, different sets of predictor coefficients would be required to describe the relationship with each response category.OLR also assumes that the relationship between the predictor variables and the natural logarithm of the odds ratio of the response (Eq. 1) is linear.Plots of partial residuals were used to confirm that the proportional odds assumption was met, to check that predictors behaved linearly, and to check for potential outliers.Lastly, a bootstrap resampling validation procedure generating Somers' D rank correlation (Somers, 1962) and R 2 index statistics was applied to assess predictive performance.For LR, R 2 is referred to as an index because the residuals in LR are always the difference between a binary value (0 or 1) and the calculated probability and, therefore, R 2 is not strictly the same in LR as in ordinary linear regression.R statistical software (R Development Core Team, 2010) was used for all statistical analyses.The response variable for the top row is the occurrence of a groundwater response.The response variable for the middle and lower rows is the duration of groundwater response.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values, except interaction variables, which were held at their respective minimum or maximum values.The timing of each period is indicated Sect.2.3.3 and Fig. 2.

Space-time patterns of water inputs to the catchment
Figure 3a shows the distribution of mean annual potential solar radiation (hereafter referred to as insolation) to highlight the catchment areas that likely experience higher rates of evapotranspiration, sublimation, and snowmelt due to greater energy inputs (Ellis et al., 2011;Jost et al., 2007Jost et al., , 2012;;Toews and Gluns, 1986).These factors influence the dynamics of soil wetness, snowpack retreat, and water input.Figure 3b shows the patterns of snowline retreat for the winters of 2007 and 2008 through snow cover mapping.The groundwater response analyses utilized 2008 data only; however, the frequency of snow cover surveys was low in 2008, so snow cover data for both years are presented.It was observed during field investigations over 3 years that the snowline retreat patterns were generally consistent from year to year, but with differences in timing.It is acknowledged that the snow cover mapping shows the distribution of potential melt (i.e.binary information) only, and not the amount or rate of melt.However, other studies have shown the important connection between topography and snowmelt processes (Ellis et al., 2011;Jost et al., 2007Jost et al., , 2012;;Toews and Gluns, 1986), including results from the CCEW.The general spatial pattern of spring snowline retreat and, thus, the spatial shifting of the lower extent of snowmelt input to the soil progresses as follows: (1) south-facing, low-elevation clear-cut areas; to (2) south-facing, middleelevation forest and clear-cut areas; to (3) south-facing, high-elevation forest areas; and north-facing, low-and middle-elevation forest and clear-cut areas; to (4) northfacing, high-elevation forest areas.Locations along the valley bottoms receive low amounts of solar energy, regardless of the slope aspect, due to hillslope shading.

Groundwater response occurrence
A water table was never observed in the well (i.e.no groundwater response) at 13 of the 50 hillslope sites.The probability of a groundwater response occurring increased with increasing upslope drainage area or increasing melt period insolation, but decreased with increasing slope gradient or Table 4. Ordinal logistic regression models for the occurrence of response over the annual period; the duration of response over the annual period, melt period, and periods 1 through 8; and the timing of first and maximum responses over the annual period.For modeling, the measured response values were replaced with the class values of 0, 1, or 2, as provided in Table 3 4 for the OLR models and Fig. 4 for the partial probabilities).Based on the main effects for the predictor variables, deep-soil K s was most important in determining the probability of occurrence, followed by upslope drainage area, insolation, and slope gradient, in order of importance (refer to Fig. 5 for the strengths of the effects).However, accounting for variable interactions, slope gradient had a highly negative effect on the probability of occurrence among low-insolation sites and a weakly positive effect among high-insolation sites (Figs. 4 and 5).Insolation had a strongly positive effect on the probability of occurrence among high slope gradient sites and a slightly weaker negative effect among low slope gradient sites.Inspection of the groundwater response data shows that unresponsive sites (indicated by the "x" symbol in Fig. S1a of the Supplement; represented as class 0 in Table 3) tended to be at middle or upper hillslope locations among areas with planar surface curvature and on ridges, which is consistent with the model results.Manual comparison between the distributions of unresponsive sites and deep-soil K s (i.e.75 cm K s ) shows consistency with the model result that deep-soil K s was important in determining the probability of occurrence.These results verify that the statistical models are consistent with the hydrologic processes.The absence of unresponsive sites within the south-facing clear-cut area suggests that forest cover might also be an important variable, but none of the forest cover parameters was significant in the model, possibly due to their being overwhelmed by the influence of insolation or due to statistical limitations associated with the sample size.

Groundwater response duration
The OLR models for groundwater response duration included two or more of the following variables: upslope drainage area, slope gradient, deep-soil K s , maximum tree diameter, and insolation (Table 4; Figs. 4 and 6).Several models also incorporated interactions between slope gradient and insolation.Other variables tested (Table 2) were either not significant, explained less variance than the selected variables, or had effects that were not physically meaningful.Well depth was not significant in any model.

Annual and melt periods
Over the annual period, the probability of a persistent response increased with increasing upslope area or insolation, but decreased with increasing slope gradient, maximum tree diameter, or deep-soil K s (Fig. 4).In order of importance, the main effects were strongest for maximum tree diameter and upslope area, somewhat weaker for deep-soil K s and slope gradient, and much weaker for insolation (Fig. 5).Interaction effects showed that slope gradient had a stronger negative effect among sites with low insolation, and a weakly positive effect among sites with high insolation (Figs. 4 and 5).Insolation had a much stronger positive effect among sites with high slope gradient and a strongly negative effect among sites with low slope gradient.Inspection of the groundwater response data shows that sites with persistent annual responses (indicated by symbol size in Fig. S1a of the Supplement) tended to be located near streams or in hillslope hollows, particularly those with low slope gradients.Moreover, 36 % of the clear-cut sites experienced a 0.75-1.0response (i.e. a water table was measured within the well 75 to 100 % of the time; largest symbol size in Fig. S1a  for the variables in the logistic regression models predicting groundwater response duration.For each effect class, the strength of the effect was (1) 1-5, (2) 5-10, (3) 10-50, (4) 50-100, and (5) > 100, as outlined in Table 5.The effects were calculated as the exponential of the following product: the coefficient for the predictor variable multiplied by the range in the data between the 25th and 75th percentile values for the variable.Interaction terms were ignored for calculating main effects.For calculating interaction effects, the corresponding interaction variables were held at their respective minimum or maximum values.The streamflow periods are (1) fall transition, (2) winter low-flow, (3) early-melt transition, (4) rising limb, (5) peak-flow, (6) falling limb, (7) post-melt transition, and (8) summer low-flow (see Sect. 2.3.3 and Fig. 2 for details).
whereas only 11 % of the forested sites experienced a 0.75-1.0response.The variables in the OLR model for the melt period, their interactions, and their signs (i.e.positive or negative effects) were the same as for the annual model, but the order of importance varied (Figs. 4 and 5).In particular, slope gradient and deep-soil K s had stronger effects than both maximum tree diameter and upslope area.Moreover, the overall strengths of the main and interaction effects were stronger for slope gradient and weaker for maximum tree diameter and upslope area.Compared to the annual data, some sites distant from the stream network and not in well-defined hillslope hollows experienced more persistent responses (Figs.S1a  and b of the Supplement).

Individual streamflow periods
For the individual streamflow periods, upslope drainage area and slope gradient were first or second in importance in all periods, except periods 3 (early-melt transition) and 4 (rising limb) during the early phases of the spring freshet (Fig. 5).The main effects of upslope area and slope gradient were positive and negative during all periods, respectively, and were at their lowest importance in periods 3 and 4, respectively (Figs. 5 and 6).
Deep-soil K s was third in importance in periods 1 through 3 and 8, and was the most important variable in period 4 (rising limb), but was not significant in periods 5 through 7 once the catchment reached its wettest condition.Maximum tree diameter was fourth in importance in periods 1, 2, and 4, and was the most important variable in period 3 (early-melt transition).The effects of deep-soil K s and maximum tree diameter were always negative.Insolation was significant only in periods 3 through 5 when snowmelt was widespread throughout the catchment.The main effect of insolation was weakly positive in period 3, moderately positive in period 4, and moderately negative in period 5, showing the spatial shifting of active water inputs from highinsolation areas in early spring to low-insolation areas later in the spring.Moreover, the effect of insolation was strongly negative among low slope gradient sites and strongly positive among high slope gradient sites in periods 3 through 5 (Figs. 5 and 6).During the same periods, slope gradient had a strongly negative effect among low-insolation sites, but the effect among high-insolation sites increased from weakly positive to moderately positive.Interestingly, sites with low slope gradient and high insolation, and sites with high slope gradient (regardless of insolation) had drier shallow soils at the start of the spring melt compared to sites with low slope gradient and low insolation (Fig. 7).
By examining the relative positions of the main effects along the respective x axes in Fig. 6, one can observe that the relations shift to higher or lower values of the predictor variables sequentially between periods.To investigate this variation in more detail, the 0.1 and 0.5 probabilities of a persistent response are plotted for each predictor variable (except insolation since the direction of its effect, i.e. positive versus negative, changes between period 4 and period 5) and period in Fig. 8. Compared to the 0.1 probability, the 0.5 probability had a higher upslope area and lower slope gradient, deep-soil K s , and maximum tree diameter in any given period.These relative positions are consistent with the positive main effect of upslope area versus the negative main effects of slope gradient, deep-soil K s , and maximum tree diameter on the duration of groundwater response (Figs. 5 and 6).The minimum upslope area required to generate persistent groundwater responses reached a maximum in period 2 or 3 and a minimum in period 4. In period 4, 0.1 and 0.5 probabilities were associated with upslope areas of 0.009 and 0.11 ha, respectively.In contrast, the corresponding upslope areas were For calculating the box plot statistics, the hillslope sites were stratified into four groups combining low or high slope gradient (less than or greater than the mean gradient) with low or high solar radiation (less than the lower quartile or greater than the upper quartile of mean potential solar radiation throughout the snow-free season).
1.6 and 140 ha in period 3, respectively.Similarly, slope gradient reached a maximum in period 4 and a minimum in periods 2 or 3.In period 4, 0.1 and 0.5 probabilities were associated with slope gradients of 68 and 32 %, respectively, whereas the corresponding slope gradients were 20 and 9 % in period 3.This spatial expansion of groundwater response to higher slope gradients and smaller upslope drainage areas (i.e. to planar hillslopes and ridges) followed by its contraction can be observed in maps showing the 0.75-1.0duration sites (i.e.groundwater response during 75 to 100 % of the period) (Fig. S2 of the Supplement).The most widespread distribution occurred in period 5.Moreover, from period 3 through period 5, moderate-and high-elevation sites were sequentially added to the sites experiencing persistent responses while a small number of low-elevation sites stopped contributing.
The limiting effect of deep-soil K s on groundwater response (Fig. 8) was greatest in period 1 (i.e.only sites with very low values of K s were likely to experience persistent groundwater responses) when the catchment was relatively dry and streamflow was low following the summer drought, and became less limiting (i.e. higher K s value) in subsequent periods.Deep-soil K s was the least limiting (i.e. even sites with high values of K s were likely to experience persistent groundwater responses) in period 4 during widespread snowmelt and was insignificant in the models for periods 5 through 7.In period 1, 0.1 and 0.5 probabilities were associated with deep-soil K s values of 9 × 10 −9 m s −1 and 5 × 10 −10 m s −1 , whereas the corresponding deep-soil K s values were 4 × 10 −5 m s −1 and 4 × 10 −6 m s −1 in period 4, respectively.
Maximum tree diameter was most limiting (i.e.only sites with small or no trees were likely to experience persistent responses; Fig. 8) in period 1 when the catchment was relatively dry and streamflow was low.Maximum tree diameter was least limiting (i.e.sites with small through large diameter trees were likely to experience persistent responses) in period 4 during widespread snowmelt and was insignificant in the models in periods 5 through 8.In period 1, forested sites had probabilities of experiencing persistent responses that were less than 0.1, whereas probabilities of 0.1 and 0.5 were associated with maximum tree diameters of 72 and 33 cm in period 4, respectively.The limiting effects of forest cover on melt-related groundwater response can be observed for period 3 when sites within the low-elevation south-facing clear-cut area experienced persistent responses, but sites in adjacent forested areas of similar or lower elevation and similar insolation did not respond or responded minimally (Fig. S2c of the Supplement).For periods 4 and 5, sites with persistent responses were distributed throughout both forested and clear-cut areas.

Groundwater response timing
The first response occurred earlier with increasing upslope drainage area or with more upslope logging, but later with increasing elevation, deep-soil K s , or silt fraction (refer to Fig. 9 for the change in probability along each variable gradient and Table 5 for the strengths of the variable effects).The main effects were strongest for upslope drainage area followed by elevation, silt fraction, deep-soil K s , and upslope logging (Table 5).No interaction effects were significant.Inspection of the spatial distribution of first response (indicated by symbol size in Fig. 1c of the Supplement) shows that sites with an early first response tended to be located near Table 5.Effect class, strength of the effect, direction of the effect (positive or negative), and effect rank (indicated numerically to the right of the variable name) for the variables in the logistic regression models predicting groundwater response timing (first and maximum responses over the annual period).The effects were calculated as the exponential of the following product: the coefficient for the predictor variable multiplied by the range in the data between the 25th and 75th percentile values for the variable.The maximum response occurred earlier with increasing insolation and clear-sky fraction, but later with increasing silt fraction (Fig. 9).Clear-sky fraction had the strongest effect, with much weaker effects from insolation and silt fraction, in order of importance (Table 5).No interaction effects were significant.Inspection of the maximum response data (Fig. S1d of the Supplement) shows that sites with an early maximum response tended to be located in the low-elevation south-facing clear-cut.

Controls on groundwater response occurrence
Deep-soil K s (measured at 75 cm depth) was the most important variable for predicting whether or not a site experienced any detectable groundwater response (i.e.occurrence of response) within the range of soil depths that were monitored (i.e.groundwater well depths ranged between 0.50 and 1.64 m).Among high slope gradient sites, a groundwater response is more likely at high-insolation sites than at lowinsolation sites, despite having negligible differences in shallow soil wetness at the start of the spring melt.These findings, coupled with peak snowmelt intensities likely being greater on high-insolation sites than on low-insolation sites (Ellis et al., 2011;Jost et al., 2007Jost et al., , 2012;;Toews and Gluns, 1986) suggest that water input intensity is an important influence on the occurrence of groundwater response.These results are consistent with the percolation-excess runoff generation mechanism described by Redding andDevito (2008, 2010), who found that water input intensity was a first-order control on the occurrence and amount of lateral flux in glacial till soils due to its influence on vertical versus lateral flux partitioning.Similar results were found for the UEC catchment, as described by Smith (2011).In the statistical models for this study, deep-soil K s likely accounts for variation in the depth of the percolation-limiting layer.
Upslope drainage area and slope gradient were also important in determining whether or not a site experiences a detectable groundwater response.This finding suggests that some sites with large upslope areas will experience a shallow groundwater response regardless of the K s of the surficial soils due to high rates of flow accumulation.Jencso et al. (2009) also found a positive relation between the occurrence of groundwater response and upslope drainage area for sites in a snowmelt-dominated catchment in Montana.However, notwithstanding these findings, the current study suggests that the underlying geology and the various physiographic influences on snowmelt dynamics might be equally or more important than topographic convergence in determining the spatial distribution of responsive sites.

Controls on groundwater response duration
During the early phases of the spring freshet, while the catchment is relatively dry (except along riparian corridors) and snow covered, increasing energy inputs begin to generate localized melt in low-elevation, high-insolation areas, similar to patterns observed by DeBeer and Pomeroy (2010).Under these conditions, the OLR models suggest that the vertical controls of localized energy and mass inputs, and vertical versus lateral flux partitioning expressed by maximum tree diameter and deep-soil K s , respectively, dominate the patterns of groundwater response duration.Once snowmelt expands throughout the catchment and hydrologic response spans large gradients in insolation, forest cover, and topography, the lateral controls of hydraulic gradient and flow path convergence expressed by slope gradient and upslope contributing area, respectively, begin to dominate.Lateral controls dominate throughout the peak-flow period and throughout the summer, autumn, and winter low-flow periods while the catchment drains.These findings are supported by those of Jencso et al. (2009), Jencso and McGlynn (2011), and Kuras et al. (2008) for other snowmelt-dominated montane catchments, and by those of Szeftel (2010) for the CCEW.These findings also corroborate the applicability of the generally accepted relations between soil wetness and various topographic indices (Quinn et al., 1995;Beven and Kirkby, 1979), as well as the importance of topographic position as a controlling factor in runoff-generation dynamics (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000), except during the early-melt and rising limb periods of the spring freshet.It highlights a temporal shift in the relative importance of typology and topography in controlling the persistence of shallow runoff response (Buttle, 2006). .Partial probability of groundwater response for predictor variables in the logistic regression models.The response variables for the top and bottom rows are the timing of maximum response and the timing of first response, respectively.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values.At any value along the variable gradient, the timing of response (i.e.early, middle, or late in the melt period) with the highest probability is the most likely outcome.
The contrast in the dominance of vertical versus lateral controls is also highlighted by the positive influence of insolation on the duration of groundwater response at high slope gradient sites during the early-melt, rising limb, and peak-flow periods compared to the negative influence at low slope gradient sites during the same periods.The same patterns exist for the probability of groundwater response occurrence and the duration of groundwater response over the annual and melt periods.Among low slope gradient areas, low-insolation sites are generally wetter at the start of the spring melt than high-insolation sites (Fig. 7), likely due to lower rates of pre-melt evapotranspiration (particularly before snowpack development), making low-insolation sites more responsive to water inputs.This finding is corroborated by field observations that north-facing areas are generally wetter throughout the snow-free season and have a higher density of streams compared to south-facing areas.In contrast, among high slope gradient areas, rapid soil drainage likely limits the influence of evapotranspiration on antecedent soil wetness at the start of the spring melt, resulting in greater influences from snowmelt rates.In essence, controls on antecedent wetness and flow path convergence overwhelm controls on water input intensity and vertical versus lateral flux partitioning among low slope gradient sites, but not among high slope gradient sites.
Forest cover was an important variable in several models.We believe this finding was related to the typically negative influence of forest cover on snowpack accumulation (i.e. total water input) and melt intensity, and positive influence on evapotranspiration (i.e.antecedent wetness).In particular, maximum tree diameter consistently explained more variance than other forest cover parameters.We believe this finding was related to the largest diameter trees in a stand (and their relatively wide crowns) likely having disproportionate influences on snowpack shading and evapotranspiration.

Controls on groundwater response timing
The timing of first groundwater response is controlled by parameters influencing the upslope hydrologic conditions and lateral redistribution (e.g.upslope drainage area, percentage of upslope logging), the local soil hydraulics (e.g.deepsoil K s , silt fraction), and the localized energy inputs and/or snowpack depth (e.g. as influenced by elevation).The importance of upslope drainage area shows that when the upslope area is large, even small amounts of melt can generate a response, likely due to the influence of lateral redistribution of soil water on antecedent wetness.The importance of upslope logging is likely related to the influence of forest cover removal on snowmelt timing and on evapotranspiration and, thus, soil wetness prior to winter snowpack development.A high value of deep-soil K s means that more storage capacity must be satisfied (i.e.due to greater depth to the percolation-limiting layer) before the first groundwater response can occur.
The timing of maximum groundwater response appears to be controlled primarily by parameters influencing localized energy inputs (e.g.clear-sky fraction, insolation) with less control by parameters influencing local soil hydraulics (e.g.only silt fraction is important) and negligible control related to lateral redistribution.Neither upslope drainage area nor slope gradient are important controls on the timing of maximum response, suggesting that vertical controls dominate lateral controls, which is consistent with the maximum response being controlled primarily by surface processes.The importance of clear-sky fraction and insolation in controlling the timing of maximum response illustrates the importance of localized energy and mass flux dynamics on the differential timing, intensity, and quantity of snowmelt and their subsequent influences on the timing of response.

Implications for runoff response and catchment modeling
Sites with high values of deep-soil K s that do not generate a shallow groundwater response should experience deep percolation and likely generate runoff via slow response pathways, resulting in continual drainage throughout the recession and low-flow periods.These findings are supported by two points: (1) precipitation exceeds actual evapotranspiration in the UEC catchment and, therefore, all sites must experience runoff (ignoring the influences of wind, which is negligible under a forest canopy, on snowpack distribution and, thus, the water balance); and (2) rapid response pathways within the deep subsoil, such as deep-soil cracks in clay or bedrock (Montgomery et al., 1997(Montgomery et al., , 2002;;Trompvan Meerveld et al., 2007), are likely limited in abundance and extent since the soils are typically 2-8 m or more in depth with only small amounts of clay and minimal bedrock.Thus, the spatial distribution of surficial soil K s is an important control on the distribution of sites that generate shallow rapid runoff versus sites that generate deep slow runoff.Consistent with our findings, Kuras et al. (2008) found that runoff patterns during low-flow periods were not explained well by surface topography and suggested that deep, inefficient flow paths dominate runoff generation during these periods.Jencso and McGlynn (2011) found that increasing proportions of permeable geology underlying relatively wet landscapes were correlated with decreased streamflow yield during wet periods and increased yield during dry periods.Moreover, sites with large upslope drainage areas and high values of soil K s that also experience a shallow groundwater response would be capable of transmitting water to the stream network at a high rate and, thus, would be critical in the connectivity of shallow runoff response areas to streams.The amount of incident solar radiation controls the distribution of shallow rapid runoff response areas during the early-melt period and, to a lesser extent, during the rising limb and peak-flow periods.However, its influence weakens as active snowmelt zones shift into more shaded (via topography and/or forest cover) locations.Peak streamflow occurs when snowmelt and shallow runoff are being generated throughout most areas of the catchment, including locations with low insolation, with mature forest cover, with high slope gradients, with convex topography (i.e.ridges), and with high elevations (Figs. 8 and S2 of the Supplement).This expansion and contraction of shallow runoff areas likely results in variable runoff efficiency at the catchment scale.It is broadly analogous to the patterns described by conceptual rainfall runoff models (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000), and is consistent with the findings of Jencso et al. (2009), Kendall et al. (1999), andKuras et al. (2008).
Topography-based indices (Beven and Kirkby, 1979;Quinn et al., 1995) and the various rainfall runoff conceptual models (Dunne and Black, 1970a, b;Freeze, 1972;Hewlett andHibbert, 1963, 1967;Sidle et al., 2000) should be reliable predictors of shallow runoff areas in snowmeltdominated montane catchments during peak-flow, recessionflow, and low-flow periods since upslope area and slope gradient were dominant in the statistical models during these periods (Jencso and McGlynn, 2011;Jencso et al., 2009;Kuras et al., 2008;Szeftel, 2010).However, they would likely be poor predictors of shallow runoff areas during early phases of the spring freshet without adequately addressing the spacetime variability of water input intensity (i.e.controls on snowpack conditions and surface energy fluxes).Kuras et al. (2008) found that differential snowmelt timing between clear-cuts and forested areas was responsible for generating different streamflow peaks.Moreover, for glacial till catchments with spatially variable soil K s profiles or for catchments with varying soil depths, catchment models must address the spatial distribution of controls on vertical versus lateral flux partitioning in the soil coupled with the distribution of controls on the rate of percolation to adequately explain runoff patterns during all periods, including differentiation of shallow rapid runoff response areas from areas dominated by deep flow paths that contribute primarily to sustaining lowflows (Redding andDevito, 2008, 2010).
This latter point has important implications for the spatial discretization of a catchment for hydrologic modeling.Models that represent lateral variability in a fully distributed manner, such as DHSVM, are capable of representing the complex and temporally variable controls on the spatial variability of groundwater response and runoff generation.However, semi-distributed models are desired for many applications that require computational efficiency, particularly for generating uncertainty estimates using GLUE (generalized likelihood uncertainty estimation) or similar procedures (e.g.Freer et al., 2010).Many semi-distributed models, such as HBV and its variants (e.g.Hamilton et al., 2000;Lindstrom et al., 1997), use a grouped response unit (GRU) approach that does not explicitly represent the effects of lateral flow between the GRUs and, thus, are not capable of representing the catchment dynamics revealed by our observations.Consequently, we suggest that an appropriate scheme might utilize separate discretizations for above-and belowsurface processes.For example, a traditional GRU approach could be used to model processes such as canopy interception and snowpack development, but a more complex approach would be required for subsurface hydrology.Further research should focus on developing and testing such schemes and comparing their predictive abilities to traditional GRU and fully distributed approaches.

Conclusions
The spatial controls on the occurrence, duration, and timing of shallow groundwater response in glacial till montane catchments that are snowmelt-dominated are complex and vary not only between seasons, but also intra-seasonally.The K s of the soil at 75 cm depth was found to be a firstorder control on the distribution of sites that generate shallow groundwater response versus sites that experience only deep percolation.Moreover, the study findings suggest that sites with highly permeable surface soils, large upslope contributing areas, and low slope gradients would be important in connecting shallow rapid runoff response areas to streams.Upslope contributing area and slope gradient are first-order controls on the duration of shallow groundwater response during peak-flow, recession-flow, and low-flow periods.Shallow runoff areas expand and contract throughout these periods coincident with catchment wetting and drying, and follow the general spatial patterns of topographic convergence.However, controls on the differential timing, intensity, and quantity of snowmelt and controls on vertical versus lateral flux partitioning in the soil overwhelm the influence of topographic convergence during early spring freshet periods.In essence, a temporal shift occurs in the relative importance of typology and topography in controlling the spatial patterns of shallow runoff response.These findings suggest that various topographic indices and topography-based rainfall runoff models would likely be poor predictors of runoff patterns during early phases of the spring freshet without adequately addressing controls on water input intensity and on vertical versus lateral flux partitioning in the soil and rate of percolation.However, their performance would likely improve as the freshet and post-freshet periods proceed.
The Supplement related to this article is available online at doi:10.5194/hess-18-1835-2014-supplement.

Figure 1 .
Figure 1.Location and study sites of the Upper Elk Creek catchment.

Figure 3 .
Figure 3. Mean annual potential solar radiation (a) and snow cover extent during the spring melt periods of 2007 and 2008 (b).

Figure 4 .
Figure 4. Partial probability of groundwater response for predictor variables in the logistic regression models.Streamflow periods are indicated on the far right side of the rows.The response variable for the top row is the occurrence of a groundwater response.The response variable for the middle and lower rows is the duration of groundwater response.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values, except interaction variables, which were held at their respective minimum or maximum values.The timing of each period is indicated Sect.2.3.3 and Fig.2.

Figure 5 .
Figure5.Effect class, direction of the effect (positive or negative), and effect rank (indicated numerically to the left of the data point) for the variables in the logistic regression models predicting groundwater response duration.For each effect class, the strength of the effect was (1) 1-5, (2) 5-10, (3) 10-50, (4) 50-100, and (5) > 100, as outlined in Table5.The effects were calculated as the exponential of the following product: the coefficient for the predictor variable multiplied by the range in the data between the 25th and 75th percentile values for the variable.Interaction terms were ignored for calculating main effects.For calculating interaction effects, the corresponding interaction variables were held at their respective minimum or maximum values.The streamflow periods are (1) fall transition, (2) winter low-flow, (3) early-melt transition, (4) rising limb, (5) peak-flow, (6) falling limb, (7) post-melt transition, and (8) summer low-flow (see Sect. 2.3.3 and Fig.2for details).

Figure 6 .Figure 7 .
Figure6.Partial probability of groundwater response for predictor variables in the logistic regression models.For slope gradient, mean of upslope and downslope gradients versus downslope gradient only is indicated in the panels.Streamflow periods are indicated on the far right side of the rows.Response variables for all rows are the duration of groundwater response.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values, except interaction variables, which were held at their respective minimum or maximum values.The timing of each period is indicatedSect.2.3.3 and Fig. 2.

Figure 8 .
Figure 8. Variation in the predicted values of (a) upslope drainage area, (b) slope gradient (downslope gradient in black, mean slope gradient in gray), (c) deep-soil K s , and (d) maximum tree diameter for the 10 and 50 % partial probabilities of groundwater response duration.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values.The streamflow periods are (1) fall transition, (2) winter low-flow, (3) early-melt transition, (4) rising limb, (5) peak-flow, (6) falling limb, (7) post-melt transition, and (8) summer low-flow (see Sect. 2.3.3 and Fig. 2 for details).

+
Silt fraction (3) streams and in hillslope hollows, and in the low-elevation south-facing clear-cut.
Figure 9. Partial probability of groundwater response for predictor variables in the logistic regression models.The response variables for the top and bottom rows are the timing of maximum response and the timing of first response, respectively.For calculating partial probabilities, all other predictor variables were held at their respective mean (geometric mean for K s ) values.At any value along the variable gradient, the timing of response (i.e.early, middle, or late in the melt period) with the highest probability is the most likely outcome.

Table 1 .
Data collection infrastructure at the streamflow, hillslope, and lysimeter sites.

Table 2 .
Physiographic parameters and corresponding transformations applied for the logistic regression analyses.Parameter symbols are defined only for parameters retained in the final models (Table4).

Table 3 .
Range of groundwater responses for the response variables modeled using ordinal logistic regression (i.e. for modeling, the measured response values were replaced with the class values of 0, 1, or 2).The response variables are the occurrence of response over the annual period; the duration of response over the annual period, melt period, and periods 1 through 8; and the timing of first and maximum responses over the annual period.Units for duration are the fractional portion of the period.Units for timing are decimal day of the year.