A hydropedological approach to quantifying groundwater recharge in various glacial settings of the mid-continental USA

: Data from six monitoring stations were combined with a soil-water dynamics model (HYDRUS 1D) to achieve physically-based estimates of shallow water-table recharge in representative hydropedological settings of the glaciated midwestern U.S.A. Calibration involved inverse modeling that yielded optimized hydraulic parameters. Root mean square errors for modeled versus measured soil moisture contents were generally within 3% for all soil layers at the six study sites. The optimized models also accurately simulated recharge values that corresponded to observed water-table ﬂ uctuations. Optimized parameter values were consistent with estimates from a pedotransfer function, lab analyses, and ﬁ eld experiments. Forward modeling indicated that shallow water-table recharge in mid-continent glacial settings is approximately 35% of precipitation, but interannual and monthly variability is signi ﬁ cant. Soil parent materials and horizon characteristics in ﬂ uence recharge primarily through their control on K s with clay-rich till parent materials producingvalues aslow as16% and coarse-grainedoutwash parent materials producingvaluesashigh as58% ofprecipitation.During the three-year study period, distinct seasonality of rechargewas observed with most recharge occurringin the winter (seasonal meanof all sites was 66% of precipitation) and lesser but interannually stable amounts in the spring (44%), summer (13%), and autumn (16%). This research underscores the importance of incorporating pedological information into models of soil-water dynamics and groundwater recharge. © 2015 The Authors. Hydrological Processes published by John Wiley & Sons Ltd.


INTRODUCTION
Groundwater resource planning requires quantitative understanding of water fluxes in the vadose zone, and there is a need for critical-zone research that integrates pedological characteristics with hydrologic models to better our understanding of water movement through the soil profile to the water table (Lin et al., 2008). Recent studies have identified the role of hydropedological units in the hydrologic function of high-relief glaciated catchments (e.g. Soulsby et al., 2000;Tetzlaff et al., 2014), but those findings may not be directly applicable in low-relief areas that were subjected to continental glaciation. In low-relief landscapes with humid continental climates such as the Great Lakes region (GLR), USA, the water table is commonly less than 5 m below ground surface, so percolating soil water readily enters the groundwater-flow system as recharge. Gannon et al. (2014) noted the uniqueness of groundwater response in different areas of an alpine glaciated catchment and identified the need for experimental data on vadose zone and groundwater dynamics in the different hydropedological units. In this research, we distinguished hydropedological units in continental glacial settings on the basis of topographic setting, soil parent material, and degree of weathering as indicated by profile development. Because of the diversity of hydropedological settings, existing published recharge estimates for regional aquifers in the GLR may not represent the true range of values that exist. The variability may be important because perched water tables are common and these are critical to the function of surface water features like intermittent streams and groundwater-fed wetlands. Additional research on the magnitude and variability of recharge is all the more important when one considers the uncertainty associated with many estimation techniques (Scanlon et al., 2002;Gates et al., 2014).
Recharge estimation approaches in glaciated regions are diverse and often quite indirect. For example, Walton (1965) combined hydrograph separation and water balance approaches to estimate recharge to aquifers in Illinois, USA, with results ranging between 6 and 24% of annual precipitation. Arnold and Allen (1996) extended the hydrograph approach of Walton (1965) and determined recharge rates of 10 to 28% of annual precipitation for three additional watersheds in Illinois. Nolan et al. (2007) combined chloride-tracer and Darcy pedotransfer methods to achieve recharge values between 0.3 and 63% of annual P in the GLR. Delin et al. (2007) estimated recharge rates through glacial deposits in Minnesota, USA, using various forms of the water table fluctuation method (e.g. Healy and Cook 2002) and determined that recharge ranged between 16 and 26% of annual precipitation. Fragalà and Parkin (2010) calculated recharge rates between 11 and 14% of annual P using chloride mass balance and nitrate tracer methods in both glacial till and alluvium and further evaluated recharge in different slope positions in the rainfall-dominated Upper Eden Valley, Cumbria, UK. Note that none of these studies used direct physically based methods to arrive at groundwater recharge estimates and all integrate recharge over the year without considering the seasonality of recharge mechanisms. In contrast, Cuthbert et al. (2010) used a more direct approach based on field observations of groundwater levels and soil tension, along with numerical modelling to derive recharge estimates between 6 and 17% of annual precipitation for fractured lodgement till in Shropshire, UK. Cuthbert et al. (2010) were also able to infer recharge seasonal dynamics using the relatively fine temporal resolution of their data set.
Any direct physically based approach to determine recharge through unsaturated materials must account for porous-media variability due to soil horizon development. Soil horizon variability is significant in the GLR, where Alfisols are the most common soil type; the Bt horizons of Alfisols have translocated clays that produce lower permeability relative to overlying A horizons and most underlying parent materials. Additional alterations may exist in clay-rich till deposits because of desiccation and freeze-thaw cycles among other processes (McKay et al., 1993;Brockman and Szabo 2000;Cuthbert et al., 2010). Biological alterations may also affect vadose zone water fluxes via worm burrows and root traces, and varying land management practices such as tillage can impact the pedological character of A horizons. These alterations, along with the thickness, texture, and organic content of A horizons, determine the distribution of flow-active pores that govern surface water allocations between runoff and infiltration (Lipiec et al., 2006). Franzmeier (1991) and Franzmeier et al. (1991) quantified soil characteristics such as porosity, bulk density, water content, and saturated hydraulic conductivity in soils in the GLR (Indiana) using 'lithomorphic' class (a combination of soil texture, development, and parent material) as an indexing method. Parameter estimates and generalized soil-classification schemes presented in such studies are essential first steps for incorporating pedological heterogeneity into direct methods for determining recharge.
Subsurface variability greatly influences recharge in humid settings where diffuse recharge is the dominant recharge mechanism (Scanlon et al., 2002). Shallow groundwater is recharged when percolation is in excess of soil moisture deficits and evapotranspiration demand (ET; Lerner et al., 1990;de Vries and Simmers 2002;Sophocleous 2004). However, Rushton (1997) cautions that surface water and vadose zone analyses can yield only potential recharge estimates, and emphasizes the importance of groundwaterlevel data as a necessary element in the determination of actual recharge. Collectively, the diverse findings on recharge to date indicate the following guidelines for direct modelling of water table recharge in the GLR: (1) include data representing distinct hydropedological units, (2) deploy monitoring stations designed to quantify soil moisture storage changes in conjunction with reference evapotranspiration and water table fluctuations, and (3) incorporate continuous data into numerical models of unsaturated groundwater flow. These are all important elements of what Lin (2003) described as the emerging field of hydropedology.
Research herein integrates soil laboratory characterization with measurements of vadose zone water content and pertinent micrometeorological data to support a physicallybased analysis of shallow groundwater recharge and watertable dynamics in six hydropedological settings. This data-driven approach provides a more direct basis for identifying recharge mechanisms that are temporally dynamic and spatially diverse-essential insights under emerging climate scenarios typified by intra-annual variability in precipitation and evapotranspiration. HYDRUS 1D (Šimůnek et al., 2005) is used as a parameter estimation tool to derive hydraulic parameters through an inverse method but more importantly as a modelling tool to improve our understanding of vadose zone dynamics (including water table recharge). To constrain and validate hydraulic parameters used in forward models, parameter selection approaches are also compared: pedotransfer values based on soil bulk density and texture, field and lab values based on in situ or lab testing, and best-fit values based on inverse modelling. The combined field-based and modelling experiments are practical for characterizing hydropedological units in the GLR and as an initial step for future catchment-scale investigations.

