Pore water isotope fingerprints to understand the spatiotemporal groundwater recharge variability in ungauged watersheds

Reliable groundwater recharge quantification at the regional scale (e.g., watershed or subwatershed) is fundamental to sustainable water resource management. Although modeling at the watershed scale is gaining wide support, the long‐term monitoring needed for model calibration is often not readily available, as many watersheds worldwide remain ungauged. In response to this situation, we propose a new approach to estimate groundwater recharge at the watershed scale. This approach is fast and accurate and takes into account the existing variability without requiring long‐term monitoring. Only a single field campaign to acquire soil water content and pore water isotopic composition depth profiles is needed. The principle is to extend a physically based, one‐dimensional unsaturated zone flow model from the local (i.e., profile) to the watershed scale, using an index method for distributed recharge based on a GIS. The methodology was validated in a gauged watershed, where previous studies have estimated recharge using a spatialized water balance model calibrated using long‐term discharge monitoring data. Scaling was investigated by comparing recharge values obtained using the local‐scale approach at 10 study sites within the watershed with coinciding values obtained at the watershed scale. Recharge values were similar in terms of both dynamics and quantity. Using the pore water isotopic fingerprint of ungauged watersheds is therefore confirmed to be a suitable approach for understanding spatiotemporal recharge variability.


INTRODUCTION
Groundwater recharge is usually considered to be the downward movement of water below the maximum rooting depth through the unsaturated zone (Simmers, 1988). The quantification of aquifer recharge is fundamental

Core Ideas
• Coupling 1D unsaturated zone model with GISbased method allows watershed-scale recharge estimation. • Pore water isotopes were used to determine soil hydraulic parameters without long-term monitoring. • This approach improves the understanding of time-space variability of recharge in ungauged watershed.
At the profile scale, and based on a single field campaign, pore water stable isotopes have been shown to be valuable tools to obtain insights into vertical water flow under various climatic, geological, and land cover conditions (Barbecot et al., 2018;Sprenger, Seeger, Blume, & Weiler, 2016;Stumpp, Maloszewski, Stichler, & Maciejewski, 2007;Stumpp, Stichler, Kandolf, & Šimůnek, 2012). Easy to sample in the field, and fast to measure in the laboratory, thanks to the development of H 2 O liquid -H 2 O vapor equilibration methods (Mattei et al., 2019;Wassenaar, Hendry, Chostner, & Lis, 2008), pore water isotopes are also relevant to calibrating unsaturated zone models (Sprenger, Volkmann, Blume, & Weiler, 2015), even without temporal monitoring (Mattei et al., 2020). In ungauged watersheds, pore water stable isotope fingerprints can provide valuable information to estimate soil hydraulic parameters accurately and quickly, and thus recharge at a lower cost when meteorological conditions are known.
Here, we propose to extend the application of physically based, 1D unsaturated zone flow modeling from the profile scale to the watershed scale using an index method for distributed recharge based on a GIS. A new methodology based on the acquisition of both water content and pore water stable isotope composition along vertical soil profiles was developed to determine soil hydraulic parameters using inverse modeling at large scales. We applied this methodology to the Raquette River watershed (150 km 2 ), located in southern Quebec, Canada. To validate this approach, the resulting recharge map was compared with that obtained in a previous study, using a spatialized water balance model calibrated using measured discharge (Larocque et al., 2015). Also, comparisons of recharge values obtained using the local approach (i.e., at the profile scale) and the previous watershed-scale approach were performed in terms of both quantity and dynamics in order to investigate issues linked to the transfer of information across scales, known as scaling (Bierkens, Finke, & Willigen, 2001) The proposed methodology allows the potential natural diffuse recharge, defined as the amount of effective rainfall that percolates and infiltrates under the root zone in the vadose zone, to be estimated. Indeed, flow processes were investigated only in the first meter of depth of the vadose zone, and not down to the water table. As Besbes and De Marsily (1984) demonstrated, actual and potential recharge are identical on average and can be distinguished only by the transit time in the unsaturated zone, except if infiltrated water is blocked or diverted by semi-impervious materials and drains (de Vries & Simmers, 2002); thereafter, for the Raquette River watershed, we did not distinguish between potential and actual recharge. Also, as in Quebec, O and H isotope values showed a strong linear correlation; only δ 2 H values are described and discussed here.

