Vadose zone modeling to identify controls on groundwater recharge in an unconfined granular aquifer in a cold and humid environment with different meteorological data sources

Groundwater recharge (GR) is a complex process that is difficult to quantify. Increasing attention has been given to unsaturated zone modeling to estimate GR and better understand the processes controlling it. Continuous soil-moisture time series have been shown to provide valuable information in this regard. The objectives of this study were to (i) analyze the processes and factors controlling GR in an unconfined granular aquifer in a cold and humid environment and (ii) assess the uncertainties associated with the use of data from different sources. Soil moisture data monitored over three years at three experimental sites in southern Quebec (Canada) were used to calibrate the HYDRUS-1D model and to estimate ranges of possible GR in a region where groundwater is increasingly used as a source of fresh water. The simulations identified and quantified important factors responsible for the near-surface water balance that leads to GR. The resulting GR estimates from 2016 to 2018 showed marked differences between the three sites, with values ranging from 347 to 735 mm/y. Mean GR for the three sites was 517 mm/y for 2016–2018 and 455 mm/y for the previous 12-year period. GR was shown to depend on monthly variations in precipitation and on soil textural parameters in the root zone, both controlling soil-water retention and evapotranspiration. Monthly recharge patterns showed distinct preferential GR periods during the spring snowmelt (38–45% of precipitation) and in the fall (29% of precipitation). The use of different meteorological datasets was shown to influence the GR estimates.


Introduction
Sustainable groundwater management is needed to limit adverse anthropogenic impacts and ensure the long-term availability of the water resource. It begins with the assessment of current and future groundwater renewal rates via groundwater recharge (GR). This question has been extensively studied under semi-arid and arid conditions (i.e., Wang et al. 2009a;Assefa and Woodbury 2013;Turkeltaub et al. 2015;Baalousha et al. 2018), where low GR can significantly limit groundwater uses. Sustainable groundwater management has also been addressed in cold and humid climates where GR is often considered to be abundant and non-limiting (i.e., Dripps et al. 2007;Rivard et al. 2014;Boumaiza et al. 2020). However, even when water is abundant overall, increased temperatures due to climate change may trigger increased evapotranspiration rates during the summer when water needs are the greatest (e.g., agriculture and recreational activities) and this can have large impacts on groundwater renewal rates (Turkeltaub et al. 2015;Saha et al. 2017;Guerrero-Morales et al. 2020).
Under a wide variety of conditions, GR is still considered to be one of the most challenging water-balance components to quantify (Dripps and Bradbury 2007) because it is difficult to measure directly and it is strongly influenced by climate, soil heterogeneity and land use. A wide range of methods (Scanlon et al. 2002) has been shown to be useful in assessing the quantity of water actually replenishing an aquifer, including isotopic profiles (Barbecot et al. 2018;Mattei et al. 2020a, b), remote sensing techniques (Jackson 2002), chloride mass balance (Szilagyi et al. 2011), base flow separation (Hung Vu and Merkel 2019), unsaturated zone modeling (Hu et al. 2019;Mattei et al. 2020a), water balance models (Dubois et al. 2021) and water table fluctuations (Crosbie et al. 2005). The reliability of these methods depends on a variety of factors, including the availability of field data to describe site characteristics, such as soil type, vegetation and depth to the water table (Scanlon et al. 2002). The lack of continuous meteorological datasets or the sparse distribution of climate stations can also be limiting factors. To resolve this issue, several gridded climate products (i.e., from interpolation methods or simulated predictions coupled with data assimilation) have been developed recently, providing different meteorological variables at varying spatio-temporal resolutions (Mesinger et al. 2006;Hutchinson et al. 2009). Because different datasets may be more or less reliable depending on the study area, their impact on GR simulations needs to be evaluated (Choi et al. 2009;Langlois et al. 2009).
In recent years, unsaturated zone models have gained popularity because of their demonstrated capability to provide reasonable GR estimates (Jiménez-Martínez et al. 2009 (HYDRUS-1D); Assefa and Woodbury 2013 (HYDRUS-1D); Thoma et al. 2014 (HYDRUS-1D); Xie et al. 2018 (WAVES); Hu et al. 2019 (HYDRUS-1D)). These models can also be used to assess how different factors, such as hydraulic properties and meteorological variables, control groundwater renewal rates (Wang et al. 2014;Min et al. 2015;Turkeltaub et al. 2015;Wang et al. 2016;Hu et al. 2019). The work of Turgeon et al. (2018) has underlined the importance of including unsaturated zone processes to better represent watershed reactions to rain events and snowmelt.
The objective of this research was to analyze the processes and factors controlling GR in an unconfined granular aquifer in a cold and humid climate and to assess the uncertainties associated with the use of data from different sources. The HYDRUS-1D model was used to simulate GR rates in the Vaudreuil-Soulanges region in southern Quebec (Canada) (Fig. 1) where groundwater use is estimated to represent 30% of GR (Larocque et al. 2015). Soil moisture monitoring at three experimental sites within unconfined sandy aquifers over a three-year period (2016-2018) allowed for model calibration. Long-term GR was estimated the preceding 12 year period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) with the calibrated model and allowed for a detailed analysis of how different gridded climate products and soil parameters control GR.

Study area
The study area is located in the Vaudreuil-Soulanges region in southern Quebec (Canada), near the Ontario and USA borders (Fig. 1). This region is part of the Groundwater Recharge Research Infrastructure (Infrastructure de recherche sur la recharge des eaux souterraines -IRRES), which consists of GR monitoring stations distributed throughout the southern regions of Quebec (Canada). In the Vaudreuil-Soulanges region, it has been estimated that 20 × 10 6 m 3 of water is consumed each year, 54% from groundwater and 46% from surface water (Larocque et al. 2015). This percentage is significantly higher than the 20% of groundwater use estimated for the province of Quebec as a whole (MELCC 2019b).

Geology and hydrogeology
Silty clay Quaternary deposits inherited from the Champlain Sea are dominant in the lowland portion of the Vaudreuil-Soulanges region, which is mostly occupied by agriculture. This silty clay can reach a thickness of up to 30 m and significantly limit GR for the whole region. Using a spatially distributed surface water budget, Larocque et al. (2015) estimated that GR in the Vaudreuil-Soulanges region can range from 0 mm/y in the clay-covered lowlands to 440 mm/y in areas where substantial sand deposits are present. At the locations of the Saint-Telesphore esker and Saint-Lazare and Hudson hills (Fig. 1), well-drained and thick fluvioglacial deposits represent highly permeable granular aquifers and preferential recharge zones (Larocque et al. 2015). It is estimated that the Hudson and Saint-Lazare hills receive 41% of all regional recharge. Other fluvioglacial deposits of a smaller extent, including the unconfined part of the Saint-Telesphore esker and a recharge zone on Mount Rigaud, receive a mean annual recharge of 256 mm/y. It is known that aquifer recharge occurs in episodic events, mostly in the fall (October to December) and during the spring snowmelt (March to mid-April), but the exact distribution throughout the year has not been quantified yet.
Three sites were monitored in the current study: the Saint-Telesphore (STEL) site and two sites located on the Saint-Lazare hill (SLZA and SLZB) are located on sandy fluvioglacial deposits reaching a thickness of 40 m (STEL) and 86 m (SLZA and SLZB) above the regional bedrock aquifer (Fig. 1). The deposits are composed of thick and complex fluvioglacial deposits that were reworked at the surface and are now characterized as deltaic and littoral glaciomarine deposits. Depending on the energy regime and the progression of the ice margin, the fluvioglacial deposits on SLZ hill are either silty beds or have blocks within the top 1 m. These thick fluvioglacial deposits lay on discontinuous and variable beds of till and ancient Quaternary sand before reaching the sandstone bedrock. This highly variable spatial distribution of sediments leads to alternating confined and unconfined conditions for the deeper aquifer on SLZ hill. The aquifer at STEL site is locally unconfined.