Monitoring sites and approach
Data from six monitoring locations representing hydrogeologic settings and soil parent materials that are common within the GLR have been collected over the past 3.5 years (2011-2015; Table I, Figure 1). Soils at each study site were described in the field using standard methods described in Birkeland (1999). Cores were also collected 1595 HYDROPEDOLOGICAL APPROACH TO QUANTIFYING RECHARGE IN GLACIAL SETTINGS from trench walls for laboratory analyses, and electronic sensors were installed to continuously monitor moisture content, soil water tension, and temperature. The trenches were carefully backfilled after installing the monitoring devices to restore the solum (weathered soil profile) back to its original order. Vegetation over the excavated area was also restored to the surrounding land use (Table I).
Cylindrical core samples (5 cm diameter × 10 cm length) collected at depth increments of 0.3 m were used to determine bulk density (ρ s , Table II) using standard methods described by Grossman and Reinsch (2002). Volumetric water content (VWC or θ) was determined from cores using the approach of Topp and Ferre (2002). Water retention relations for a draining scenario were determined for one subsoil core sample from each site using a modified version of the pressure cell method presented by Dane and Hopmans (2002); the modified method entailed rinsing the bottom of the porous plates with deaired water between pressure intervals instead of using a burette to remove air from below the porous plate. Soil samples collected at 0.3-m depth increments were also subjected to particle-size analysis using the simplified clay fraction procedure described by Gee and Or (2002). The resulting soil texture data (Table II) have been normalized by removing an equivalent weight ratio representing the gravel size fraction (>2 mm) for each depth so that the US Department of Agriculture (USDA) soil-classification scheme can be applied for pedotransfer functions. Particle densities were determined at each site using a gas pycnometer, and these data were used to determine porosity (ϕ , Table II) through a method outlined by Flint and Flint (2002). All of the monitoring sites are located in low-relief terrain, including terraces, slightly dissected hill crests, and a floodplain ( Table I). Three of the sites (AL, SGT, and OT) are underlain by soils having parent materials from a proglacial depositional environment. The site located at the edge of a floodplain is composed of siltdominated overbank alluvial (AL) deposits in the upper 2.4 m that are underlain by glaciolacustrine clay. Site OT is an outwash terrace (~10 m above the current river level) and has the coarsest sediments with the highest sand contents. Site SGT is a terrace composed of supraglacial till likely deposited by debris flows; the parent materials contain a slightly higher percentage of clay (23-29%) compared with that of typical outwash. Soils from the remaining sites (GM, EM1, and EM2) are glacial sediments associated with moraines. Site GM (ground moraine) is located in a till plain with parent materials that were deposited subglacially. Sites EM1 and EM2 (end moraine) are on relatively high-relief landforms that were deposited at the terminus of a glacier, and both sites are within 0.5 km of existing wetlands, which were more widespread in these settings prior to modern agricultural drainage practices. The till parent materials have higher clay contents (Table II) and have bulk densities and porosities that are also generally higher than parent materials at the other sites. Site EM1 is unique relative to the other sites because it contains a 1.2-m-thick layer of loess at the surface, whereas loess occurrence tends to be focused in the upper 50 cm at other sites. The predominant soil horizons identified during trenching at all sites (Table II) include the following: Ap (organic-rich horizon disturbed by plowing), Bt (strong structure with significant accumulation of translocated clay), Bw (clay-rich but with low translocated clay content and weak structure), Cu (unweathered parent material), and Cox (oxidized parent material).
Continuous monitoring was initiated at five of the six sites in late summer, 2011, whereas monitoring at site GM began in early October 2012. Dataloggers were programmed to record from micrometeorological and vadose zone instruments at an hourly time interval. Soil tension was measured using matric potential blocks installed at 0.15-, 0.3-, 0.6-, and 1.2-m depths (heatdissipation probes were used to measure tension at the GM site). Temperature sensors and water content reflectometers with 30-cm-long rods were installed at 0.3-m intervals to continuously measure thermal and moisture gradients, respectively. Dielectric permittivity from water content reflectometers were related to volumetric water content using the equation of Topp et al. (1980)  except at GM, EM1, and EM2, which were subjected to site-specific calibration due to high clay contents and bulk densities. The following micrometeorological instruments were also installed and programmed to record hourly data at each site: pyranometer, anemometer, tipping bucket rain gauge, and relative humidity sensor. Winter precipitation data were taken from local weather stations that are part of the National Centers for Environmental Information using heated rain gauges during periods when snowfall was recorded. The reference stations also include snow depth observations; these were used to develop simple snowmelt estimates where precipitation was summed over a period when snow pack observations indicated continuous snow cover and re-distributed evenly over those days when a snow pack reduction was reported (presumed to be melt days). The winter snowmelt precipitation results were used as precipitation inputs for hydrological modelling during the winter season. Shallow monitoring wells (Table III) were installed near the micrometeorological and vadose zone instruments at four sites between October 2012 and April 2014. Wells were constructed of 5.1 cm diameter PVC casing with screen intervals that were 30-60 cm in length and surrounded by gravel filter packs. The wells at GM and OT existed prior to this study and are in the same hydropedological settings as the vadose zone stations (same parent material, terrain, and weathering regime). All continuous groundwater-level data were measured and recorded using pressure transducers with internal dataloggers.
The micrometeorological data (precipitation, insolation, ambient air temperature, relative humidity, and wind speed) were used to estimate surface energy and water fluxes. Reference evapotranspiration (ET 0 ) using daily time step and a short grass surface were calculated following the method described by Allen et al. (1998).