Study area and study site descriptions
The Raquette River watershed is located 50 km west of Montreal, in the Vaudreuil-Soulanges area (Quebec, Canada; Figure 1). It contains eight municipalities and has a population of > 76,000 inhabitants over an area of 133 km 2 . The current water supply is principally groundwater from public and private wells (Larocque et al., 2015).
The Raquette River watershed is within the St. Lawrence Lowlands geological region. Like the rest of Quebec, the geomorphology of the study area is marked by glaciationdeglaciation phases, with unconsolidated sediments of glacial and postglacial origin overlying the fractured bedrock to a thickness ranging from 5 to 86 m (Larocque et al., 2015).These unconsolidated sediments largely control the hydrogeological response of the underlying fractured bedrock aquifers. A detailed description of the sediment stratigraphy is available in Larocque et al. (2015).
The studied soils in the Raquette River watershed are classified by texture (gravel, sand loam, till, and clay) and by drainage (from well to poorly drained) (Larocque et al., 2015). Clayey soils are dominant (63%) and are found everywhere in the plain ( Figure 2). Sandy soils also represent a major coverage and are present on the Hudson and Saint-Lazare hills, as well as in the area of Sainte-Justine-de-Newton. Till deposits are mostly on the flanks of Mount Rigaud. Along the Raquette River, bedrock outcrops in only a few locations (< 1%).
The majority of the watershed consists of agricultural lands (51%), located in the clay plains ( Figure 3). Cereal agriculture dominates (75%) (Larocque et al., 2015). Urban areas represent < 5% of the watershed and are also mostly located in the clay plains. As soon as the topography increases, forested area dominates (43%). Forest areas are mainly hardwoods (72%) composed of maples (Acer spp.) and birches (Betula spp.) (Larocque et al., 2015). Mixed forests occupy 23% of the forest area, whereas conifers make up 5% of the forest cover. Less than 3% of the study area is occupied by grassland. There are very few wetlands in the watershed.

Sampling and analysis
Eleven vertical soil profiles were sampled in the study area between 12 Oct. 2018 and 1 Nov. 2018, one at each of the Vadose Zone Journal F I G U R E 3 Land cover class (top) and infiltration capacity (bottom) distributions in the Raquette River watershed sites described above, for water content and pore water isotope composition measurements. Soil samples were taken with a spatula, from between 5 cm from the soil surface to 30-cm depth, and then from 10-to 200-cm maximum depth (see Table 1 for details). Soil samples were stored in airtight glass jars until stable isotope analysis in the laboratory by direct equilibration (Wassenaar et al., 2008) with an improved setup (Mattei et al., 2019) on a triple isotopic water analyzer (Los Gatos Research). Gravimetric water content was measured on the unconsolidated samples and converted into volumetric water content, assuming a bulk density of 1,500 kg m −3 . From January 2015 to December 2018, the isotopic composition of precipitation was monitored at Mont-Saint Bruno (75 km from the study area). Rain was sampled twice a month with a rain gauge consisting of a funnel (16-cm diam.) connected to a plastic bottle (5 L) with a permanent connection to the outside air via a 15-m-long soft polypropylene plastic tube (5-mm i.d.) (Gröning et al., 2012). Snow was collected on a monthly basis using a   Based on the long-term performance of the triple isotopic water analyzer, measurement precision for δ 2 H was 0.2‰ for liquid water samples and 0.3‰ for pore water samples.
Daily weather data (minimum and maximum temperature, precipitation, and relative humidity) were obtained from North American Regional Reanalysis (NARR), with a 32-km grid spacing (Mesinger et al., 2006). Since the distance between the different study sites is < 32 km, the atmospheric conditions are considered to be similar for all of them.

Soil flow and transport model
Water flow and δ 2 H transport were modeled using the flow and mass transport code, METIS (1D), developed by the Geosciences Department of MINES ParisTech (Goblet, 2010). A detailed description of the model setup and parameterization is available in Mattei et al., 2020. The sections below will provide a brief overview.