Study sites
The two stations at the Saint-Lazare site were used to evaluate the effect of vegetation on soil water content and GR. SLZA is located in a small clearing area within the forest and covered with a sparse grass, while SLZB, located 15 m away, is located in a jack pine forest. The STEL is located close to a sand quarry and is dominated by dense prairie grass. The instruments were located 10 m from a forest. At the three locations, well-drained sandy soils were dominant at the surface. The thickness of the unsaturated zone varied between 4.8 and 7.4 m, and the average water table depth was 5.5 and 6.8 m below the surface at SLZ (site A and B; one piezometer between the two sites) and STEL, respectively.

Meteorological conditions
The 1981-2010 climate normal means for the three weather stations in the region (MELCC 2019a) report a mean annual air temperature of 6.1 °C with a maximum mean of 25.7 °C in July and a minimum mean of −15.2 °C in January. The mean total annual precipitation is 980 mm, of which 16% falls as snow between November and March. During the study period, the total annual precipitation was 984 mm/y (2016), 1206 mm/y (2017), and 829 mm/y (2018). The maximum daily mean temperature was 26 °C (2016 and 2017) and 29 °C (2018), while the minimum daily mean temperature was −24 °C for the three years. The year 2018 is considered to be a dry year, while 2017 is characterized as a wet year based on climate normals. The latter was coupled with an early snow melt (mid-February), which led to major flooding events in the southern regions of Quebec (NASA 2017).

Site instrumentation
The model input data are the volume of liquid water available for infiltration and the potential evapotranspiration at the chosen time step (here both calculated from a meteorological dataset; see section 'Model setup and parameterization') in addition to soil moisture as volumetric water content (VWC) for calibration. Soil moisture was measured in situ starting in August 2015 (SLZA), December 2015 (STEL) and October 2016 (SLZB) ( Table 1). Soil volumetric water content was measured with capacitance sensors (model EC-5, Campbell Scientific Inc.) located at 10, 20, 50 and 100 cm depths at sites SLZA and SLZB and at 25, 50, 75 and 100 cm depths at site STEL. VWC was recorded every 15 min and then averaged to daily values. The VWC sensors were not calibrated and the generic calibration accuracy of 0.03 cm 3 /cm 3 was used.
The weather station at STEL was equipped with various devices to collect basic meteorological data every 15 min (all at a height of 2 m, except for the anemometer, which was located at 3 m). Incoming solar radiation was measured with a Kipp & Zonen Pyranometer SP LITE2 and net radiation with a Kipp & Zonen Net Radiometer Sensor NR-LITE2, both with a temperature range of between −40 to 80 °C. Wind speed and direction were available through a Young Wind Monitors 05103-45 anemometer with an accuracy of ±0.3 m/s. Air temperature and relative humidity were recorded with a Campbell Scientific HMP60 probe with an accuracy of ±0.6 °C at an operating temperature range of between −40 and + 60 °C, and an accuracy varying between 3 and 7% for relative humidity, depending on the air temperature. Groundwater levels at the three sites were monitored hourly by automatic pressure transducers (Solinst Leveloggers 3001) in PVC piezometers (5.08 cm inside diameter) installed in the unconfined granular aquifer using a cone penetrometer at STEL (depth of 12 m and strainer length of 3 m) and a hollow auger with continuous sampling at SLZ (depth of 12 m and strainer length of 1.5 m). Finally, liquid precipitation was measured using a tipping bucket rain gauge (Campbell Scientific Inc. TB4).
Piezometers instrumented with level loggers (Solinst) were installed at the two sites. The borehole piezometer at SLZ was drilled in September 2015 and was 8.5 m deep. The sediments were characterized as medium to fine sand uniformly distributed through the soil column. At STEL, the piezometer borehole was drilled in August 2013 and the sediments were fine to medium sand from 0 to 4.1 m, silty clay between 4.1 and 6.5 m, and again fine to medium sand from 6.5 to 10 m. Below 10 m, the sediments varied between coarse and fine medium sand down to 21.6 m.

Site characterization
Soil cores were sampled at different depths at each site. Soil samples were manually collected using a shovel in spring 2019 at STEL from the 10, 20, 40, 60 and 90 cm depths, and at SLZ in summer 2018 at the 20, 40, 55, 90 and 115 cm depths for grain size analysis by laser sieving. The physical description of the textural group to which the samples belong is based on Folk (1954). The gravel fraction was under 1% for the most part, with a maximum of 4% at STEL from the soil sample at 60 cm depth. The organic matter content was estimated by ignition loss (in an oven at 500 °C) (Dean 1974;Heiri et al. 2001). Hydraulic conductivity was measured at the 20 and 50 cm depths at SLZA and SLZB, and at the 30, 60, and 90 cm depths at the STEL site using a Guelph permeameter (Soil Moisture Equipment Corp 2013).

Water-table fluctuation method
To evaluate the reliability of simulated GR, it was compared with the results obtained from the water-table fluctuation (WTF) method . The WTF method requires knowledge of the aquifer specific yield (S y ) and is based on the assumption that rising groundwater level in an unconfined aquifer is caused by percolating water reaching the water table and recharging the aquifer. The S y parameter is treated as a storage term that accounts for the instantaneous change in water storage upon a change in total head (Healy 2010). This assumes that the amount of available water in a soil column of unit surface area is S y times the height of water in the column. Because this method does not rely on meteorological data as an input, it is especially useful for a comparison with the simulated GR.