Description of the HYDRUS 1D numerical models
Hydropedological settings for this study were chosen on plains, terraces, or ridge tops to be consistent with a 1-D flow approximation through the unsaturated zone. Soil horizons at the study sites are nearly horizontal features as indicated in trenches. Ground slopes at the six sites based on LiDAR derived digital elevation models varied between 0.1°and 2.4°(mean 1.1°). Numerical experiments using a 2-D cross-sectional model with dipping model layers and representative parameters revealed vertical to horizontal flow ratios that exceeded 17 : 1 for the two sites having slopes greater than 1.5°s upporting the vertical flow assumption. Therefore, it was concluded that employing HYDRUS 1D to simulate soil water flow and evapotranspiration from the soil column was appropriate for all six GLR sites.
Many other researchers have used HYDRUS 1D to examine both recharge (Jiménez-Martínez et al., 2009;Zhu and Sun 2009;Lu et al., 2010;Leterme et al., 2012) as well as solute and contaminant transport (Cheviron and Coquet 2009;Rozemeijer et al., 2009;Rieckh et al., 2014). HYDRUS 1D version 4.16 (Šimůnek et al., 2005) solves the Richards equation for vertical unsaturated groundwater flow, and the models in this study apply the van Genuchten (1980) and Mualem (1976) equations for soil hydraulic parameters. Model layers were assigned to each site based on the analysis of soil horizon designations in Table II. Initial soil-layer hydraulic parameters were assigned using: (1) lab ϕ as saturated moisture content (θ s ), (2) field-saturated hydraulic conductivity (K s ) for upper horizons measured using a Guelph permeameter, and (3) the Rosetta pedotransfer function (Schaap et al., 2001) for residual-moisture content (θ r ) and van Genuchten fitting parameters n and α.
In HYDRUS, plant water uptake is simulated using the methods of Feddes et al. (1978), where uptake is partitioned into each layer according to the depth-specific root density based on the vegetation. The depth-specific root density was determined based on the maximum root depth noted in trench descriptions for the turfgrass (EM1 and GM) and prairie/conservation sites (AL, EM2, OT) with transient root depths fluctuating with leaf area index (LAI; the maximum root depth coincides with the maximum annual LAI). Maximum corn (0.93 m) and soybean (0.8 m) root depths for the corn/soybean crop rotation site (SGT) were based on rooting depths in till reported by Dwyer et al. (1988) with transient conditions fluctuating with LAI. The water stress response function used in Feddes et al. (1978) was estimated based on surfacelayer-vegetation type and was taken from the database contained in HYDRUS 1D (Šimůnek et al., 2005). Water stress response parameters for turfgrass were used for the EM1 and GM sites (lawns); pasture values were used for the AL, EM2, and OT sites (prairie/conservation); whereas bean and vegetative-stage corn parameters were used for the SGT site (corn/soybean rotation). LAI values for turfgrass settings (EM1 and GM) were based on a maximum LAI of 2.25 and seasonal growth patterns described for the GLR in Gelernter and Stowell (2005). LAIs for corn and soybeans were interpolated based on data obtained from the Bondville Ameriflux site (Wilson and Meyers 2007), whereas prairie/conservation LAIs were interpolated based on data from the Fermi Ameriflux site (Gomez-Casanovas et al., 2012).