Water flow
Transient water flow within the unsaturated soil profile was simulated by numerically solving the Richards equation using the finite-element code, METIS (1D), developed by the Geosciences Department of MINES ParisTech (Goblet, 2010): where K s is the saturated hydraulic conductivity [L T −1 ], K r is the relative hydraulic conductivity [-], H is the hydraulic head [L], F is the porosity [L 3 L −3 ], C(h) is the specific retention capacity specified by S the saturation ratio [L 3 L −3 ], and h is the matric head [L]: The behavior of the unsaturated medium depends on two essential characteristics of the material: the relative permeability curve and the suction curve. These relationships are determined by the residual and maximal saturation ratio (S min [L 3 L −3 ] and S max [L 3 L −3 ], respectively), with three empirical parameters defining the shape of the retention curves: α [L −1 ], n [−], and m [-] (with m = 1 − 1/n), following the Mualem-van Genuchten model (Mualem, 1976;van Genuchten, 1980): In order to account for root water uptake, a sink term was defined according to Feddes, Kowalik, Kolinska-Malinka, and Zaradny (1976). The required matric potential that describes the minimal pressure head for the root water uptake was set to −1,000 cm. Root water uptake was parameterized by the site-specific rooting depth and assuming a uniform root distribution.

δ 2 H transport
The δ 2 H transport through the soil profile was simulated according to a modified version of the widely used advection-dispersion model (assuming a single porous medium), in order to account for isotopic enrichment due to fractionation during evaporation (Mattei et al., 2020). Dispersivity, λ [L], is the only parameter describing the advection-dispersion equation.