Available meteorological data
Meteorological data were available from three different sources (Table 2). First, the automated weather station located at the STEL site records local air temperature, precipitation, solar radiation and wind speed. Secondly, observed daily minimum and maximum temperature (T min and T max ) and precipitation data interpolated on a 10 km grid by Natural Resources Canada (NRCan) based on the Australian National University Spline (ANUSPLIN) interpolation method (Hutchinson et al. 2009;Hopkinson et al. 2011;McKenney et al. 2011) are also available until 2017. Finally, daily temperature, precipitation, radiation, albedo and wind speed data are available from the North American Regional Reanalysis (NARR) atmospheric and land surface hydrology dataset, which uses the very high-resolution NCEP ETA model together with the Regional Data Assimilation System (Mesinger et al. 2006). The added value of using NARR and ANUSPLIN data was investigated through hydrological modeling at the STEL station. Two different runs with varying configurations of the ANUSPLIN and NARR datasets as input weather data were performed for the 2016-2018 period. The first run, referred to as the ANUSPLIN scenario, employs T min , T max and precipitation data from the ANUSPLIN database. The missing radiation data were estimated from temperature or were set as constant (wind speed was set to 2 m s −1 and albedo set to 0.23 for grass, according to Allen et al. (1998)). The second run, denoted as combined scenario, incorporates data from both datasets (i.e., T min , T max , and precipitation from ANUSPLIN and all other variables from NARR). The latter represents a more complete dataset that maximises the number of available variables for GR modeling. Because the ANUSPLIN data have not been produced since 2017, temperature and precipitation data were from the STEL weather station for 2018 because those two variables are very strongly correlated in the data sources. A significant bias was observed for NARR precipitation data, producing a substantial underestimation of the total precipitation (similar results were observed with data from a nearby grid point). Therefore, NARR data were only used in the combined scenario to assess the effect of using albedo, radiation and wind speed from this dataset instead of using constants or estimating these terms.
Because of temporal discontinuities in the STEL dataset, the weather station data were only used for comparative purposes to evaluate the accuracy of the NARR and ANSUPLIN datasets. The comparison was made against measurements from a weather station that includes the main meteorological variables (temperature, precipitation, radiation and wind speed) needed to calculate potential evapotranspiration (PET). Because of the high spatial variability of precipitation, the closest grid point in the ANUSPLIN dataset to both SLZ and STEL sites was selected and used in the calculation to better represent the meteorological conditions at each location.

Model setup and parameterization
The physically based model HYDRUS-1D (Simunek et al. 2015) for variably saturated media was used. Using Richards' equation, the model simulates soil moisture dynamics within a soil column representing the unsaturated zone. The drainage through the root zone leaving the base of the soil column is considered to be GR. The van Genuchten- Mualem model (VGM;Mualem 1976;van Genuchten 1980) as chosen to represent the water retention and hydraulic conductivity characteristics of the soil samples with continuous mathematical functions: where θ [L 3 /L 3 ] is the volumetric moisture content; h [L] is the pressure head; θ r and θ s are the residual and saturated moisture content respectively; K [L/T] and K s [L/T] are unsaturated and saturated hydraulic conductivity, respectively; and S e = (θ -θ r )/(θ s -θ r ) [−] is the saturation degree. The fitting factor α [1/L] is inversely related to the pressure at the inflexion point of the retention curve, while n [−] is related to the pore size distribution of the soil with m = 1-1/n, and l [−] is a parameter accounting for pore tortuosity and connectivity.
The initial VGM parameters (θ r , θ s , K s , α and n) were estimated using the ROSETTA software (Schaap et al. 2001). ROSETTA uses pedotransfer functions to predict the VGM parameters using the soil textural distributions obtained from laboratory analyses for the three study sites. K s values from field tests and the predicted ones from ROSETTA were in the same range, except for the layer comprised between 16 to 30 cm depth at SLZA (Table 3). The hydraulic conductivities of the two layers at STEL, comprised between 41 and 300 cm, were higher than those predicted from ROSETTA, but they exhibited a similar increase with depth. It was hypothesized that field measurements were more representative of field conditions and local heterogeneities of the unsaturated zone than the ROSETTA predictions. For this reason, K s values for the HYDRUS model were set to fieldmeasurements at the available depth. The pore-connectivity parameter l was set to an initial value of 0.5, corresponding to an average value for many soils (Mualem 1976). Free drainage was set at the base of the soil column as a lower boundary condition. The length of the modeled soil columns was set to 3 m, with a total of 301 nodes evenly distributed between the surface and bottom.
The top boundary condition corresponds to daily vertical inflow (VI) and potential evapotranspiration (PET) used as daily values to drive the unsaturated zone model. VI is calculated from daily snowmelt values added to the liquid fraction of daily precipitation. The selected surface boundary condition allows surface runoff to occur in the model when the surface layer becomes saturated. This excess water leaves the system as runoff and is not available for infiltration. Due to the highly permeable soils at the three stations, saturation was never reached in this study and thus no runoff was simulated, but this could not be verified with field data. The simulated GR is thus considered to be an upper limit for GR, considering that runoff could occur in situ during high intensity precipitation events. Because precipitation data is available as total water equivalent, the calculation of VI as input to the snow model included separation between solid and liquid precipitation, as suggested by Turcotte et al. (2007): where T min and T max are the daily maximum and minimum temperatures (°C) and SnowFrac is the snow fraction for daily precipitation events.
A degree-day model was used to assess daily snowmelt available for infiltration: where Melt is the daily snowmelt (mm/day), C melt represents the snowmelt rate (mm/°C/day), T air is the mean daily air temperature (°C), and T melt is the temperature at which the snow starts to melt (was set to 0 °C).
The snowpack density and depth were retrieved from a nearby governmental weather station (MELCC 2019a) to calibrate the snowmelt coefficient from the degree-day model and simulate the evolution of the snowpack during the winter seasons from 2000 to 2017. The daily calculated snowmelt values were added to the liquid fraction of daily precipitation to generate vertical inflow values (VI). Daily PET was calculated using the Penman-Monteith equation (Allen et al. 1998). Beer's law (Ritchie 1972) was used to partition PET into potential evaporation (E p ) and transpiration (T p ) directly in the model: where k is an extinction coefficient and LAI is the leaf area index (L 2 /L 2 ). LAI data were obtained from the MODIS_MCD15A3H dataset, with a spatial resolution of 500 m × 500 m at 4-day intervals (Myneni et al. 2015). Daily LAI data for each site were obtained by linear interpolation between these intervals, coupled with a 30-day window moving average. LAI was used as a primary control of PET among different ecosystems in the same ecozone, such as the forest and pasture (Zha et al. 2010). Therefore, LAI data at the SLZ site were extracted from two points close to the site but each more representative of their respective ecosystem (forest and grassland).
The root water uptake was computed using the Feddes et al. (1978) model: where α(h) is a dimensionless function varying between 0 and 1, depending on soil matric potential, which corresponds to the force with which water is held within the soil matrix, and S p [1/T] is the potential root water uptake and assumed to be equal to T p .
The distribution of S p over the root zone depends on root density distributions (between 0 and 1) attributed to each site and was selected based on literature descriptions of the vegetation's physiological characteristics (Rudolph 1985;Wang et al. 2009b). In the absence of in situ measurements of root density, a rectangular distribution (homogeneous density distribution with depth) or triangular (linear density decreases with depth) profile are generally recommended. For example, 3-m-deep sandy soil cores analyzed for root  Wang et al. (2009b) showed that in sandy soil cores from a grassland environment, 60-70% of the total root biomass occurred in the top 20 cm depth. At the STEL and SLZA sites, the root density was assumed to be linearly distributed between 1 at the surface to 0 at a depth of 30 cm. Jack pine trees generally develop a lateral root system, and the bulk of the root system is largely confined to the upper 45 cm of the soil, and mostly in the top 15 cm (Rudolph 1985). Therefore, at the SLZB site, the root density was set equal to 1 from 0 to 30 cm depth and a density ranging from 1 to 0 was used between 30 and 45 cm. The root water uptake was then assumed to be equal to actual transpiration. The actual evapotranspiration (AET) was calculated based on the water availability from PET and was the sum of the actual soil evaporation and actual transpiration rates. The soil matric potential values for delineating root water uptake were taken from the database integrated into the HYDRUS-1D model and assigned as alfalfa. Precipitation interception (by plant canopies before reaching the soil) was considered negligeable.