HYDRUS 1D boundary conditions
Upper and lower model boundary conditions were established using methods described by Lu et al. (2010). Surface infiltration is simulated when precipitation exceeds evapotranspiration, as long as the pressure head at the soil surface exceeds a critical minimum value calculated based on daily relative humidity and air temperature (not to exceed À15 000 cm in this study). Once the critical value is reached during periods of high evaporation, actual evaporation is decreased as a timevariable flux because the soil is too dry to lose water to the atmosphere at the potential rate. When the surface becomes saturated during wet periods, surface runoff occurs when precipitation exceeds infiltration.
Treatment of the lower boundary condition is important for this study of shallow groundwater systems because we define groundwater recharge as infiltrating water that percolates below the root zone and increases storage in the underlying saturated zone. Consistent with Lu et al. (2010), the lower boundary conditions were defined as either variable head (with head established by the water table position) or free drainage (a unit vertical gradient at the base of the domain) based on the water table depth relative to the soil profile at each site (discussed further in the results section). Downward flux at the base of the HYDRUS model domain is commonly linked to groundwater recharge (e.g. Lu et al., 2010); however, because measured groundwater rose to within the lower portion of our model domains at five sites, water table flux was derived from the variable head models. This was accomplished by taking the first occurrence (starting deeper in the profile) of downward flux between 1.8 m and the rooting zone for each model (1.8 m is the lowest model output node used for comparison with measured data). If downward flux occurred only within the rooting zone (which occurred when rain coincided with a period of high evapotranspiration), then water-table flux and recharge were assumed to be zero. Hereinafter, modelled groundwater recharge values are determined using the water table flux estimation approach and referred to as 'recharge'.
Initial estimates of soil hydraulic parameters were refined by optimization via an inverse method. Saturated moisture content (θ s ) was restricted to within 10% of labderived ϕ, but other parameters were fit without such constraints. Minimization of the difference between simulated and observed VWC for each model layer was employed as the objective criterion. Parameters were adjusted in HYDRUS between iterations using the Levenberg-Marquardt procedure (Marquardt 1963) until a minimum was reached between measured and simulated VWC. Following the model optimization phase, forward model runs were initiated using simple data assimilation. This was accomplished by the direct insertion method, using observed VWC data as initial soil moisture conditions. Sites AL, OT, and SGT were initiated at the beginning of each water year (1 October) while the moraine sites (EM1 and EM2) were started at the onset of high water table conditions (see data in results section) to achieve the best representation of the water table for each setting.

