# 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

Sabrina Bruneau, Florent Barbecot, Marie Larocque, Viorel Horoi, Yves Coquet, Sophie Guillon

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sabrina Bruneau, Florent Barbecot, Marie Larocque, Viorel Horoi, Yves Coquet, et al.. 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. Hydrogeology Journal, 2022, 30, pp.653 - 672. 10.1007/s10040-021-02429-6. hal-04153174

HAL Id: hal-04153174

https://agroparistech.hal.science/hal-04153174

Submitted on 6 Jul 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 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

Sabrina Bruneau<fnref n="1" />,<fnref n="2" />* · Florent Barbecot<fnref n="1" />,<fnref n="2" /> · Marie Larocque<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> · Viorel Horoi<fnref n="1" /> · Yves Coquet<fnref n="4" /> · Sophie Guillon<fnref n="5" />

Received: 25 April 2021 / Accepted: 17 November 2021

© The Author(s) 2021

## Abstract

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.

**Keywords** Groundwater recharge · Numerical modeling · Unsaturated zone · Soil moisture content · Canada

## 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

<figure><img src="image_3.png" /></figure>

<footnote>Sabrina Bruneau
sabrina.bruneau.f@gmail.com</footnote>

<footnote n="1">¹ Département des sciences de la Terre et de l'atmosphère, Université du Québec à Montréal, C.P. 8888, Succ. Centre-Ville, Montréal, QC H3C 3P8, Canada</footnote>

<footnote n="2">² GEOTOP Research Center - Université du Québec à Montréal, C.P. 8888, Succ. Centre-Ville, Montréal, QC, Canada</footnote>

<footnote n="3">³ GRIL Research Center - Université du Québec à Montréal, C.P. 8888, Succ. Centre-Ville, Montréal, QC, Canada</footnote>

<footnote n="4">⁴ MR 1402 ECOSYS, AgroParisTech, INRAE, Université Paris Saclay, 78850 Thiverval-Grignon, France</footnote>

<footnote n="5">⁵ MINES ParisTech, Centre de Géosciences, PSL Research University, 35 rue Saint Honoré, 77300 Fontainebleau, France</footnote>
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–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 \times 10^6 \text{ 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

<figure><img src="image_4.png" /></figure>
are now characterized as deltaic and littoral glaciomarine
deposits. Depending on the energy regime and the progres-
sion 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 evalu-
ate 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 domi-
nant 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).

**Materials and methods**

Site instrumentation