Model calibration
The model was used to simulate daily VWC and those were compared to measured values at the four depths where soil moisture sensors are located. The models were calibrated between 2015 and 2018 (depending on the monitoring period of each site), using the year 2015 as a warm-up period. Because of the limitations of the various meteorological observations, multiple data sources have been combined to calibrate the unsaturated zone model (Shen et al. 2010;Maggioni et al. 2014). The soil hydraulic parameters (θ r , θ s , α, n, l and K s ) were calibrated using the default Marquardt-Levenberg type parameter optimization algorithm (Marquardt 1963) implemented in the HYDRUS-1D model to reproduce measured daily VWC at different soil depths at the three sites. The hydraulic parameters of each soil material were optimized successively using the field monitoring data, here starting with the top layer downwards until there was no further improvement in the objective function (Hopmans et al. 2002). The calibration period excludes the winter months (January to March) because of uncertainties regarding the quality of the VWC measurements during the freezing period. To evaluate how the choice of different meteorological datasets can influence GR estimates, the calibrated model at the STEL site was run for the three meteorological scenarios described above. In many studies, θ r is not calibrated (Thoma et al. 2014;Turkeltaub et al. 2015) because GR estimates are found to be insensitive to this parameter (Simunek et al. 1998;Scharnagl et al. 2011). In the current study, the VWC at the three sites varied mostly within the dry range (VWC closer to θ r than θ s ) which made it difficult to calibrate θ s at a daily time step (rapid drainage and saturation state might not be properly captured). This parameter was thus calibrated only for the first soil layer at STEL and SLZB, where higher VWC were measured (see Fig. 4 and Fig. S2 of the electronic supplementary material (ESM)), but θ r was calibrated because simulated VWC remained close to this value. In the literature, the pore connectivity parameter l in the VGM model is often set to a value of 0.5 (Mualem 1976;Zhu et al. 2013;Turkeltaub et al. 2015) considered to be an average value for many soils (Simunek et al. 2015). Flow in the unsaturated zone was shown not to be very sensitive to this parameter in general (Wang et al. 2009a). However, because sensitivity was demonstrated for very low VWC similar to those encountered in the current study (Vrugt et al. 2001), this parameter was calibrated here as suggested by Scharnagl et al. (2011). The upper and lower calibration bounds for all parameters (Table S3 of   were estimated based on the observed VWC at each site or on literature values for coarse materials. The calibrated models were used to simulate recharge between 2004 and 2015, with the year 2003 as the spin-up period. Meteorological data from the combined scenario that incorporates both NARR and ANUSPLIN datasets were used for the three sites. These simulations assumed that land use did not change during this period and were used to study the long-term GR processes and for comparison with other published results.

Comparison of meteorological datasets
Because of temporal discontinuities in the STEL dataset and the short temporal coverage, the weather station data were only used for comparative purposes to evaluate the accuracy of the gridded meteorological datasets NARR and ANUSP-LIN. Comparing ANUSPLIN and NARR data with those of the STEL weather station showed a strong correlation between the datasets for temperature, with R 2 = 0.99 and slope = 1 for ANUSPLIN and R 2 = 0.95 and slope = 0.98 for NARR ( Fig. 2a and b). The precipitation data from ANUPS-LIN had a very similar distribution to the STEL data in 2017 ( Fig. 2c and d) and exhibited similar results for 2016 (data not shown). A significant bias was observed for NARR precipitation data, which significantly underestimated the total precipitation (similar results were observed with data from a nearby grid point). The reason for this remains unclear, but could be attributable to the fact that there has been no assimilation of precipitation with observed data over Canada since 2003 (Mesinger et al. 2006). From 2016 and 2018, the average difference between total annual precipitation measured at the STEL station compared to the NARR dataset was 332 mm/y. However, this difference is an underestimation because there were days without measurements at STEL in 2016 and 2017 and the comparison was made on equivalent datasets (days without measurement at STEL were also removed from NARR dataset). The agreement between net radiation from the NARR dataset and that measured at STEL was somewhat scattered but reasonably good (R 2 = 0.82 and slope = 1.01), with lower agreement for the calculation of net radiation from temperature values (R 2 = 0.77 and slope = 0.95) (Fig. 2f and e). Both net radiation results showed an overestimation, particularly during the summer months. The overestimation of radiation in the reanalysis dataset may be explained by the uncertainty linked to the small scale of the study area compared to the larger scale (32 km grid) of the NARR dataset (e.g., cloud cover fraction and forest cover limiting solar exposure).
An overview of the hydroclimatic conditions across the study area is given as a statistical summary of the calculated annual mean VI and PET for each site (Table 4). VI values are similar at the three sites, except in 2016, when they were 99 mm lower at STEL (935 mm) than at the SLZ station (1034 mm; different ANUSPLIN point grid) and considered closer to normal conditions than the other two years. Also in 2016, the PET for the SLZ and STEL sites were 680 and 684 mm, respectively. The year 2017 can be characterized as a wet year, with a significant amount of VI (1208-1205 mm for SLZ and STEL). The PET during 2017 was 677 and 681 mm for SLZ and STEL, respectively. The year 2018 can be considered to be a dry year, with 829-830 mm of VI and a PET of 701 mm for SLZ and STEL. PET was the highest in 2018 for all sites because it was the hottest year (maximum daily mean temperature 26 °C in 2016 and 2017, and 29 °C in 2018).

Soil hydraulic characteristics
A textural analysis was performed on samples extracted from up to 90 cm depth at STEL and 120 cm at SLZ (Table 3). For modeling purposes, the results of the deepest layer were considered representative of the conditions below that layer down to the bottom of the soil column (300 cm). None of the samples contained a clay fraction. The STEL site showed the most drastic change in soil texture with depth. The first layer, from 0 to 40 cm, was rich in silt (13.2% compared to 7.4% for SLZ sites), with a K s of 362 cm/day obtained with the Guelph permeameter. The grain size shifted to a silt fraction of 5.3% and a K s of 1296 cm/day between 41 and 65 cm depths. This is consistent with field observations of a very dense and fine soil horizon between approximately 20 and 50 cm. The two prairie sites showed the least organic matter, with 2.7% at SLZA in the first 15 cm depth and 2.5% in the first 40 cm depth at STEL. In the forest, organic matter content decreased from 6.9% in the first 15 cm depth to 3% in the layer from 31 to 60 cm depth. At SLZ, the soil characteristics of sites A and B are considered to differ only in their percentage of organic matter content and root distribution, here based on their proximity (15 m) and field observations. Therefore, the same textural analysis results (but not the same percentage of organic matter) were used to simulate the two sites. An important difference can be seen in the measured K s at the 20 cm depth for both sites, with a value of 1036 cm/day in the prairie and 240 cm/day in the forest. This difference illustrates the influence of organic matter content and vegetation on infiltration at shallow depths. It is consistent with field observations, where unconsolidated medium sand at SLZA and a rich organic matter horizon in the forest were observed, enhancing water holding capacity at shallow depths.
Because organic matter content appears to play a significant role in hydraulic properties but is not considered in ROSETTA's pedotransfer functions, the organic matter content was integrated into the soil textural distribution as a silt fraction. This reflects the fact that moisture storage capacity is positively correlated with organic matter and silt (Bot and Benites 2005). Wang et al. (2009b) found that increasing soil organic matter induced lower hydraulic conductivity in the sandy soils of the Nebraska Sand Hills, suggesting that soil carbon should be considered when estimating the hydraulic conductivity. In the current study, the modification of the textural distribution to include organic matter contents resulted in initial parameters better representing the hydraulic properties of the upper soil layer.

Measured soil moisture contents
The VWC measured by the uppermost sensor at SLZA (Fig. 3a) exhibited a weak dynamic behavior, with a few small variations. This contrasted with the VWC measured at the same depth at SLZB, which was highly dynamic, covering a large range of moisture contents (Fig. 4a). Measured VWC at 10 cm depth in the forest site (SLZB) varied between 0.04 and 0.40 cm 3 /cm −3 (mean 0.15 cm 3 / cm 3 ), while in the prairie site (SLZA), it varied between 0.02 and 0.19 cm 3 /cm 3 , (mean 0.12 cm 3 /cm 3 ) (Table S1 of the ESM). At the STEL site, two very low VWC periods were recorded by the uppermost sensor in August 2016 and August 2018 (Fig. 5a). Generally, the mean soil moisture content measured by the uppermost sensor (25 cm) remained high, between 0.24 and 0.37 cm 3 /cm 3 except for the two periods of very low VWC. At this site, the soil needed particularly dry conditions for many consecutive days to reach a low water content. The measured VWC of the deeper sensors generally remained within the same narrow range of variation at every site, between 0 and 0.15 cm 3 /cm 3 at SLZB, between 0 and 0.10 cm 3 /cm 3 at SLZA, and between 0.04 and 0.08 cm 3 /cm 3 at STEL. There was no noticeable interannual difference in the range of variability of VWC values during the study period, particularly for the deeper sensors at the three sites. An exception was the uppermost sensor at SLZB (Fig. 4a), which showed highly variable VWC throughout all of spring and summer 2017, while in 2018 (dry year), the VWC was rapidly in the dry range at the beginning of the spring and showed less variability until fall (no observed data available for summer 2016).
At the three sites, the VWC generally decreased with depth, and becomes increasingly lower, less temporally variable and closer to the residual water content (θ r ). A seasonal trend of VWC between April-May and August can be observed at SLZB down to 50 cm depth in 2017 and 2018. A similar trend can be observed at 25 cm between April-May and August in 2016 and 2018 at STEL. This pattern was absent at SLZA and at deeper locations.

Simulated soil moisture contents
The calibrated VGM parameters for the three sites are reported in Table S4 of the ESM, while the initial values and calibration bounds are presented in Tables S2 and S3 of the ESM. The simulated daily soil moisture contents for the three sites (Figs. 3, 4, and 5) represented the measured values at the four instrumented depths relatively well. The mean absolute error (MAE) values were below the accuracy of the soil moisture sensors (0.03 cm 3 /cm 3 ), indicating that the simulated long-term soil water storage matched the observed ones well (see Figs. 3b, 4b and 5b). The mean error (ME) was nearly 0 at all sites, demonstrating that there was no tendency for over-or under prediction. The root mean square error (RMSE) was almost equal to half the standard deviation (sd) at SLZB (sd/2 = 0.030; RMSE = 0.032), showing an acceptable agreement with the experimental data (Singh et al. 2004). The model at STEL exhibited the best fit with a RMSE of 0.014 compared with 0.054 for half the standard deviation, while at SLZA, both values were equal (0.018). The R 2 values showed the same tendency, with the best goodness-of-fit at STEL with a value of 0.98, followed by SLZA (R 2 = 0.75) and SLZB (R 2 = 0.71), showing acceptable agreement between the observed and simulated moisture content.
The observed and simulated VWC at the 10 cm depth at SLZA for the calibration period showed that the simulated values overestimated the measured ones in 2016 and underestimated them in 2017 (Fig. 3). This might be attributable to the evolving compaction of the sand after the installation of the sensors. The simulated VWC values at this site showed more rapid variations than the measured values at 10 and 20 cm depths. The simulated soil column at this site were very reactive to the atmospheric boundary conditions, which is typical of coarse sediments coupled with low organic matter content (Table 3). This might indicate the need to integrate an additional superficial fine soil layer in the first centimeter depth to dampen the infiltration rate and better represent natural conditions with vegetation. Thus, the natural soil water dynamic could not be described only by the estimates of the ROSETTA-calibrated pedotransfer functions as reported by Saïd Ahmed and Coquet (2018). The VWC simulated at SLZB covered the large range of measured moisture contents (Fig. 4) and reproduced relatively well the moisture dynamic and seasonal trend in measured soil moisture content. The simulated VWC from the calibrated model at STEL showed good fit with the experimental data, though the magnitude of the two August droughts was not completely simulated in 2016 and 2018 (Fig. 5). This might be caused by underestimated evaporation in the topsoil layers in the modeled soil column or misrepresentation of root density with depth. In general, the resulting simulated soil moisture contents corresponded well with the dynamic of the VI inputs and were able to capture most of the observed soil moisture peaks and drainage events.

Calibration period simulation
Interannual variability in GR was observed at each site, in agreement with the variation in annual VI and AET. The simulated recharge rates ranged from 344 (SLZB 2018) to 703 mm/y (SLZA 2017) (Table 4). This large difference is caused by the annual variation in precipitation and the contrasted effect of vegetation and organic matter content on GR and AET. Mean GR recharge for 2016-2018 at the three sites was 517 mm/y. The STEL and SLZB sites were found to have similar GR, with a maximum of 575 (STEL) and 590 mm/y (SLZB) in 2017 and a minimum of 374 mm/y (STEL) and 344 mm/y (SLZB) in 2018. SLZA produced a slightly higher range of GR values, with a maximum of 703 mm/y (2017) and a minimum of 486 mm/y (2018). Similar observations of interannual variability can also be made for the simulated AET, which varied from 368 mm/y (SLZA 2018) to 620 mm/y (SLZB 2017). The maximum AET was obtained in 2017 (wettest year) at every site. The water not used for AET or GR corresponded to increased storage in the soil column. The GR computed for STEL is lower than that for the SLZA site. This was unexpected because these sites have similar land cover, but can be explained in part by the higher AET at STEL (493 to 620 mm/y) than at SLZA (368 to 460 mm/y) for the calibration period. The lower GR simulated for STEL in 2016 than for the two other sites can be explained by the precipitation values being 99 mm lower at this site in 2016 (different meteorological grid point; see Table 4).

Recharge estimated with the water table fluctuation method
One of the main uncertainties in the GR estimates from the WTF method arise from the difficulty in estimating S y . Loheide et al. (2005) showed that coarse-grained sediment S y are generally similar to the difference between θ s and θ r . The calibrated VGM parameters were thus used Fig. 2 Comparison of meteorological variables from the Australian National University Spline (ANUSPLIN) dataset in the left column and the North American Regional Reanalysis (NARR ) dataset in the right column against data monitored at the STEL site. a and b show linear regressions between mean daily temperature, c and d cumulative precipitation in 2017, e and f linear regressions between the net radiation from NARR and that calculated from temperature from the ANUSPLIN dataset to estimate S y , leading to values of 0.33 and 0.34 for the SLZ and STEL sites, respectively. These values are comparable with those from the literature: between 0.32 and 0.39 for sandy soils (Loheide et al. 2005), from 0.11 to 0.36 (Cohen 1965), and from 0.20 to 0.35 for a coarse sand (Johnson 1967). Lower S y values are also reported in the literature (e.g., 0.17 for a gravelly sand; Stephens et al. 1998). To better represent the possible range of S y values and assess the resulting uncertainty in the WTF method, GR was also calculated with a S y of 0.15 and 0.38 (Fig. 6). The results showed that simulated GR represented the higher bound of WTF GR estimates.
The water levels at both sites exhibited a well-defined seasonal cycle of rapid rise in spring followed by a gradual decrease throughout the remainder of the year, although the water level at STEL is more subject to daily variations (see Figs. S1 and S2 of the ESM). Field observations at STEL showed a silty soil layer with traces of sand and clay near the water Table (4 to 6.5 m). This lower permeability layer could explain the more dynamic behavior of the water table at the STEL site (average water table depth of 5.5 and 6.8 m at SLZ and STEL, respectively).
The GR results obtained by WTF show the highest recharge rate for both sites in 2017 and the lowest in 2018 (Fig. 6). This observed trend follows the variations in annual precipitation. Therefore, when the two sites are compared, the GR is higher at STEL in 2016 and 2017, even if precipitation at this site in 2016 was lower by 99 mm. GR for SLZ is higher than that computed for STEL in 2018. The water table levels were the lowest at both sites in 2018 (Figs. S1 and S2 of the ESM), probably because of drier conditions.

Simulated recharge between 2004 and 2015
The calibrated HYDRUS model was used to simulate GR between 2004 and 2015 with the combined scenario (mean VI of 975 and 954 mm/y at SLZ and STEL). The

Fig. 3
Observed and simulated soil volumetric water content (VWC) during calibration at different soil depths at the SLZA site and statistical results of the calibration performance. The red and blue lines represent fitted and observed daily mean soil moisture content values, and the gray areas indicate the accuracy of the sensors results showed mean long-term GR rates of 547, 410, and 408 mm/y at SLZA, SLZB, and STEL, respectively (inter-site average 455 mm/y). These values represented the mean ratios of GR/VI of 56% for SLZA and 42% for SLZB and STEL, respectively, ranging from 26 to 64% of annual VI following meteorological variations. The mean AET/VI ratios for this period were of 42, 56, and 55% for SLZA, SLZB and STEL, respectively.
The monthly GR rates showed strong annual and seasonal variations (Fig. 7). The main recharge peaks occurred during the spring, with a maximum value in April (means between 99 and 108 mm/month), before declining rapidly throughout the summer period and were the lowest in August and September (means between 4 and 30 mm/month) as evapotranspiration increased and the antecedent soil moisture content decreased. The magnitude of the estimated recharge in March and April fluctuated because of variations in the start of the snowmelt period and in the winter snow cover. These observations indicated the significant role of the snowpack on spring GR. Another distinct but less important recharge period was observed between October and December (means between 37 and 69 mm/month). At the three sites, spring GR represented a mean of 37% for SLZA, 43% for SLZB and 44% for STEL, while in fall, it represented 29% at SLZA and 31% at SLZB and STEL. It should be noted that these ratios were highly variable between years. Recharge events also occurred during the winter months of January (means between 40 and 43 mm/month) and February (means between 16 and 17 mm/month).
The minimum mean GR was obtained in July, August or September for the three sites (Fig. 7) even if PET and VI in May and June were in similar ranges. Hence, those five months should have similar GR variability, especially in May and September because of their very similar PET and VI values. However, the variability in monthly GR was significantly higher in May for every site (see Fig. 7). This is likely due to the antecedent soil moisture content Fig. 4 Observed and simulated soil volumetric water content (VWC) during calibration at different soil depths at the SLZB site and the statistical results of the calibration performance. The red and blue lines represent fitted and observed daily mean soil moisture content values, and the gray areas indicate the accuracy of the sensors during the preceding month. Late spring and early summer GR were preceded by the two wettest months (March and April), while September was preceded by the two driest months. Recharge generally increased in October when PET became significantly lower than VI. No significant correlation between annual GR values and VI values was found (not shown here). This strongly suggests that the recharge fluxes were mainly influenced by relatively recent precipitation patterns and antecedent moisture content of the one to two preceding months.

Soil properties controlling volumetric water content
The contrasted moisture dynamics between SLZA and SLZB observed at shallow depths cannot be caused by the soil textural distributions alone. Therefore, the high organic matter content observed in the forest was considered to be an important soil characteristic for representing the soil moisture behavior in the top layers of these permeable sediments. The organic matter content or silt fraction in the uppermost soil layers can explain in part the higher and more variable VWC measured over time within the first 30 cm at every site and the consistently lower and less variable VWC with increasing depth. These observations are the result of the enhanced soil water holding capacity because of soil characteristics (fine fraction and presence of organic matter at shallow depths), in addition to the tighter coupling between soil moisture and land surface processes at shallow depths that directly receive precipitation and are subject to more evapotranspiration (Martínez-Fernández and Ceballos 2003;Guber et al. 2008). At the study sites, the sand fraction increased with depth (Table 3), promoting rapid percolation and lower water holding capacity.
The relationship between moisture and soil characteristics was also demonstrated by the measured VWC dynamic at The red and blue lines represent fitted and observed daily mean soil moisture content values, and the gray areas indicate the accuracy of the sensors each site at shallow depths. The lower variability of VWC at SLZA (Fig. 3) is typical of coarse sediments characterized by low holding capacity and favoring rapid drainage. It contrasted with the VWC measured at SLZB, which was highly variable and covered a large range of moisture content (Fig. 4). These observed conditions supported the integration of the organic matter into the silt fraction to account for field differences not reflected in the soil textural distribution alone.
At the STEL site, the soil needed particularly dry conditions for many consecutive days to reach low water content and create drought conditions, as observed by the first sensor in 2016 and 2018 (Fig. 5a). Hence, a higher AET demand was necessary to evacuate the pore water in the denser and finer soil at STEL. These soil characteristics allowed for a higher water holding capacity, comparable to the SLZB site, but with less dynamic variations. Also, the higher evapotranspiration associated to the vegetation present at SLZB must have played an important role in this dynamic behavior of VWC. Therefore, the higher natural silt fraction in the first soil layer at STEL permitted higher water retention and slow drainage. The higher fraction of organic matter in the forest also enhanced the water holding capacity, similarly to STEL, but allowed the soil to drain more easily when mixed with mostly medium sand. Other studies have also shown that many factors influence soil hydraulic characteristics, including bulk density and pore size distribution (Schaap et al. 2001;Wang et al. 2009b). This points to the importance of adequately characterizing the soil properties to simulate VWC and GR with an unsaturated zone model.
The model was able to capture the major processes in the soil moisture dataset obtained from an important portion of the unsaturated zone, extending from the ground surface to well below the root zone. Factors that might contribute to the deviation between observed and simulated VWC may come from a mismatch between the recorded precipitation or evapotranspiration rates, from uncertainties in the observed soil moisture data, or from the unaccounted heterogeneity of the soil layers, such as the presence of macropores. For example, the simulated VWC from the calibrated model at SLZB showed the lowest calibration statistics, which could be caused by a greater presence of natural heterogeneity in the forest ground, which was more likely to create minor local discrepancies between the measurements and simulation. Simulation errors can also be attributed to the significant spatial variability of soil moisture and degree to which the sensors were in contact with the soil material, especially in coarse soil. Nonetheless, the calibration statistics obtained (RMSE between 0.014 and 0.032; R 2 between 0.71 and 0.98; ME of nearly 0) were within the general ranges of values found in the literature (Jiménez-Martínez et al. 2009;Assefa and Woodbury 2013;Min et al. 2015;Wang et al. 2016), indicating good model performance.

Soil properties controlling recharge
Soil characteristics were shown to influence soil moisture dynamics and, subsequently, recharge processes. The higher GR and lower AET at SLZA were attributed to the coarser nature of the sediments in the root zone and to the vegetation type, which were both associated with a low soil waterholding capacity. Also, preferential flow was not considered in the model and could had an influence on the GR as water can reach the water table faster. With preferential flow, precipitation arriving at the surface typically infiltrates almost instantaneously because of the low organic matter content and can percolate rapidly through the root zone to recharge the underlying aquifer. Water flow of this nature limits plant uptake. The lowest GR was computed for SLZB, which was expected because the site was under forest cover and should take up more water in the soil by transpiration than prairie environments because of the systematically higher AET. The highest GR was computed at SLZA site (703 mm/year) and was probably out of bound. The GR computed at STEL was significantly lower than at SLZA, which was surprising because they had similar land cover. Because of the finer soil texture at STEL in the first 40 cm of depth, water flow was slowed, and the water-holding capacity was increased, permitting higher AET rates compared with the other prairie site. This explained the similar GR estimated at STEL and SLZB, which exhibited similar soil characteristics in the first Fig. 6 Groundwater recharge estimation results from HYDRUS-1D and the water-table fluctuation (WTF) method for the calibration period at the three study sites centimeter depth in the model (Table 3; the higher organic matter content at SLZB here was transformed in silt fraction in the model). Thus, the soil textural distribution in the root zone was an important factor controlling the flow processes in the unsaturated zone.
The analysis of water level data to quantify GR using the WTF method yielded higher GR rates than those simulated by HYDRUS at STEL in 2016 and 2017 (Fig. 6). The interpretation of the numerous small increases as recharge events at STEL in 2016-2017 (wettest years) might have resulted in uncertainty propagation in estimating S y or recharge events, leading to exaggerated GR rate estimations for those years. Some artifacts can arise in the water level time series from logger inaccuracy and from the Lisse effect, which affects the water level in unconfined aquifers by the entrapment of air between the water table and a wetting front (Weeks 2002;Guo et al. 2008). The Lisse effect is normally prevalent in fine-textured soils with shallow water tables <1-1.3 m (Weeks 2002), which was not the case here. The WTF method is also prone to errors when lateral groundwater inflows from upgradient areas are different from groundwater downflows to downgradient areas, causing potential artificial variations of the water table level that are not directly due to local vertical recharge. The water table time series at SLZ showed a less dynamic behavior, probably leading to fewer misinterpreted recharge events during wet years. This induced a lower estimate GR for STEL with the WTF method. It should also be noted that the SLZA site only represented a clearing area in a large forested zone where PET was probably high. This could explain the similar GR estimated by the WTF and model in the forest because the groundwater level fluctuations were representative of regional hydrogeological conditions. It could also be attributed to an underestimation of the recession curve used in the WTF method because of steady infiltration process in the coarser soil at SLZ. Also, the results from WTF obtained for STEL gave a similar, although slightly lower GR estimate than the HYDRUS model at the other prairie site (SLZA) in 2016-2017. This might be an indication of an overestimation of the fine-textured soil control in simulated GR in top soil layer at STEL.
Despite the uncertainties associated with S y and the interpretation of recharge events, the WTF method captured the whole interannual range of GR variability simulated by the model (Fig. 6). The year 2017 was the wettest, leading to the highest GR, and 2018 was the driest, with the lowest GR for each site. Consequently, this approach represents a useful tool for estimating the possible range of GR independently of meteorological data or when only few data are available. The range of WTF GR values for different S y (0.15 to 0.38) framed relatively well the simulated GR but illustrated that GR could be overestimated both by the WTF and by the HYDRUS model. Results should thus be considered as possible maximum ranges of GR. It should also be noted that the baseline S y was calculated with two optimized parameters from HYDRUS and for this reason WTF GR estimates are not completely independent from those from HYDRUS. The HYDRUS-1D model provides Fig. 7 Monthly median, and 25th and 75th percentile results for HYDRUS-simulated groundwater recharge from 2004 to 2015 at the three study sites: a SLZA, b SLZB and c STEL. The blue and the red lines represent monthly mean vertical inflows (VI) and potential evapotranspiration (PET), respectively, for the simulation period. The black points represent outlier values insights into processes taking place in the unsaturated zone at the site scale, while the WTF method provides information at the regional scale and does not rely on meteorological data. Therefore, GR estimates from the two methods should be used with caution to estimate local GR rates at other locations within the area. The main advantage of using HYDRUS is that it is based on flow equations in the unsaturated zone, a reservoir that is still poorly understood, yet crucial for recharge. Another advantage of this type of model is that is uses data from soil moisture sensors that are relatively easy to install and maintain at low cost in a wide network. As pointed out in many studies, it is important to combine different GR estimation methods to reduce the uncertainty on GR estimation (Scanlon et al. 2002;Rivard et al. 2014;von Freyberg et al. 2015).

Long-term groundwater recharge simulation
Based on simulations from a spatially distributed surface water budget model (1989( ), Larocque et al. (2015 found that the areas of unconfined granular aquifers in the region received a mean GR of 331 mm/y, reaching a maximum of 440 mm/y. This maximum was lower than the mean GR rate found by HYDRUS-1D. The spatially distributed surface water budget model partitioned the water budget into multiple components not accounted for in the current study (runoff and subsurface flow) and considered negligible at the local scale because of the highly permeable soils (between 84 and 100% sand) and flat topography. These characteristics promoted rapid infiltration and prevented runoff and subsurface flow during rain events at the local scale. However, this assumption might have overestimated GR values because most of the infiltrated water reached the bottom of the soil column. Larocque et al. (2015) estimated that the runoff and subsurface flow components accounted for 56% of VI, corresponding to a mean of 540 mm/y for the whole region. However, the model from Larocque et al. (2015) was spatially distributed over the entire Vaudreuil-Soulanges region where topography was not flat and highly permeable sediments were not present everywhere. In the same study, GR from March to May was found to represent 38% of the annual recharge, and in 44% of the annual GR during the fall (October to December). These water budget components generated estimates that fell within the range of GR and GR/VI ratios computed by HYDRUS-1D for the three stations. Larocque et al. (2015) showed for the study area that mean annual AET represented 39% of VI using a surface water budget, and 46% of VI with the conceptual MOHYSE model . These ratios were lower than the findings of the present study but still within the same range of magnitude. January and February GR estimates from Larocque et al. (2015) represent 9 and 5 mm/month, respectively. For the same year, the present study simulated GR ranges from 22 to 30 mm/month in January and from 10 to 14 mm/month in February for the three sites. These simulated GR ranges appeared to be relatively important (especially for January). This could be due to the extremely local simulation of GR in HYDRUS-1D model, compared to the large-scale spatially distributed water budget model of Larocque et al. (2015) which did not considered explicitly soil characteristics. Barbecot et al. (2018) estimated GR at STEL from 2010 to 2013 with a hydro-isotopic water budget and found a GR of 200-372 mm/y occurring exclusively during the snowmelt period, with an AET of 608-517 mm/y. At STEL, for the same period, a mean GR of 400 mm/y and AET of 542 mm/y were calculated with HYDRUS-1D. The results from the hydro-isotopic water budget and from the model used in the present study were within the same order of magnitude and in good agreement. Mattei et al. (2020a) estimated GR in the same region using water isotopic profiles. The estimated GR at SLZ for 2017 varied between 380 and 395 mm/y and from 158 to 214 in 2018. They used the NARR meteorological dataset, which, as shown here, significantly underestimated precipitation values (an average of 331 mm/y as calculated in this study from 2016 to 2018). This might explain the lower GR estimation in Mattei et al. (2020a).

Timing of recharge and evolution under climate change
As seen in Fig. 7, the interannual variability of monthly GR was largest during the spring and fall when PET was negligible and with snow cover conditions being very different from year to year. Fall evapotranspiration and the timing of the spring snowmelt controlled the timing of the preferential recharge periods, while precipitation mainly controlled the magnitude and interannual variability of GR rates. This seasonality of GR processes might be an important factor to consider in the southern areas of Quebec, where precipitation is expected to increase under climate change (Ouranos 2015). It is possible that increased precipitation will be counteracted by increased PET resulting from higher temperatures, especially if higher precipitation occurs during the summer period. Conversely, if more precipitation occurs in spring, when AET is equal to PET and groundwater levels are high, climate change could lead to more runoff and eventually more flooding events. Milder winters might also be an important issue to consider in a climate change context because this will influence the actual seasonality of GR (Allen et al. 2010) by producing episodic recharge events during the winter period and reducing spring recharge volume (Rivard et al. 2009).

Impact of the meteorological dataset
Two different model runs using varying configurations of the ANUSPLIN and NARR datasets as input data were performed at the STEL site to assess the impact of their utilization on GR estimates. These runs were performed at the STEL site with the model calibrated using the combined scenario as a sensitivity analysis. The GR computed with the ANUSPLIN scenario (398, 535, and 292 mm/y) was systematically lower than that calculated with the combined scenario (446, 575, and 374 mm/y), even if the same VI was used as an input. This was because of the PET calculation in the ANUSPLIN scenario, which was computed using the temperature data as a proxy for radiation terms. Using the NARR dataset in the combined scenario enabled the calculation of net radiation without approximating it from temperature data. It has been shown that the radiation term has the most impact on PET calculation (Sentelhas et al. 2010). Even though radiation data from NARR were overestimated during the summer compared with the measurements made at the STEL station, they yielded lower PET values (closer to those calculated with STEL measurements), resulting in higher GR rates for the combined scenario. The NARR data significantly underestimated precipitation and were therefore not used as a single input meteorological data source for GR simulation. Similar results were obtained by Wong et al. (2017) who intercompared several gridded precipitation products over 15 terrestrial ecozones in Canada for different seasons from 1979 to 2012. The overall reliability of NARR was found to be low, probably because the absence of precipitation assimilation over Canada since 2003 (Mesinger et al. 2006). They reported that ANUSPLIN performed well in capturing the key timings and in minimizing the error magnitude of the precipitation. Langlois et al. (2009) ran three snow models with NARR data in the southern regions of Quebec during three winters between 2004 and 2008. All the models gave accurate results compared with the simulations driven by data from a ground-based station or measurements of snow water equivalent. They demonstrated that the regional reanalysis can be used in snow models to predict snow water equivalent in Quebec. Choi et al. (2009) also demonstrated the applicability of NARR temperature and precipitation data for modeling hydrological processes in the northern regions of Manitoba. Therefore, for areas without direct measurements of radiation or AET, the NARR dataset may provide valuable additional information for the calculation of PET, and the underestimation of precipitation might differ by location and period of study.

Conclusion
The HYDRUS-1D model of vertical flow in the unsaturated zone was used to simulate GR rates at three experimental sites in the preferential recharge areas of the Vaudreuil-Soulanges region of southern Quebec (Canada). The results obtained were in agreement with those from the WTF method and those from previous researches which used different scales and approaches. In addition, the results demonstrated the strong influence of climate and soil characteristics over recharge processes in the Canadian humid continental climate. Overall, the model was able to capture the major phenomena revealed in the monitoring data and with the induced vertical inflow dynamics. The grain size analysis could not on its own account for specific field conditions particularly in the forested site, where more organic matter content was present and was hypothesized to have a significant influence on soil moisture. The inability of the model to fully describe the experimental data stemmed from the gap between model representation and a lack of incorporation of more complex processes and structures present in nature. The results also showed the importance of soil properties in controlling volumetric water contents, percolation and GR.
The lack of complete and reliable meteorological datasets is often a limiting factor in hydrological studies around the world. The accuracy of such datasets is not systematically questioned prior to using them for hydrogeological modeling and was found to have a major impact on GR estimation. For areas without direct measurements of radiation or AET, the NARR dataset might provide additional valuable information for the calculation of PET, and the underestimation of precipitation might differ by location or time period considered. Comparing NARR and ANUSPLIN data with those of a local weather station should be carried out with caution. In other regions, the quality of the NARR dataset could be better than suggested here. These results provide information to potential users regarding the performance of different meteorological datasets in the southern regions of Quebec for hydrogeological modeling specifically.
The current study offers guidance to groundwater modeling practitioners by highlighting the main controls of recharge and the processes occurring in the unsaturated zone in a Canadian context. The unsaturated zone is omitted in many studies but was found to be of great interest for understanding important processes herein. As in every model, the unsaturated zone model that was used could not fully consider every process and mechanism occurring in the natural environment (e.g., preferential flow) but could still be suitable in other geological or climatic contexts. Processes as surface and subsurface runoff, even in nearly flat topography might play an important role in GR assessment by limiting the amount of water recharging the aquifer.