Monitoring data
The 2012 water year (10/2011-9/2012; WY 2012) included drought conditions during the summer (sites averaged 14 cm of precipitation over 3 months compared to 31 cm normally received during the same period), while near-normal amounts of precipitation were experienced throughout the 2013 water year (10/2012-9/2013; WY 2013), and slightly above normal precipitation occurred during the 2014 water year (10/2013-9/2014; WY 2014) (Figure 2). Reference evapotranspiration equaled or exceeded precipitation at three of five sites during WY 2012, and the measured soil moisture storage, calculated by integrating volumetric water content over the measured profile, decreased at these three sites by an average of 5.5 cm during the 12-month period. In contrast, precipitation exceeded ET 0 at all six sites during WY 2013 by 20% and was within 12% of the 30-year normal precipitation on average. Therefore, monitoring data for WY 2013 were used for inverse simulations and parameter estimation. Soil moisture and groundwater-level data for all six sites ( Figure 3)  zone hydrologic regimes at each site. The site with alluvium over glaciolacustrine sediment (AL) includes a thick sequence of silt loam and soil moisture down to 1.2 m decreased significantly during the WY 2012 drought indicating moderately well-drained soil conditions. The water table response at this site and the site underlain by supraglacial till soil parent material (SGT) indicates that groundwater levels respond abruptly to infiltrating precipitation and then return to base-level conditions. The outwash terrace site (OT) is underlain by a thick unsaturated zone relative to the other sites (groundwater levels are >5 m below the surface) and water-level rise occurs with long intermittent periods at near base-level conditions. The coarse subsoil at this site yields low water contents and very minor fluctuations throughout the year, which is characteristic of a welldrained site and soils having steep soil water retention curves. On the other hand, the moraine sites (GM, EM1, and EM2) have shallow water tables (<4 m) and prolonged periods of up to 6 months when water levels are high with less variability. Subsoil moisture ranges (indicated by the 1.2-and 1.8-m VWC data) are subtle, but do show that the profiles drain mostly in conjunction with the lowering water table following periods of profile saturation (1.2-m VWC decreases followed by 1.8-m VWC). The measured hydrologic regimes at each site provided an important context for developing the 1-D models used for groundwater recharge estimation. Because groundwater rose to within 2 m of the ground surface at five sites ( Figure 3), a variable head lower boundary condition was used for HYDRUS models, established by groundwater depth. Groundwater-level data were not available at some of the sites during WY 2012 and WY 2013 so estimated or 'synthetic' groundwater-level data were calculated using a multiple regression model that related deep (>1.2 m) soil moisture and soil tension data to measured groundwater levels (see dotted lines in Figure 3); the average R 2 for the regression models at five sites was 0.67. Groundwater data for the sixth site (OT) was available for the entire study period.