The model input data are the volume of liquid water availa-
ble 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 cali-
brated and the generic calibration accuracy of 0.03 cm³/cm³
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 temper-
ature. Groundwater levels at the three sites were monitored
hourly by automatic pressure transducers (Solinst Levelog-
gers 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
(deep 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 (Healy and Cook 2002). 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).

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Location of the study sites, the groundwater recharge research infrastructure network (IRRES), and the position of the extracted meteorological datasets in the Vaudreuil-Soulanges region</figcaption></figure>

<figure><img src="image_6.png" /></figure>
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⁻¹ 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 ANUSPLIN 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:

$$ \theta(h) = \begin{cases} \theta_r + \frac{\theta_s - \theta_f}{(1 + |\alpha h|^m)^m}, & h < 0 \\ \theta_s, & h \ge 0 \end{cases} \quad (1) $$

$$ K(S_e) = K_s \times S_e^l \times [1 - (1 - S_e^{1/m})^m]^2 \quad (2) $$

where $\theta$ [L³/L³] is the volumetric moisture content; $h$ [L] is the pressure head; $\theta_r$ and $\theta_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 = (\theta - \theta_r)/(\theta_s - \theta_f)$ [-] is the saturation degree. The fitting factor $\alpha$ [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 ($\theta_r$, $\theta_s$, $K_s$, $\alpha$ 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 field-measurements 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

<table><caption>**Table 1** Recharge monitoring site characteristics (VWC stands for volumetric water content)</caption><thead><tr><td>Station</td><td>Land cover</td><td>Depth of VWC sensors</td><td>Time period</td></tr></thead><tbody><tr><td>SLZA</td><td>Grassland</td><td>10, 20, 50 and 100 cm</td><td>Aug 2015 to Dec 2018</td></tr><tr><td>SLZB</td><td>Jack pine forest</td><td>10, 20, 50 and 100 cm</td><td>Oct 2016 to Dec 2018</td></tr><tr><td>STEL</td><td>Prairie grassland</td><td>25, 50, 75 and 100 cm</td><td>Dec 2015 to Dec 2018</td></tr></tbody></table>
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):

$$ \text{if } T_{\text{max}} \le 0^{\circ}\text{C, SnowFrac} = 1 \quad (3) $$

$$ \text{if } T_{\text{min}} \ge 0^{\circ}\text{C, SnowFrac} = 0 \quad (4) $$

$$ \text{else SnowFrac} = 1 - \frac{T_{\text{max}}}{T_{\text{max}} - T_{\text{min}}} \quad (5) $$

where $T_{\text{min}}$ and $T_{\text{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:

$$ \text{Melt} = \begin{cases} C_{\text{melt}} \times (T_{\text{air}} - T_{\text{melt}}), & T_{\text{air}} > T_{\text{melt}} \\ 0, & T_{\text{air}} < T_{\text{melt}} \end{cases} \quad (6) $$

where Melt is the daily snowmelt (mm/day), $C_{\text{melt}}$ represents the snowmelt rate (mm/°C/day), $T_{\text{air}}$ is the mean daily air temperature (°C), and $T_{\text{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:

$$ E_p(t) = \text{PET}(t) \times e^{-k \times \text{LAI}(t)} \quad (7) $$

$$ T_p(t) = \text{PET}(t) - E_p(t) \quad (8) $$

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:

$$ S(h) = \alpha(h) \times S_p \quad (9) $$

where $\alpha(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

<table><caption>**Table 2** Description of daily meteorological variables used in this study, where *T* stand for air temperature and *P* for precipitation</caption><thead><tr><th>Dataset name</th><th>Source</th><th>Type</th><th>Variable</th><th>Time period</th><th>Coverage</th></tr></thead><tbody><tr><td>STEL</td><td>IRRES station at Saint-Teles-phore site</td><td>In situ weather station</td><td><i>T</i>, <i>P</i>, longwave and shortwave radiation, wind speed</td><td>2016 to 2018</td><td>Local</td></tr><tr><td>ANUSPLIN</td><td>Natural Resources Canada (NRCan)</td><td>Interpolated</td><td><i>T</i><sub>min</sub>, <i>T</i><sub>max</sub>, <i>P</i></td><td>1950 to 2017</td><td>10-km grid across Canada</td></tr><tr><td>NARR</td><td>NOAA's National Center for Atmospheric Prediction (NCEP)</td><td>Reanalysis</td><td><i>T</i><sub>min</sub>, <i>T</i><sub>max</sub>, <i>P</i>, longwave and shortwave radiation, wind speed, albedo</td><td>1979 to 2018</td><td>32-km grid across North America</td></tr></tbody></table>
biomass by 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 negligible.

## 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 ($\theta_r$, $\theta_s$, $\alpha$, $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, $\theta_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 $\theta_r$ than $\theta_s$) which made it difficult to calibrate $\theta_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 $\theta_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 the ESM)

<table><caption>**Table 3** Soil layer characterization at the three study sites from field and laboratory analyses (Field $K_s$ stands for saturated hydraulic conductivity calculated from Guelph permeameter measurements, and Predicted $K_s$ stands for saturated hydraulic conductivity predicted by ROSETTA)</caption><thead><tr><th>IRRES site</th><th>Depth (cm)</th><th>Sand (%)</th><th>Silt (%)</th><th>Organic matter (%)</th><th>Field K<sub>s</sub> (cm/day)</th><th>Predicted K<sub>s</sub> (cm/day)</th></tr></thead><tbody><tr><td rowspan="4">SLZA</td><td>0–15</td><td>92.6</td><td>7.4</td><td>2.7</td><td>—</td><td>375</td></tr><tr><td>16–30</td><td>92.6</td><td>7.4</td><td>0</td><td>1036</td><td>549</td></tr><tr><td>31–70</td><td>96.7</td><td>3.3</td><td>0</td><td>950</td><td>959</td></tr><tr><td>71–300</td><td>100</td><td>0</td><td>0</td><td>—</td><td>950</td></tr><tr><td rowspan="4">SLZB</td><td>0–15</td><td>92.6</td><td>7.4</td><td>6.9</td><td>—</td><td>216</td></tr><tr><td>16–30</td><td>92.6</td><td>7.4</td><td>5.4</td><td>241</td><td>257</td></tr><tr><td>31–60</td><td>96.7</td><td>3.3</td><td>3</td><td>820</td><td>641</td></tr><tr><td>60–300</td><td>100</td><td>0</td><td>0</td><td>—</td><td>950</td></tr><tr><td rowspan="4">STEL</td><td>0–40</td><td>86.8</td><td>13.2</td><td>2.5</td><td>362</td><td>184</td></tr><tr><td>41–65</td><td>94.7</td><td>5.3</td><td>0</td><td>1296</td><td>736</td></tr><tr><td>66–80</td><td>94.7</td><td>5.3</td><td>0</td><td>—</td><td>736</td></tr><tr><td>81–300</td><td>98.8</td><td>1.2</td><td>0</td><td>1900</td><td>1243</td></tr></tbody></table>
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.

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² = 0.99 and
slope = 1 for ANUSPLIN and R² = 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 pre-
cipitation 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 meas-
ured 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 equiva-
lent 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² = 0.82
and slope = 1.01), with lower agreement for the calcula-
tion of net radiation from temperature values (R² = 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 Ks of 362 cm/day obtained with the Guelph permeameter. The grain size shifted to a silt fraction of 5.3% and a Ks 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 Ks 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 sig-
nificant role in hydraulic properties but is not considered in

<figure><img src="image_7.png" /></figure>
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.

## Temporal dynamic of soil moisture contents

### 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 \text{ cm}^3/\text{cm}^{-3}$ (mean $0.15 \text{ cm}^3/\text{cm}^3$), while in the prairie site (SLZA), it varied between 0.02 and $0.19 \text{ cm}^3/\text{cm}^3$, (mean $0.12 \text{ cm}^3/\text{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 \text{ cm}^3/\text{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 \text{ cm}^3/\text{cm}^3$ at SLZB, between 0 and $0.10 \text{ cm}^3/\text{cm}^3$ at SLZA, and between 0.04 and $0.08 \text{ cm}^3/\text{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 ($\theta_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 \text{ cm}^3/\text{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
<figure><img src="image_8.png" /><figcaption>**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</figcaption></figure>

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.

## Estimated recharge

### 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

<table><caption>**Table 4** Simulated water-budget components of the model at the three study sites for the calibration period (mm/y)</caption><thead><tr><th>IRRES site</th><th>Water-budget component</th><th>2016</th><th>2017</th><th>2018</th></tr></thead><tbody><tr><td rowspan="4">SLZA</td><td>GR</td><td>619</td><td>703</td><td>486</td></tr><tr><td>VI</td><td>1034</td><td>1208</td><td>829</td></tr><tr><td>PET</td><td>680</td><td>677</td><td>701</td></tr><tr><td>AET</td><td>387</td><td>460</td><td>368</td></tr><tr><td rowspan="4">SLZB</td><td>GR</td><td>514</td><td>590</td><td>344</td></tr><tr><td>VI</td><td>1034</td><td>1208</td><td>829</td></tr><tr><td>PET</td><td>680</td><td>677</td><td>701</td></tr><tr><td>AET</td><td>493</td><td>620</td><td>500</td></tr><tr><td rowspan="4">STEL</td><td>GR</td><td>446</td><td>575</td><td>374</td></tr><tr><td>VI</td><td>935</td><td>1205</td><td>830</td></tr><tr><td>PET</td><td>684</td><td>681</td><td>701</td></tr><tr><td>AET</td><td>474</td><td>603</td><td>461</td></tr></tbody></table>

GR stands for groundwater recharge, VI for vertical inflows, PET potential evapotranspiration and AET actual evapotranspiration

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 $\theta_s$ and $\theta_r$. The calibrated VGM parameters were thus used
to estimate S<sub>y</sub>, 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<sub>y</sub> 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<sub>y</sub> values and assess the resulting uncertainty in the WTF method, GR was also calculated with a S<sub>y</sub> 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

<figure><img src="image_9.png" /><figcaption>**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</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**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</figcaption></figure>

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 5** Observed and simulated soil volumetric water content (VWC) during calibration at different soil depths at the STEL 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</figcaption></figure>

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.

Discussion

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 moist-
ture 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 char-
acteristics (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 percola-
tion and lower water holding capacity.

The relationship between moisture and soil characteristics
was also demonstrated by the measured VWC dynamic at
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² 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 water-holding 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

<figure><img src="image_13.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_14.png" /></figure>
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

<figure><img src="image_15.png" /><figcaption>**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</figcaption></figure>

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
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 it 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–2009), 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 (Fortin and Turcotte 2007). 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).

<figure><img src="image_16.png" /></figure>
## 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.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10040-021-02429-6.

**Acknowledgements** This work was funded by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant of Florent Barbecot (award number RGPIN-2015-04966) and by a MITACS scholarship. The IRRES network was funded by the Canadian Foundation for Innovation (CFI). The authors would like to thank the private landowners for making their property available for this study.

## Declarations

**Conflict of interest** On behalf of all authors, the corresponding author states that there is no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Allen RG, Pereira LS, Raes D, Smith M (1998) Crop evapotranspiration-guidelines for computing crop water requirements - FAO irrigation and drainage paper 56. FAO, Rome

Allen DM, Cannon AJ, Toews MW, Scibek J (2010) Variability in simulated recharge using different GCMs. Water Resour Res 46(10). https://doi.org/10.1029/2009wr008932

Assefa KA, Woodbury AD (2013) Transient, spatially varied groundwater recharge modeling. Water Resour Res 49(8):4593-4606

Baalousha HM, Barth N, Ramasomanana FH, Ahzi S (2018) Groundwater recharge estimation and its spatial distribution in arid regions using GIS: a case study from Qatar karst aquifer. Modelling Earth Syst Environ 4(4):1319-1329. https://doi.org/10.1007/s40808-018-0503-4

Barbecot F, Guillon S, Pili E, Larocque M, Gibert-Brunet E, Hélie J-F, Noret A, Plain C, Schneider V, Mattei A (2018) Using water stable isotopes in the unsaturated zone to quantify recharge in two contrasted infiltration regimes. Vadose Zone J 17(1), https://doi.org/10.2136/vzj2017.09.0170

Bot A, Benites J (2005) The importance of soil organic matter: key to drought-resistant soil and sustained food production. Food & Agriculture Org, Rome

Boumaiza L, Chesnaux R, Walter J, Stumpp C (2020) Assessing groundwater recharge and transpiration in a humid northern region dominated by snowmelt using vadose-zone depth profiles. Hydrogeol J 28(7):2315-2329. https://doi.org/10.1007/s10040-020-02204-z

Choi W, Kim SJ, Rasmussen PF, Moore AR (2009) Use of the north American regional reanalysis for hydrological modeling in Manitoba. Canadian Water Resour J. 34(1):17-36

Cohen P (1965) Water resources of the Humboldt River valley near Winnemucca. Nevada, US Geological Survey Water-Supply Paper

Crosbie RS, Binning P, Kalma JD (2005) A time series approach to inferring groundwater recharge using the water table fluctuation method. Water Resour Res 41(1), W01008. https://doi.org/10.1029/2004WR003077

Dean WE (1974) Determination of carbonate and organic matter in calcareous sediments and sedimentary rocks by loss on ignition; comparison with other methods. J Sediment Res 44(1):242-248

Dripps WR, Bradbury KR (2007) A simple daily soil-water balance model for estimating the spatial and temporal distribution of groundwater recharge in temperate humid areas. Hydrogeol J 15(3):433-444

Dubois E, Larocque M, Gagné S, Meyzonnat G (2021) Simulation of long-term spatiotemporal variations in regional-scale groundwater recharge: contributions of a water budget approach in southern Quebec. Hydrol Earth Syst Sci [preprint]. https://doi.org/10.5194/hess-2021-71 In review 2021

Feddes RA, Kowalik PJ, Zaradny H (1978) Simulation of field water use and crop yield. Jonh Wiley & Sons, New York

Folk RL (1954) The distinction between grain size and mineral composition in sedimentary-rock nomenclature. J Geol 62(4):344-359

Fortin V, Turcotte R (2007) Le modèle hydrologique MOHYSE (the MOHYSE hydrological model). Research report, Environnement Canada

Guber A, Gish T, Pachepsky Y, van Genuchten MT, Daughtry C, Nicholson T, Cady R (2008) Temporal stability in soil water content patterns across agricultural fields. Catena. 73(1):125-133

Guerrero-Morales J, Fonseca CR, Goméz-Albores MA, Sampedro-Rosas ML, Silva-Gómez SE (2020) Proportional variation of potential groundwater recharge as a result of climate change and land-use: a study case in mexico. Land. 9(10):364

Guo H, Jiao JJ, Weeks EP (2008) Rain induced subsurface airflow and Lisse effect. Water Resour Res 44(7), W07409. https://doi.org/10.1029/2007WR006294

Healy RW (2010) Estimating groundwater recharge. Cambridge University Press, UK

Healy RW, Cook PG (2002) Using groundwater levels to estimate recharge. Hydrogeol J 10(1):91-109

Heiri O, Lotter AF, Lemcke G (2001) Loss on ignition as a method for estimating organic and carbonate content in sediments: reproducibility and comparability of results. J Paleolimnol 25(1):101-110

Hopkinson RF, McKenney DW, Milewska EJ, Hutchinson MF, Papadopol P, Vincent LA (2011) Impact of aligning climatological day on gridding daily maximum-minimum temperature and precipitation over Canada. J Appl Meteorol Climatol. 50(8):1654-1665

Hopmans J W, Simunek J, Romano N, Durner W (2002) Inverse methods. In methods of soil analysis, part 4: physical methods. 5:963-1008

Hu W, Wang YQ, Li HJ, Huang MB, Hou MT, Li Z, She DL, Si BC (2019) Dominant role of climate in determining spatio-temporal distribution of potential groundwater recharge at a regional scale. J Hydrol 578:124042. https://doi.org/10.1016/j.jhydrol.2019.124042

Hung Vu V, Merkel BJ (2019) Estimating groundwater recharge for Hanoi, Vietnam. Sci Total Environ 651:1047-1057. https://doi.org/10.1016/j.scitotenv.2018.09.225

Hutchinson MF, McKenney DW, Lawrence K, Pedlar JH, Hopkinson RF, Milewska E, Papadopol P (2009) Development and testing of Canada-wide interpolated spatial models of daily minimum-maximum temperature and precipitation for 1961-2003. J Appl Meteorol Climatol 48(4):725-741
Jackson TJ (2002) Remote sensing of soil moisture: implications for groundwater recharge. Hydrogeol J 10(1):40–51

Jiménez-Martínez J, Skaggs TH, van Genuchten MT, Candela L (2009) A root zone modeling approach to estimating groundwater recharge from irrigated areas. J Hydrol 367(1–2):138–149

Johnson A I (1967) Specific yield – Compilation of specific yields for carious materiales, U.S. Geological Survey Water-Supply Paper 1662-D

Langlois A, Kohn J, Royer A, Cliche P, Brucker L, Picard G, Fily M, Derksen C, Willemet J (2009) Simulation of snow water equivalent (SWE) using thermodynamic snow models in Québec, Canada. J Hydrometeorol 10(6):1447–1463

Larocque M, Meyzonnat G, Ouellet M A, Graveline M H, Gagné S, Barnetche D, Dorner S (2015) Projet de connaissance des eaux souterraines de la zone de Vaudreuil-Soulanges - rapport scientifique (Groundwater knowledge project in the Vaudreuil-Soulanges area - scientific report). Rapport déposé au ministère du Développement durable, de l'Environnement et de la Lutte contre les Changements Climatiques, Quebec, Canada

Loheide SP, Butler JJ Jr, Gorelick SM (2005) Estimation of groundwater consumption by phreatophytes using diurnal water table fluctuations: a saturated-unsaturated flow assessment. Water Resour Res 41(7), D02114. https://doi.org/10.1029/2005WR003942

Maggioni V, Sapiano MR, Adler RF, Tian Y, Huffman GJ (2014) An error model for uncertainty quantification in high-time-resolution precipitation products. J Hydrometeorol 15(3):1274–1292

Marquardt DW (1963) An algorithm for least-squares estimation of nonlinear parameters. J Soc Ind Appl Math 11(2):431–441

Martínez-Fernández J, Ceballos A (2003) Temporal stability of soil moisture in a large-field experiment in Spain. Soil Sci Soc Am J. 67(6):1647–1656

Mattei A, Barbecot F, Goblet P, Guillon S (2020a) Pore water isotope fingerprints to understand the spatiotemporal groundwater recharge variability in ungauged watersheds. Vadose Zone J 19(1). https://doi.org/10.1002/vzj2.20066

Mattei A, Goblet P, Barbecot F, Guillon S, Coquet Y, Wang S (2020b) Can soil hydraulic parameters be estimated from the stable isotope composition of pore water from a single soil profile? Water. 12(2):393

McKenney DW, Hutchinson MF, Papadopol P, Lawrence K, Pedlar J, Campbell K, Milewska E, Hopkinson RF, Price D, Owen T (2011) Customized spatial climate models for North America. Bull Am Meteorol Soc 92(12):1611–1622

MELCC (Ministère de l'Environnement et de la Lutte contre les changements climatiques) (2019a) Données du Programme de surveillance du climat Québec: Direction générale du suivi de l'état de l'environnement (Data from the Quebec Climate Monitoring Program: General Directorate for Monitoring the State of the Environment), http://www.environnement.gouv.qc.ca/climat/normales/index.asp. Accessed 20 Aug 2019

MELCC (2019b) Programme d'acquisition de connaissance sur les eaux souterraines (Groundwater knowledge acquisition program), http://www.environnement.gouv.qc.ca/eau/souterraines/programes/acquisition-connaissance.htm. Accessed 25 Aug 2019

Mesinger F, DiMego G, Kalnay E, Mitchell K, Shafran PC, Ebisuzaki W, Jović D, Woollen J, Rogers E, Berbery EH (2006) North American regional reanalysis. Bull Am Meteorol Soc. 87(3):343–360

Min L, Shen Y, Pei H (2015) Estimating groundwater recharge using deep vadose zone data under typical irrigated cropland in the piedmont region of the North China Plain. J Hydrol 527:305–315

Mualem Y (1976) A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour Res 12(3):513–522

Myneni R, Knyazikhin Y, Park T (2015) MCD15A3H MODIS/Terra+Aqua Leaf Area Index/FPAR 4-day L4 Global 500m SIN Grid V006 [satellite data], NASA EOSDIS Land Processes DAAC

NASA (2017) Quebec Flooding 2017, NASA Earth Science Disasters Program, https://disasters.nasa.gov/quebec-flooding-2017. Accessed 5 July 2019

Ouranos (2015) Vers l'adaptation: Synthèse des connaissances Sur les changements climatiques au Québec, 2015 (towards adaptation: synthesis of knowledge on climate change in Quebec, 2015) edn., Montreal, Quebec

Ritchie JT (1972) Model for predicting evaporation from a row crop with incomplete cover. Water Resour Res 8(5):1204–1213

Rivard C, Vigneault H, Piggott AR, Larocque M, Anctil F (2009) Groundwater recharge trends in Canada. Canadian J of Earth Sci 46(11):841–854

Rivard C, Lefebvre R, Paradis D (2014) Regional recharge estimation using multiple methods: an application in the Annapolis Valley, Nova Scotia (Canada). Environ Earth Sci 71(3):1389–1408. https://doi.org/10.1007/s12665-013-2545-2

Rudolph T D (1985) Jack pine (Pinus banksiana Lamb), For Serv

Saha GC, Li J, Thring RW, Hirshfield F, Paul SS (2017) Temporal dynamics of groundwater-surface water interaction under the effects of climate change: a case study in the Kiskatinaw River watershed. Canada, Investig of Coast Aquifers 551:440–452. https://doi.org/10.1016/j.jhydrol.2017.06.008

Saïd Ahmed H, Coquet Y (2018) Water uptake by date palm on Haplic Luvisols in the Djibouti coastal plain, Geoderma Regional. 15:e00189. https://doi.org/10.1016/j.geodrs.2018.e00189

Scanlon BR, Healy RW, Cook PG (2002) Choosing appropriate techniques for quantifying groundwater recharge. Hydrogeol J 10(1):18–39

Schaap MG, Leij FJ, van Genuchten MT (2001) Rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J Hydrol 251(3):163–176

Scharnagl B, Vrugt JA, Vereecken H, Herbst M (2011) Inverse modeling of in situ soil water dynamics: investigating the effect of different prior distributions of the soil hydraulic parameters. Hydrol Earth Syst Sci 15(10):3043–3059

Sentelhas PC, Gillespie TJ, Santos EA (2010) Evaluation of FAO pentman-Monteith and alternative methods for estimating reference evapotranspiration with missing data in southern Ontario, Canada. Agric Water Manag 97(5):635–644

Shen Y, Xiong A, Wang Y, Xie P (2010) Performance of high resolution satellite precipitation products over China. J Geophys Res Atmos 115, D02114. https://doi.org/10.1029/2009JD012097

Simunek J, van Genuchten MT, Wendroth O (1998) Parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci Soc Am J. 62(4):894–905

Simunek J, Sejna M, Saito H, Sakai M, van Genuchten M T (2015) The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media, version 4.17. (HYDRUS software series 3), Department of Environmental Sciences, University of California Riverside, USA

Singh J, Knapp HV, Demissie M (2004) Hydrologic modeling of the Iroquois River watershed using HSPF and SWAT (ISWS CR 2004–08). Illinois State Water Survey, Champaign, Illinois

Soil moisture Equipment Corp (2013) Guelph permeameter calculator spreadsheet, http://www.soilmoisture.com/operating.html. Accessed 14 Nov 2017

Szilagyi J, Zlotnik VA, Gates JB, Jozsa J (2011) Mapping mean annual groundwater recharge in the Nebraska Sand Hills, USA. Hydrogeol J 19(8):1503–1513

Thoma MJ, Barrash W, Cardiff M, Bradford J, Mead J (2014) Estimating unsaturated hydraulic functions for coarse sediment from a field-scale infiltration experiment. Vadose Zone J 13(3). https://doi.org/10.2136/vzj2013.05.0096

Turcotte R, Fortin LG, Fortin V, Fortin JP, Villeneuve JP (2007) Operational analysis of the spatial distribution and the temporal
evolution of the snowpack water equivalent in southern Quebec,
Canada. Nordic Hydrol 38(3):211–234

Turgeon F, Larocque M, Meyzonnat G, Dorner S, Bourgault M-A
(2018) Examining the challenges of simulating surface water-
groundwater interactions in a post-glacial environment. Canadian
Water Resour J 43(2):262–280. https://doi.org/10.1080/07011784.
2017.1414635

Turkeltaub T, Kurtzman D, Bel G, Dahan O (2015) Examination of groundwater recharge with a calibrated/validated flow model of the deep vadose zone. J Hydrol 522:618–627

van Genuchten MT (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1. Soil Sci Soc Am J. 44(5):892–898

von Freyberg J, Moeck C, Schirmer M (2015) Estimation of ground
water recharge and drought severity with varying model com-
plexity. J Hydrol 527:844–857. https://doi.org/10.1016/j.jhydrol.
2015.05.025

Vrugt JA, Bouten W, Weerts AH (2001) Information content of data for
identifying soil hydraulic parameters from outflow experiments.
Soil Sci Soc Am J 65(1):19–27

Wang T, Zlotnik VA, Simunek J, Schaap MG (2009a) Using pedotransfer functions in vadose zone models for estimating groundwater recharge in semiarid regions. Water Resour Res 45(4). https://doi.org/10.1029/2008WR006903

Wang T, Wedin D, Zlotnik VA (2009b) Field evidence of a negative correlation between saturated hydraulic conductivity and soil carbon in a sandy soil. Water Resour Res 45(7), W07503. https://doi.org/10.1029/2008WR006865

Wang T, Franz TE, Zlotnik VA (2014) Controls of soil hydraulic char-
acteristics on modeling groundwater recharge under different cli-
matic conditions. J Hydrol 521:470–481

Wang T, Franz TE, Yue W, Szilagyi J, Zlotnik VA, You J, Chen X, Shulski MD, Young A (2016) Feasibility analysis of using inverse modeling for estimating natural ground water recharge from a large scale soil moisture monitoring network. J Hydrol 533:250–265

Weeks EP (2002) The Lisse effect revisited. Groundw. 40(6):652–656

Wong JS, Razavi S, Bonsal BR, Wheater HS, Asong ZE (2017) Inter-comparison of daily precipitation products for large-scale hydro-climatic applications over Canada. Hydrol Earth Syst Sci 21(4):2163–2185

Xie Y, Cook PG, Simmons CT, Partington D, Crosbie R, Batelaan
O (2018) Uncertainty of groundwater recharge estimated from a
water and energy balance model. J Hydrol 561:1081–1093

Zha T, Barr AG, van der Kamp G, Black TA, McCaughey JH, Flanagan LB (2010) Interannual variation of evapotranspiration from forest and grassland ecosystems in western Canada in relation to drought. Agric For Meteorol 150(11):1476–1484

Zhu Y, Ren L, Horton R, Lü H, Chen X, Jia Y, Wang Z, Sudickt E (2013) Estimating the contribution of groundwater to rootzone soil moisture. Hydrol Res 44(6):1102–1113

**Publisher's note** Springer Nature remains neutral with regard to
jurisdictional claims in published maps and institutional affiliations.