Model setup
The upper boundary condition for water flow was defined as a time-variable flux boundary conditions by atmospheric data at a daily time step. Air temperature measurements were used to estimate the potential evapotranspiration via the Hargreaves and Samani (1982) formula. The leaf area index (LAI) obtained from MODIS (Myneni, Knyazikhin, & Park, 2015) at a 500-m grid spacing (MCD15A3H) and a canopy radiation extinction factor of 0.463 were used to divide potential evapotranspiration into potential evaporation and potential transpiration according to Sutanto, Wenninger, Coenders-Gerrits, and Uhlenbrook (2012), since evaporation modifies the isotopic composition of infiltrated water, while transpiration does not (Souchez, Lorrain, & Tison, 2002). Snow accumulation and snowmelt were accounted for using the two-parameter semidistributed Snow Accounting Routine model, CemaNeige (Valéry, 2010), calibrated on snow depth data from Coteau-du-Lac (20 km from the study site and supplied by Ministère de l'Environnement et Lutte contre les changements climatiques Québec). A degree day snow melting constant of 0.5 cm d −1 C •−1 and a cold content factor of 0.1 were determined through this calibration. A zero-potential condition was assumed at the bottom of the model. The upper boundary conditions for δ 2 H transport were defined as a time-variable solute flux condition by prescribed concentration data at the surface at a monthly time step. Initial conditions were defined through a steady-state simulation with a given surface water isotope composition, representing the weighted average precipitation composition (−85‰ for δ 2 H). Absolute δ 2 H values were used for the calculations. The modeled period extends from 1 Jan. 2016 to 31 Dec. 2018. The influence of the initial conditions on the calibration are minimized, as a spin-up period of almost 2 yr was simulated (January 2016 to December 2017) prior to the start of the calibration period. For all study sites, the depth of the soil profiles was set to 600 cm and discretized into 600 elements of 1-cm length (i.e., 601 nodes). Each profile was divided into two horizons according to the field observations (Horizon 1 and Horizon 2, Table 1).

Model parameterization
The model parameters were determined through inverse modeling with a global optimization procedure using soil moisture and pore water isotope composition profiles from a single sampling time simultaneously as a fitting target. Indeed, according to Mattei et al., 2020, when modeling is used in order to estimate the annual groundwater recharge, it is possible to calibrate the model without continuous monitoring data if using water content and pore water isotopic composition profiles. The best combination of parameters S min , S max , α, n, K s , λ, and F (porosity) was searched within the boundaries of the parameter space, determined by expert knowledge (Sprenger et al., 2015). Each soil horizon was represented by its individual parameter combination, but to reduce the number of parameters, a single value of λ was calibrated over the entire profile. Even though an F value was determined for each horizon, we assumed a unique value of bulk density for all horizons, set to 1,500 kg m −3 . Thirteen parameters (i.e. α 1 , n 1 , S min1 , S max1 , K r1 , F 1 , α 2 , n 2 , S min2 , S max2 , K r2 , F 2 , and λ) were therefore calibrated for each study site to describe their soil hydraulic and transport properties, according to Mattei et al., 2020. The objective function for the inverse modeling procedure was based on the Kling-Gupta efficiency (KGE) (Gupta, Kling, Yilmaz, & Martinez, 2009). This dimensionless index compares simulated and observed data with regards to their correlation (r), the ratio of their mean values (bias ratio, β), and the ratio of their coefficient of variation (variability ratio, γ) as follows: For each site, in order to simultaneously take both the soil moisture profile (KGE θ ) and the δ 2 H profile (KGE c ) into account, a multi-objective function (KGE tot ) was defined as the average of KGE θ and KGE c . For parameter combinations that did not lead to a numerical convergence of the METIS code, a low value of the objective function was assigned (−100). The set of parameters leading to the best KGE tot was retained.

Potential groundwater recharge estimation
In order to investigate the effects of scaling, potential recharge was evaluated at two scales: the local (i.e., profile) and the large (i.e., watershed) scales. For both, potential groundwater recharge is defined here as the downward water flux [L T −1 ] at 100-cm depth.

Local-scale potential recharge determination
The application of physically based, 1D unsaturated zone flow models, such as the model described above, is classically performed to determine potential groundwater recharge locally (Mattei et al., 2020;Sprenger et al., 2015;Stumpp, Stichler, Kandolf, & Šimůnek, 2012). In such a local-scale approach, all physical and geographical data required by the model (LAI, root zone depth, and horizons depth) are derived from site observations and are site dependent. Soil hydrodynamic parameters are also site specific and were determined using the measured depth profiles of both pore water isotope composition and water content through the inverse modeling approach described above.

Large-scale potential recharge determination
To expand potential recharge estimation from the local to the large scale (e.g., watershed proposed here), an index method based on a GIS was used to identify soil map units, which are areas homogeneous in terms of soil properties and effective infiltration. In each of these units, a physically based, 1D unsaturated zone flow model was developed to estimate potential recharge. First, following Galvão, Hirata, and Conicelli (2018), the soil types and land cover maps were reclassified to permeability and capacity for infiltration maps, respectively. A slope map was generated and used to derive the effective infiltration capacity. Second, hydrodynamic parameters were determined for each permeability class (i.e., soil profiles belonging to a given permeability class were aggregated). Soil hydraulic parameters obtained using the largescale approach are independent of those of the soil profiles of the same permeability class obtained using the local-scale approach. Third, modeling was performed at the scale of the identified spatial units, using the capacity of infiltration class and permeability class parameters. With scaling, the assumption of zero surface runoff made at the local scale, because no runoff was observed in the field, is no longer valid at the large scale. The amount of precipitation that flows on the surface was therefore determined, for each soil map unit, using permeability, capacity for infiltration, and slope maps, following the methodology proposed by the American Society of Civil Engineers (1969). Finally, a potential recharge map was generated. These methods will be described in further detail below.

Soil map unit definition
The study area was discretized into 25-mš × 25-mš cells, corresponding to the soil map units. Following Galvão et al. (2018), map units were defined as a function of three criteria: • Permeability. The permeability of a soil represents the ability of water to move through it. It is primarily controlled by soil texture and soil organic matter content (Leeper & Uren, 1993). Usually, the finer the soil texture, the lower the permeability. The permeability map was obtained by reclassifying the soil type map of the Raquette River watershed (Figure 2) into three different types: low (clayey soils and bedrock outcrop), moderate (silty soils and till), and high permeability (sandy soils and gravelly soils). • Infiltration capacity. Soil infiltration capacity is the maximum rate at which a soil in a given condition will absorb water. Several studies have pointed out the link between land cover and use and infiltration capacity (Benavides et al., 2018;Cameron, Shaykewich, DeJong, Chanasyk, & Green, 1981;Fu, Chen, Ma, Zhou, & Wang, 2000;Tromble, Renard, & Thatcher, 1974;Wood & Blackburn, 1981). Indeed, vegetation can slow runoff, allowing more time for water to seep into the ground (Allen, 1997). Impervious surfaces such as urban areas act as a shortcut for rainfall, draining the water directly into streams (Fetter, 1994). Agriculture and tillage also change the Vadose Zone Journal infiltration patterns of a landscape. Water that would infiltrate directly into the soil under natural conditions now runs into streams. An infiltration capacity map was obtained by reclassifying the land cover map of the Raquette River watershed (Figure 3) into five major types: cultivated, urban, forest, grassland (including all of bare soil, wild land, quarry), wetland, with three subtypes for forest (hardwood, mixed, and softwood). • Slope. The slope map was generated using elevation data extracted from the digital elevation model (Figure 1). Three slope intervals were defined: < 2% (low slopes), 2-7% (moderate slopes), and > 7% (high slopes).
Combining permeability, infiltration capacity, and slope maps but excluding wetlands, where recharge processes are different, potential groundwater recharge estimates were carried out for 54 map unit types (Table 2).

Runoff coefficient determination
The runoff coefficient, C, is a dimensionless ratio that represents the fraction of rainfall converted to runoff. In practice, C can be calculated from tables found in the literature. One of the most used is that proposed by the American Society of Civil Engineers (1969) ( Table 3). The C values for each type of map unit were determined according to the equation = 1 − ( ′ 1 + ′ 2 + ′ 3 ), and C 1 , C 2 , and C 3 values presented in Table 3. A runoff coefficient was thus calculated for each of the 54 map units considered (Table 2). We assumed that the amount of effective precipitation that runs off does not infiltrate into the soil (never, nowhere in the study area). This is equivalent to removing this amount of water from the water input to the model.

Potential recharge computation
To compute potential recharge, soil hydrodynamic parameters, root zone depth, soil horizon depths, and LAI are required.
At the soil map unit scale, root zone depth, first soil horizon depth, and LAI were estimated for each infiltration capacity class by averaging observed data from the field. For urban areas, we did not have any field observations. However, since urban areas consist mostly of lawns in the study area, we assumed that the soil horizon depths and root zone depth were identical to those for grassland. The LAI was estimated following Sutanto et al. (2012), assuming a crop height of 0.10 m when snow is totally melted. For each infiltration capacity class, Table 4 summarizes the root zone depth, first soil horizon depth, and LAI range applied to compute potential groundwater recharge.
Soil hydrodynamic parameters were determined for each permeability class through inverse modeling by aggregating, respectively, soil water content and δ 2 H depth profiles observed over the study area by permeability class. This aggregation approach consists of grouping several soil profiles within the same objective function. For the high-permeability class, seven soil profiles were aggregated (1-7); for the moderate-permeability class, three (8-10) were aggregated; and as only one depth profile was made in the low-permeability class area (11), the soil hydrodynamic parameters were directly estimated from it (see Table 1 for details). Map unit hydrodynamic parameters were attributed as a function of their permeability class. Aggregation of the data by permeability class was chosen because hydraulic parameters are directly linked with soil textural classes through pedotransfer functions (Vereecken et al., 2010). However, Gribb, Forkutsa, Hansen, Chandler, and McNamara (2009) found that hydraulic properties estimated by inversion from in situ measurements provided better descriptions of the soil water dynamics than pedotransfer functions. We therefore estimated soil hydrodynamic parameters by inverse modeling and not using pedotransfer functions. The methodology used for inverse modeling is the same as that described in the model parameterization section.
Finally, daily potential recharge was calculated for each map unit using METIS according to the parameters thereof (soil hydrodynamic properties, root depth, first horizon depth, LAI range, and runoff coefficient), and an annual cumulative potential recharge map was generated by taking the sum of the calculated daily values.

Local-scale approach
Annual cumulative potential recharge estimation was carried out at 11 locations across the study area through inverse modeling using a calibration based on both pore water content and δ 2 H depth profiles (see Supplemental  Table S1 for values). Water content profiles (KGE θ values ranging from 0.5 to 0.9) and δ 2 H profiles (KGE c values ranging from 0.4 to 0.8) were satisfactorily reproduced ( Table 5). The KGE tot values obtained ranged from 0.5 to 0.8. The poorer fit of the δ 2 H profiles could be explained, on the one hand, by the fact that precipitation samples for isotope measurements were not taken from the study site itself, but at a distance of almost 75 km, and, on the other hand, by the fact that sampling was carried out monthly. Field δ 2 H values were therefore representative of a monthly mean value. It was therefore not possible to reproduce δ 2 H values of the first centimeters of the soil  profile, which correspond to the last precipitation episode. Below that, the isotopic composition is the result of the mixing of several precipitation events, and thus a monthly precipitation sampling would be representative of the measured value in vadose zone. Simulated annual cumulative potential recharge values ranged from 281 to 487 mm in 2017 and from 100 to 304 mm in 2018 (Table 5). As expected, the mean potential recharge values obtained for the high-permeability class (Study Sites 1-7) are the highest: 403 mm for 2017 and 216 mm for 2018. Annual cumulative potential recharge values obtained for the moderate-permeability (Study Sites 8-10) and low-permeability (Study Site 11) classes are sim-Vadose Zone Journal F I G U R E 4 Terrain slope map, with slopes of <2%, between 2 and 7%, and >7% indicated ilar: 358 and 361 mm for 2017 and 210 and 213 mm for 2018, respectively. However, as there is only one soil profile belonging to the low-permeability class, caution must be taken when interpreting the results.

3.2
Large-scale approach

Runoff map
First, the permeability map (Figure 2, bottom) was generated based on the soil texture map (Figure 2, top). Most of the study area (64%) is covered by low-permeability soils (predominantly clayey soils) found throughout the plain. High-permeability soils (predominantly sandy and gravelly soils) also cover a major part (28%) of the area and are present on the Hudson and Saint-Lazare hills, as well as in the area of Sainte-Justine-de-Newton. Moderatepermeability soils (predominantly silty soils and till) are mainly found on the flanks of Mount Rigaud (8%). Next, the land cover map (Figure 3, top) was reclassified based on the associated capacities for infiltration to obtain the infiltration capacity map (Figure 3, bottom). The distribution of potential high-permeability areas correspond to forested areas, which cover 42% of the study area, whereas low-permeability areas correspond to anthropic and cultural land covers, which span 51 and 4% of the study area, respectively (Figure 2). Less than 3% of the study area is occupied by grassland.
Finally, the surface runoff map was obtained by integrating permeability, infiltration capacity, and terrain slope maps ( Figure 5). For the 54 map units considered, the runoff coefficient, C, ranges from 0.1 to 0.7, depending on the combination of soils features, land cover, and terrain slope (Table 2, Column 7).

Soil hydrodynamic parameters
According to Table 1, Soil Profiles 1-7 were grouped in a single objective function to determine the soil hydrodynamic parameters relevant for the high-permeability context through inverse modeling. Soil Profiles 8-10 were similarly used for the moderate-permeability context. Soil Profile 11 was the only profile used to determine the soil hydrodynamic parameters relevant for the lowpermeability context in this way. For all three permeability contexts, it was possible to identify a unique set of parameters that satisfactorily reproduces all aggregated soil profiles, both in terms of water content and δ 2 H composition. Indeed, values obtained for KGE tot range from 0.61 to 0.77, with values for KGE θ ranging from 0.5 to 0.8, and values for KGE c ranging from 0.72 to 0.75 (Table 6). Contrary to what was expected, the soil saturated permeability values obtained for the three permeability classes were quite similar ( Table 6). The presence of substantial organic matter in the low-permeability profile could explain the high fitted hydraulic conductivity value obtained there. Many reports have shown that high soil organic matter content results in increased water-holding capacity, porosity, infiltration capacity, and hydraulic conductivity (Haynes & Naidu, 1998). However, values for the Mualem-van Genuchten parameters α and n obtained for the three permeability classes differ greatly (Table 6). These two parameters were identified as the being most influential parameters in the sensitivity analysis previously carried out for the model (Mattei et al., 2020). Only a slight modification of these two parameters could lead to very dif-ferent infiltration profiles. Moreover, the sensitivity analysis highlighted the high degree of interaction between the parameters. It is therefore possible that this lack of contrast between the soil permeability values could be due only to interactions between parameters. Consequently, the values obtained for individual parameters are not relevant for calibrating another model (e.g., the soil hydraulic conductivity values are not useful to calibrate a hydrogeological model). Soil hydraulic parameters obtained through the local-scale parameterization were similar to those obtained through the large-scale approach, and permeability values are also quite similar (Supplemental Table S1). This reveals that the model mostly characterized a superficial horizon, the properties of which are less contrasted than those of the deeper horizons.

Potential recharge values
Results of annual potential recharge values for 2017 and 2018 are summarized in Table 2  Annual cumulated potential recharge maps are presented in Figure 6 for 2017 and 2018. The mean annual cumulative potential recharge for the Raquette River watershed is 163 mm for 2017 and 115 mm for 2018, which represents 18 and 20% of the water input (i.e., precipitation), respectively. Maximum annual cumulative potential recharge values for the study area are obtained for Saint-Lazare and Hudson Hill, Mont-Rigaud, and Saint-Justine-de-Newton. These results are in agreement with a previous study aimed at determining preferential recharge areas and the amount recharge of the Vaudreuil-Soulange zone, using a spatial water budget calibrated on baseflow (HydroBudget; Larocque et al., 2015). Indeed, Saint-Lazare and Hudson Hill, Mont-Rigaud, and Saint-Justinede-Newton were also identified as preferential recharge areas. Mean recharge for the period of 1989-2009 was estimated at 183 mm for the Raquette River watershed, which is 19% of the water input. The main difference observed between the two studies is the amount of recharge estimated for the cultivated clay plains. In this area, the presence of agricultural drains explains the difference in potential recharge and recharge mean values.
In our study, we estimated the potential recharge values, whereas Larocque et al. (2015) estimated actual recharge values.

Scaling issues
In order to investigate issues linked to the transfer of information across scales, we compared potential recharge values obtained using local-and large-scale approaches in the study sites, both in terms of dynamics and quantity. As Profile 11 was the only observation data used to determine potential recharge through the large-scale approach for the corresponding map unit, the comparison of local-and large-scale approaches was not performed for this study site. Figure 7 presents the cumulative potential recharge obtained using local-and large-scale approaches for the year 2018. In order to compare the amount of potential recharge estimated through the two approaches, the upper boundary condition was modified in the local-scale approach to consider runoff. The results obtained using the large-scale approach fit well those obtained using a localscale approach. Indeed, the two major recharge events that occurred in the spring and autumn were able to be correctly identified. The KGE values obtained for study sites belonging to high-permeability class map units ranged from 0.56 to 0.92. For study sites belonging to mediumpermeability class map units, KGE values ranged from 0.43 to 0.62. However, comparing cumulative potential recharge curves, we observed an overestimation of evapotranspiration using the large-scale approach at Study Sites 1-3, whereas it is underestimated at Study Sites 8-11. The averaging of vegetation properties (root zone depth and LAI) performed by infiltration capacity class and used in the large-scale approach could explain such results. The large-scale approach could be improved by considering, for example, root zone depth by both infiltration capacity class and permeability class, since root zone depth is also influenced by soil texture class (Kim & Jackson, 2012). However, more soil depth profiles would be needed. Figure 8 compares annual potential cumulative recharge estimated through the local-and large-scale approaches for the 10 high-and moderate-permeability study sites for both 2017 and 2018. Potential recharge is shown as a percentage of water input. Amounts estimated are very similar (R 2 = .87, y = 1.00x + 1.89). The maximum difference observed between estimates at the two scales is <10% (Figure 8). Large-scale estimation can therefore be considered reliable.

CONCLUSION
We propose a new methodology to allow the application of physically based, 1D unsaturated zone flow modeling to be extended from the profile scale to the watershed scale using an index method for distributed potential recharge based on GIS. This methodology, based on a single field campaign to acquire vertical soil depth profiles of both water content and pore water stable isotope composition, provided valuable estimations of potential annual cumulative recharge. These estimations compare well with results from a spatial water balance model calibrated using longterm discharge monitoring. Moreover, annual cumulative potential recharge values obtained using the local-and large-scale approaches for 10 study sites were similar in terms of both dynamics and quantity. A better description of root zone processes, and more generally vegetation characteristics, at the large scale could further improve the already successful results. Pore water isotopic fingerprints of previously ungauged watersheds can therefore be considered suitable for understanding the spatiotemporal variability of natural diffuse potential recharge. As the 1D model approach does not consider subsurface lateral flows, it would be interesting to test the calibration approach developed here in a three-dimensional deterministic hydrogeological model to improve accuracy of potential recharge estimations.

A C K N O W L E D G M E N T S
The authors would like to thank V. Horoi for his help during fieldwork and S. Bruneau for granulometric analyses. We also wish to thank the two anonymous reviewers and the associate editor for their helpful comments. This research was funded by a Natural Sciences and Engineering Research Council (NSERC) discovery grant to F. Barbecot.