Results of numerical inversion and model validation
The best-fit parameter estimates (Table IV) show that K s values are typically highest for A horizons, except at SGT and OT where deep Bw and Cu horizons are coarser (these subsoil horizons also have relatively high α and n values as would be expected for well-drained sediments). High ρ s values for Cu and Cox horizons in basal till resulted in low laboratory ϕ values at the GM, EM1, and EM2 sites, and this is mirrored in the estimated values of θ s .
Percolation below the rooting zone is crucial to water table recharge, so accurate modelling of moisture dynamics in the lower layers of the soil deserves special attention. The best matches between modelled and measured moisture content at depth occurred at sites AL, OT, and SGT ( Figure 4). These sites are similar in that they have parent materials consisting of coarse-grained sediments deposited beyond an ice front (supraglacial settings). In contrast, sites with parent materials consisting of basal till (moraine settings) had lower correlations between predicted and observed moisture content (Figure 4), with site GM showing the poorest fit of all. Forward modelling was not conducted for the GM site because the timing of water movement through the lower model layer did not appear to be accurately depicted by the HYDRUS model. The poor model performance at this site is attributed to prominent fractures in the glacial till that likely result in dual porosity and/or permeability that may require a modelling approach that is not based on the Richards equation. However, the other five sites have low-layer root mean square error (RMSE) values within 2% coupled with good visual fits suggesting that groundwater recharge magnitudes and timing are reasonably accurate. Indeed, a forward model validation phase conducted for a 258-day period when continuously measured groundwater hydrograph data were available at the sites shows that hydrograph response compares well with modelled recharge (Figure 5), especially for supraglacial settings. This correlation was tested using the nonparametric Spearman rank correlation method, showing good correlation for the supraglacial sites (AL, OT, SGT; |x| = 0.6), whereas the moraine sites (EM1, EM2; |x| = 0.2) were less well coupled. Because the groundwater response represents a broader area than the local scale response represented by the profile VWC data, it is expected that the groundwater correlation would not be entirely consistent. Nevertheless, the reasonable agreement of groundwater level versus modelled recharge provides an important additional system response for model validation. With respect to profile VWC data, the validation period confirmed the accuracy of soil moisture simulation. The RMSE results for all layers from the validation period (AL = 0.01, OT = 0.02, SGT = 0.01, EM1 = 0.03, and EM2 = 0.02) were nearly identical to the calibration period ( Figure 4).

Comparison of parameter estimation approaches
Estimated values of θ s and K s are of interest because they largely govern the role of soil drainage during recharge periods. To serve as a check on the suite of bestfit parameter values from numerical inversion, we compared values of θ s and K s to laboratory determinations from soil samples, field experiments, field-measured values for similar horizons published by Franzmeier (1991), and values based on Rosetta pedotransfer function (Table V). Because θ s is difficult to measure independently, porosity (ϕ) is often used as an estimation even though it is recognized that θ s is lower than ϕ in field settings (Smith et al., 2002). The statistical correlations between modelled θ s and lab porosity were, on average, within 3.8% of each other. Similarly, when comparing best-fit θ s values with Franzmeier's total ϕ measurements, the modelled values were 5.4% lower over the same range (ϕ ≈ 0.02), which is expected because most soils contain dead-end pores. The correlation between best-fit modelled θ s calculated from the Rosetta pedotransfer function showed modelled values greater by an average of 2.3%. Modelled θ s values compared well with pressure cell data, and the best-fit values were only 3.5% higher than the measured values. Overall, best-fit values of θ s based on numerical inversion of the Richards equation have a wider range compared to those using the alternative parameter estimation approaches. A possible explanation is that high-porosity A and B horizons may be more likely to develop macropores (e.g. worm burrows and fractures). Alternatively, in this study, C horizons having compacted basal till parent material may have higher bulk density and lower porosity when compared to estimates based on data from other sites. Thus, the values of θ s derived from numerical inversion may more accurately depict field conditions than the alternative methods based on laboratory measurements. With respect to K s , Guelph permeameter data compared poorly with inverse-derived K s (r = 0.32), whereas Rosetta pedotransfer estimates based on lab bulk density and texture inputs were in better agreement (r = 0.80) as were conductivity values from Franzmeier's (1991) field data (r = 0.81). K s values from numerical inversion were lower than values estimated using  Franzmeier (1991) for almost the entire range of modelled K s (except for the Bw horizon of site AL) but consistently within an order of magnitude.

Results of forward model simulations
Forward models were developed for each site using the optimized parameter estimates (Table IV) along with WY 2012, WY 2013, and WY 2014 forcing data. In order to test the general performance of forward models, soil moisture content was integrated for the measured profile (0-180 cm), and average RMSE was calculated between measured and modelled profile-integrated soil moisture for each site during the three water years. Consistent with Kendy et al. (2003), the RMSE values were calculated as a percentage of the average total water content at each site with the following results: AL = 4.3%, OT = 9.7%, SGT = 4.8%, EM1 = 9.7%, and EM2 = 7.4%. The values indicate a good agreement between measured and modelled soil moisture and are consistent with the model error reported by Kendy et al. (2003) who reported a 9.3% RMSE for their 1-D groundwater recharge model.
Monthly recharge for each of the five forwardmodelled sites is shown in Figure 6 along with total recharge for WY 2012, WY 2013, and WY 2014 (percentage of WY precipitation is shown in parentheses) demonstrating intra-annual and interannual variability in shallow groundwater recharge. With the exception of site SGT, all sites had the highest recharge ratio (recharge expressed as a percentage of precipitation) during the 2012 drought year, and this is attributed to large rain events in November and early December of 2011 that resulted in significant water table flux during a water year when total precipitation was low due to a subsequent summer drought. Recharge ratios were lowest at four of five sites during WY 2013 owing to reduced water table flux during the winter months when compared to the other water years in the study. Monthly recharge ratios were aggregated for all sites and lumped into seasons (Figure 7) to investigate the seasonality of shallow groundwater recharge because comparing site data between water years may not elucidate crucial temporal dynamics. Recharge during the winter months (Dec, Jan, and Feb) was highly variable, whereas the other seasons exhibited remarkably stable recharge ratios from one water year to the next. An analysis of seasonal recharge (water table flux) aggregated for sites having supraglacial soil parent materials (AL, OT, and SGT) and sites having parent materials deposited in moraine settings (EM1 and EM2) shows that a marked difference exists between the two groupings ( Figure 8). Recharge was higher for supraglacial settings in the winter and spring of all water years owing to high hydraulic conductivity and low evapotranspiration during these seasons. In contrast, shallow groundwater recharge was higher for the moraine settings during the fall (2011 and 2013) and summer (2012 and 2013) owing to higher water retention and flatter retention curves with smaller subsoil n values (the slope of the retention curve) and, therefore, more protracted recharge for these hydropedological units. Recharge was much greater for site EM1 than EM2 having distinct periods of increased water table flux ( Figure 6, WY 2012 Nov-Jan and WY 2014 Dec-May). Significant amounts of annual runoff (>5 cm) were simulated only at the end moraine sites (EM1 and EM2), with periods of runoff that occurred during all three water years at EM2 (11-23% of precipitation). Differences in the timing and rates of recharge were revealed when we compared shallow groundwater recharge between sites (and observed measured data, Figure 3), which supports the importance of properly distinguishing pedological model units (e.g. horizons) and assigning appropriate hydraulic parameters.
We examined our forward-modelled recharge results in relation to an independent calculation of recharge using groundwater hydrograph analysis (groundwater level data in Figure 5). Recharge values derived from the HYDRUS models for all sites (except for GM, where forwardmodelled recharge estimates were deemed unreliable) were compared with recharge values calculated from groundwater fluctuations using methods outlined in Healy and Cook (2002). Aquifer specific-yield values shown in Table III are based on tabled ranges from Johnson (1967) for the material textures. All recharge values from the water table fluctuation method deviated only between 1 and 11% from those derived from HYDRUS models. The best correspondence was with the supraglacial sites.

Parameter estimation
The variability of measured and modelled soil horizon parameters (Table IV) demonstrates the strong control that both parent material and pedological development have on hydraulic properties. Studies such as Boswell and Olyphant (2007) show that saturated hydraulic conductivity (K s ) is often the most sensitive parameter in unsaturated flow simulations, but these authors also noted that the quantification of K s in both field and laboratory settings is fraught with uncertainties. Some uncertainty was diminished in this study by putting greater emphasis on the definition and characterization of the primary soil horizons at each site. The best-fit values of K s were often an order of magnitude different within primary horizons of the weathering profile (Table IV). When comparing the  (Table IV) and lab data (Table II) indicated that soils having supraglacial parent materials generally showed an increase in porosity with depth (possibly from textural grading or accumulation of clay minerals with depth), whereas soils that developed within moraine settings demonstrated decreasing porosity with depth (likely owing to an increase in bulk density with depth). Inverse-derived K s (Table IV) Allred's (2000) values reported for fractured till with 0.0025-cm fracture aperture and spacing at 5 to 100 cm. More work is needed to determine accurate parameters for soils developed in till deposits and more attention may need to be placed on dual porosity/permeability modelling where highly fractured basal tills exist. Nevertheless, our findings with regard to parameter estimates reveal good correspondence among all of the methods for hydraulic parameter estimation.

Shallow groundwater recharge
Modelled recharge results (Figures 6 and 7) demonstrate the value of employing measurement-driven, deterministic models for recharge calculation. Although WY 2012 included a protracted period of summer drought, annual recharge totals ranged from 27 to 58% of precipitation because much of the recharge occurred in the preceding November and December when significant rainfall occurred. This highlights the temporally dynamic nature of vadose zone hydrology in much of the GLR. The seasonal variability of the modelled recharge ratios (Figure 7) indicates that shallow aquifers are most sensitive to potential contaminants during winter (66%) and spring (44%) months, and this has implications for water policy regarding the timing of agricultural nutrient applications. It is noteworthy that across all three water years of this study, the seasonal mean recharge ratio is remarkably stable for the spring (44%), summer (13%), and autumn (16%) seasons. Thus, despite the variability from site to site, we observed surprising stability in the seasonal mean response across all sites for all seasons except winter. These results are consistent with the findings of Cuthbert et al. (2010), who found summer recharge to be about half that of winter recharge in similar glacial sediments.
Although seasonal recharge ratios were fairly consistent from one year to the next when averaging all sites (Figure 7), there were distinct variations when comparing yearly ratios for end members such as the well-drained outwash terrace site (OT) and the poorly drained end moraine site with a shallow Cu horizon (EM2). The OT site had a much higher average recharge ratio for all 3 years (41% compared to 22%) and recharge that generally receded slowly following winter months. The EM2 site, however, experienced pulses of recharge during winter and spring that coincided with periods when the water table rose into the profile (Figure 3). Poorly drained moraine sites (GM, EM1, and EM2) typically had high residual water contents when compared to the coarser soils of supraglacial settings (AL, OT, SGT) suggesting that these hydropedological units store significant amounts of shallow groundwater in the landscape when water tables are high. Vast wetlands once covered these areas prior to modern farming and tile drainage, and although hydroperiods have been reduced or eliminated, the timing of fluxes to shallow water tables still plays a major role in controlling catchment hydrology for continentally glaciated landscapes.
In addition to parent material controls (e.g. subglacial vs supraglacial), the field data and recharge estimates demonstrate the control that pedologic characteristics (weathering) impose on the hydraulic properties (and drainage characteristics) of soil profiles in continentally glaciated terrains. Estimated recharge at the EM2 site, which has an unleached subglacial soil parent material (Cu horizon) within one metre of the ground surface, was much lower than the EM1 site, which has a leached and oxidized subglacial soil parent material (Cox horizon) with a larger K s value (Table IV). Keller et al. (1988) found similar trends when comparing two glacial tills having varying weathering characteristics in flat upland recharge area, which they attributed to variations in the groundwater flow regime that is governed by vertical flux rates, water table depth, and depth of oxidized till.
The 1-D model results and parameters presented herein reveal the importance of soil types and horizons in influencing shallow recharge and seasonal perching of groundwater, but the 1-D approach has limitations when considering watersheds as a whole. Future modelling should place more emphasis on understanding groundwater flow paths from glacial hydropedological settings to adjacent streams. Such modelling efforts are essential for understanding the fate and transport of non-point source pollution from agricultural land that is experiencing rapid increases in production demands and therefore nutrient and pesticide inputs. The dynamics of shallow groundwater systems in glaciated terrains are probably more complex than those of deeper groundwater aquifer systems. Future investigations 1606 S. NAYLOR ET AL.
would benefit from merging field data across multiple scales (e.g. Gates et al., 2014) with 3-D models to better understand the water balance and groundwater flow paths in glaciated landscapes.

CONCLUSIONS
This work provides examples of how pedological (parent material and solum) information when combined with instrumentation and a physically-based model can lead to a better understanding of soil moisture dynamics and water table recharge in a region having heterogeneous surficial geology and shallow water tables. Best-fit values of saturated hydraulic conductivity (based on a numerical inversion of the Richards equation using HYDRUS 1D) exhibited order of magnitude differences within soil profiles as well as between sites having differing soil parent materials. Modelled recharge values obtained from the optimized models were very different at five monitoring sites and ranged from 16 to 58% of P with the lowest (16-29% of P) at a site with unweathered clayrich basal till parent material and highest (R = 30-58% of P) where the parent material consists of glacial outwash. Modelled recharge through parent materials consisting of alluvium over lacustrine, thick loess over loamy till, and supraglacial till was intermediate with recharge ranging from 26 to 52% of P based on the three water years in this study, which experienced between 80 and 110% of 30-year normal P. The simulated recharge, when averaged across all sites and three water years, indicated that water table recharge was approximately 35% of P (median 34%; mean 35%). The high recharge rates provide a context for explaining how groundwater-fed wetlands and intermittent streams are perpetuated in this region of generally low relief. Distinct seasonality of recharge was observed, with most recharge being received in the winter and lesser but interannually stable amounts received in spring, summer, and fall. Further modelling refinements will ultimately be necessary to realize all the benefits of a fully hydropedological approach to studying water storage and transfer through the critical zone. More attention may need to be given to small scale details of soil morphology and structure (Lin 2003, p.3) especially in areas where prominent fractures in soils derived from basal till occur. In five of the six case studies presented here, good model fits (RMSE of VWC generally within 3%) were obtained without the need for dual porosity parameters. Finally, this study demonstrates the value of employing hydrometeorological stations in representative hydropedological settings as a basis for long-term monitoring and ultimately, more effectively representing how surface energy and moisture budgets influence the dynamics of groundwater and surface water interactions.