# Impacts of changes in groundwater recharge on the isotopic composition and geochemistry of seasonally ice-covered lakes: insights for sustainable management

Marie Arnoux, Florent Barbecot, Elisabeth Gibert-Brunet, John Gibson, Aurélie Noret

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

## ► To cite this version:

Marie Arnoux, Florent Barbecot, Elisabeth Gibert-Brunet, John Gibson, Aurélie Noret. Impacts of changes in groundwater recharge on the isotopic composition and geochemistry of seasonally ice-covered lakes: insights for sustainable management. Hydrology and Earth System Sciences, 2017, 21 (11), pp.5875-5889. 10.5194/hess-21-5875-2017. hal-04236290

HAL Id: hal-04236290

https://cnrs.hal.science/hal-04236290v1

Submitted on 12 Oct 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 CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Impacts of changes in groundwater recharge on the isotopic composition and geochemistry of seasonally ice-covered lakes: insights for sustainable management

Marie Arnoux<fnref n="1" />,<fnref n="2" />, Florent Barbecot<fnref n="1" />, Elisabeth Gibert-Brunet<fnref n="2" />, John Gibson<fnref n="3" />, and Aurélie Noret<fnref n="2" />

*Correspondence to:* Marie Arnoux (marie.arnoux@protonmail.com), Florent Barbecot (barbecot.florent@uqam.ca), Elisabeth Gibert-Brunet (elisabeth.gibert@u-psud.fr), John Gibson (jjgibson@uvic.ca), and Aurélie Noret (aurelie.noret@u-psud.fr)

Received: 29 March 2017 – Discussion started: 22 May 2017

Revised: 8 September 2017 – Accepted: 14 September 2017 – Published: 27 November 2017

**Abstract.** Lakes are under increasing pressure due to widespread anthropogenic impacts related to rapid development and population growth. Accordingly, many lakes are currently undergoing a systematic decline in water quality. Recent studies have highlighted that global warming and the subsequent changes in water use may further exacerbate eutrophication in lakes. Lake evolution depends strongly on hydrologic balance, and therefore on groundwater connectivity. Groundwater also influences the sensitivity of lacustrine ecosystems to climate and environmental changes, and governs their resilience. Improved characterization of groundwater exchange with lakes is needed today for lake preservation, lake restoration, and sustainable management of lake water quality into the future. In this context, the aim of the present paper is to determine if the future evolution of the climate, the population, and the recharge could modify the geochemistry of lakes (mainly isotopic signature and quality via phosphorous load) and if the isotopic monitoring of lakes could be an efficient tool to highlight the variability of the water budget and quality.

Small groundwater-connected lakes were chosen to simulate changes in water balance and water quality expected under future climate change scenarios, namely representative concentration pathways (RCPs) 4.5 and 8.5. Contemporary baseline conditions, including isotope mass balance and geochemical characteristics, were determined through an intensive field-based research program prior to the simulations. Results highlight that future lake geochemistry and isotopic composition trends will depend on four main parameters: lo-

cation (and therefore climate conditions), lake catchment size (which impacts the intensity of the flux change), lake volume (which impacts the range of variation), and lake G index (i.e., the percentage of groundwater that makes up total lake inflows), the latter being the dominant control on water balance conditions, as revealed by the sensitivity of lake isotopic composition. Based on these model simulations, stable isotopes appear to be especially useful for detecting changes in recharge to lakes with a G index of between 50 and 80%, but response is non-linear. Simulated monthly trends reveal that evolution of annual lake isotopic composition can be dampened by opposing monthly recharge fluctuations. It is also shown that changes in water quality in groundwater-connected lakes depend significantly on lake location and on the intensity of recharge change.

## 1 Introduction

For decades, climate change, combined with rapidly expanding urban, industrial, and agricultural water needs, has placed increasing stress on water resources and on groundwater resources in particular. Future pressure on these resources is likely to be even more pronounced, as groundwater is likely to be increasingly exploited to enhance water supply and to alleviate the worsening drought situation in some arid regions (Dragoni and Sukhija, 2008). Many studies have suggested that sustainable groundwater use has to be based

<footnote n="1">¹GEOTOP, Université du Québec à Montréal, Montréal, Québec, H3C 3P8, Canada</footnote>

<footnote n="2">²GEOPS, UMR 8148, CRNS-Université Paris Saclay/Paris-Sud, Orsay, France</footnote>

<footnote n="3">³Alberta Innovates Technology Futures, 3-4476 Markham Street, Victoria, BC V8Z 7X8, Canada</footnote>
on, among other things, a reliable assessment of recharge, which largely controls its evolution. Aquifer recharge refers to the quantity of water reaching the saturated zone of an aquifer, and therefore replenishing the water table. Unfortunately, in many parts of the world, recharge rates are often not well known at the regional scale (Rivard et al., 2013). While aquifer recharge is crucial to supporting sustainable management of regional groundwater resources, it is difficult to accurately estimate, owing mainly to limited data availability, as well as limitations inherent to estimation methods and field measurements (Rivard et al., 2013). Recharge rates are controlled by geology, soil characteristics, topography, land cover, land use, and climate (Rivard et al., 2014). Though literature reviews of the various techniques that exist to quantify groundwater recharge are provided in Scanlon et al. (2002) and Healy (2011). Many methods can be used to estimate groundwater recharge, such as water budget methods, modeling methods, tracer methods, and methods based on surface water interaction studies. The latter is based on the estimation of groundwater discharge to surface water, mainly by streambed seepage determination, streamflow duration curves, or streamflow hydrograph separation (Scanlon et al., 2002). The recharge amount (in mm yr⁻¹) is then typically obtained by dividing measured or estimated discharge flow by the surface drainage area at the measurement site. This procedure assumes that aquifer boundaries coincide with watershed boundaries, and consequently that the area of the aquifer that contributes to groundwater discharge is equal to the surface drainage area (Kuniansky, 1989; Rutledge, 1998, 2007). However, this assumption must be considered carefully, as groundwater basins and watershed boundaries can differ drastically (Tiedeman et al., 1997). Miscalculation of the aquifer contributing area leads to a proportional error in recharge estimate.

Although the groundwater inflow to streams is often taken into account in water budgets, it is less commonly considered for surface water bodies, probably due to the greater difficulty of quantifying groundwater discharge in these settings. However, in recent years some studies have proven that groundwater flow into lakes can be reliably quantified. Interactions between lakes and groundwater depend on geology, soil, and sediment properties, and also on hydraulic gradient, which is strongly dependent on climatic conditions and recharge (Winter, 1999). Therefore, variation in groundwater fluxes may indicate a change in recharge in the lake catchment (Meinikmann et al., 2013).

In Quebec (Canada), more than 10% of the surface is covered by freshwater, with more than one million lakes known to exist. In many cases, these are connected to underlying aquifers. However, lake-groundwater interactions are highly dynamic throughout the year, and, even if it is now possible to quantify groundwater inflow with a reasonable degree of confidence, it is difficult to determine how and to what extent lakes can be sensitive to changes in groundwater recharge. The lake water isotopic composition has been

proven to be particularly useful for determining water balance parameter controls under changing conditions. For example, as shown in Turner et al. (2010), lake isotopic composition can highlight that (i) reduced winter precipitation could cause snowmelt-dominated lakes to become rainfall-dominated lakes, or that (ii) during longer ice-free seasons, lakes that are mainly rainfall-dominated (but also potentially snowmelt-dominated) may turn into evaporation-dominated lakes. Moreover, among all the methods used to quantify groundwater inflow to lakes, isotopic balances appear to be especially well adapted for quantifying groundwater flux variations on seasonal and yearly timescales (Arnoux et al., 2017a). Water stable isotopes are therefore expected to be very useful for monitoring seasonal and interannual variations in the water budget under changing recharge conditions.

The impact of climate change on groundwater recharge is not easy to determine, because of the complexity of interactions and processes involved, and can vary vastly depending on regions (Rivard et al., 2014; Crosbie et al., 2013). In addition, it is predicted to shift differentially under various climate scenarios and models (Jyrkama and Sykes, 2007; Levison et al., 2014). In Canada, highly variable recharge rates have been proposed in previous studies, for example, for the 2050 horizon (mainly the period 2041–2070) relative to modern (2000–2015) or past recharge rates (1950–2010), depending on study site, scenario, and model: +10 to +53% in the Grand River watershed, Ontario (Jyrkama and Sykes, 2007), -41 to +15% in the Chateauguay River watershed, Quebec (Croteau et al., 2010), -6 to +58% in the Otter Brook watershed, New Brunswick (Kurylyk and MacQuarrie, 2013), -4 to +15% at Covey Hill, Quebec (Levison et al., 2014), +14 to +45% in the Annapolis Valley, Nova Scotia (Rivard et al., 2014), and -28 to +18% for the Magdalen Islands, Quebec (Lemieux et al., 2015).

Recharge fluctuations can also impact lake water quality by changing groundwater fluxes, which are closely linked to phosphorous (P) loading to lakes. It is known that lake water quality is mainly driven by variations in P load, since it plays a critical role in limiting lake primary productivity and algal biomass, which in turn regulate lake trophic status. Increasing P concentration in the water column is the primary factor responsible for accelerated eutrophication and associated algae blooms (Schindler, 1977; Wang et al., 2008). At sites without urban drainage or point P sources, such as sewage treatment plants, domestic waste from septic systems may represent the largest anthropogenic source of P to lakes on the Canadian Shield (Dillon and Evans, 1993). Increases in shoreline development and population, combined with groundwater fluxes variations, can clearly impact lake quality, but still remain to be quantified.

For the present study, 10 lakes in southern Quebec were sampled to quantify their yearly groundwater inflows (see Arnoux et al., 2017a for more details), and one of these lakes was sampled over the course of a year to quantify its monthly groundwater inflows (see Arnoux et al., 2017b for more de-
tails). Small kettle lakes without surface inlets located in fluvioglacial deposits, and that are most likely well connected to shallow unconfined aquifers, were specifically targeted. The two main objectives of this study are (i) to determine how future groundwater recharge changes might affect lake water balance and geochemistry and (ii) to assess whether stable isotopes might be an effective tool for identifying lakes that are susceptible to change or are undergoing changes in water balance and water quality. To address these objectives, seasonal models of water and isotopic budgets were established for several lakes, and the models were then forced with future yearly and monthly timescale climate data from predictive global models to simulate anticipated conditions. Climate outputs of the Canadian Regional Climate Model were used, based on scenarios RCP 4.5 and RCP 8.5 (Moss et al., 2010; IPCC, 2014). It was assumed that recharge fluctuation was the main parameter influencing groundwater fluxes into lakes, and thus a percentage of recharge change leads to the same percentage of change of groundwater fluxes to lakes. Different recharge scenarios, which were translated into changes in groundwater inflow, were then tested to determine changes in water budget and isotopic evolution of the lakes. Predicted changes in recharge were then compared to predicted population growth in the study areas to discuss lake quality evolution. After the evolution of the lake geochemical signature is determined, the ways in which lakes that are connected to groundwater can be used to identify changes in groundwater recharge can also be determined, as can whether or not the isotopic composition of lakes can serve as an effective indicator of change or variability.

# 2 Method

## 2.1 Study sites

The 10 lakes chosen are located in four regions of southern Quebec characterized by contrasting climatic conditions: Laurentides (LAU), Outaouais (OUT), Abitibi-Témiscamingue (AT), and Saguenay–Lac-Saint-Jean (SAG). These kettle lakes, located in coarse-grained (sand and gravel) fluvioglacial deposits, are specifically targeted in this study, because they (i) are small enough to be sensitive to environmental changes on a short timescale, (ii) do not have permanent surface inflow streams and so are largely groundwater dependent, (iii) are generally characterized by predictable and uniform geomorphological features, and (iv) are likely connected to shallow, unconfined aquifers (Arnoux et al., 2017a; Isokangas et al., 2015). Kettle lakes originate as depressions in the landscape formed following the melting of ice blocks buried in the ground after glacial retreat of the Late Glacial to Holocene transition period (from -12 to -7 kyr). These kettle holes, becoming kettle lakes when they are filled with water, are mainly found in fluvioglacial deposits, such as outwash plains, deltas, eskers, and kame terraces (Benn and Evans, 2011). Figure 1 shows the locations of the 10 lakes analyzed here. Their main characteristics are described in Table 1.

## 2.2 Lake isotopic composition

### 2.2.1 Sampling

Water samples from each lake were retrieved during two field campaigns, in June–July and October–November 2014. When physicochemical parameters, measured in situ along the water column, revealed a well-mixed lake, the lake was considered to be homogeneous, and only one sample was collected, from close to the lake bottom, at its greatest depth. Otherwise, for stratified periods, two samples were collected: one from the top of the epilimnion and one from the base of the hypolimnion, in order to obtain the complete range of isotopic composition variation. Whenever possible, groundwater was sampled from private wells located in the vicinity of the studied lakes. Untreated groundwater samples were collected from residential wells from the tap after purging approximately 3 times the well volume.

Samples were transported in a cooler, and subsequently stored at 5‰ C until analyses were performed. Water stable isotopic compositions were measured with a Laser Water Isotope Analyser (OA ICOS DLT, Los Gatos Research, now ABB) at the GEOPS Laboratory (University of Paris-Sud–Paris-Saclay, France). The reproducibility of the sample measurement is ±1‰ versus VSMOW for δ²H and ±0.2‰ versus VSMOW for δ¹⁸O. Results are reported in δ values, representing deviations in per mill (‰) from the isotopic composition of the international standard (Vienna Standard Mean Ocean Water, VSMOW), such that δ²H or δ¹⁸O = ((Rsample/RVSMOW) - 1) × 1000, where R refers to ²H or ¹⁸O/¹⁶O ratios.

One of the lakes, Lake Lacasse, was sampled in more detail throughout 2015–2016. Water samples were collected from the lake at 2-week to 1-month intervals, mainly from the deepest part of the lake, and at depth intervals of 1 to 2 m in order to monitor the vertical heterogeneity of the water column. Groundwater was sampled twice from eight private wells in the vicinity of the lake (see Arnoux et al., 2017b for more details).

### 2.2.2 Water mass balance

The lake water budget is defined as follows:

$$ \frac{dV}{dt} = I - E - Q, \qquad (1) $$

where V is the volume of the lake (m³), t is time (days); E is evaporation (m³ day⁻¹); I is the instantaneous inflow (m³ day⁻¹), corresponding to the sum of upstream surface inflow (I_S; zero for the studied lakes because they do not have surface inlets), runoff (I_R; considered negligible because of
<table><caption>**Table 1.** Main lake characteristics.</caption><thead><tr><th>Region</th><th>ID</th><th>Lake name</th><th>Lake surface area (10<sup>3</sup> m<sup>2</sup>)</th><th>Lake volume (10<sup>3</sup> m<sup>3</sup>)</th><th>Catchment Area (10<sup>3</sup> m<sup>2</sup>)</th></tr></thead><tbody><tr><td>AT</td><td>1</td><td>Clair</td><td>115</td><td>695</td><td>2646</td></tr><tr><td>AT</td><td>2</td><td>Paix</td><td>41</td><td>97</td><td>796</td></tr><tr><td>AT</td><td>3</td><td>Sauvage</td><td>44</td><td>142</td><td>89</td></tr><tr><td>OUT</td><td>4</td><td>Lachigan</td><td>33</td><td>142</td><td>336</td></tr><tr><td>OUT</td><td>5</td><td>Lacroix</td><td>236</td><td>1080</td><td>772</td></tr><tr><td>OUT</td><td>6</td><td>Lanthier</td><td>25</td><td>125</td><td>1134</td></tr><tr><td>LAU</td><td>7</td><td>Lacasse</td><td>27</td><td>67</td><td>148</td></tr><tr><td>SAG</td><td>8</td><td>Beau Portage</td><td>42</td><td>271</td><td>364</td></tr><tr><td>SAG</td><td>9</td><td>Girard</td><td>67</td><td>679</td><td>211</td></tr><tr><td>SAG</td><td>10</td><td>Ludovic</td><td>94</td><td>400</td><td>1829</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Locations of the study lakes (circles) and sources of climate data (triangles).</figcaption></figure>

the permeable nature of the sandy soils of kettle lakes), groundwater inflow ($I_G$), and precipitation on the lake surface ($P$); and $Q$ is the outflow (m$^3$ day$^{-1}$), which is the sum of surface ($Q_S$) and groundwater ($Q_G$) outflows. Under constant atmospheric and hydrologic conditions, a steady state is assumed (Gibson et al., 2016), implying that $dV/dt = 0$. Therefore $I_G = Q_S + Q_G + E - P$ for the entire lake.

### 2.2.3 Stable isotopic mass balance

Considering water stable isotopes, the lake isotopic mass balance is as follows:

$$V \frac{d\delta_L}{dt} + \delta_L \frac{dV}{dt} = I \delta_I - E \delta_E - Q \delta_Q, \quad (2)$$

where $\delta$ is isotopic composition of the lake ($\delta_L$; equals to the mean if the lake is stratified – see Arnoux et al., 2017a and b, for more details about lake isotopic compositions used in the model); total inflow ($\delta_I$), which includes runoff ($\delta_R$), precipitation ($\delta_P$), surface inflow ($\delta_S$), and groundwater inflow ($\delta_G$); and total outflow ($\delta_Q$), which includes surface ($\delta_{QS}$)

and groundwater ($\delta_{QG}$) outflows. The isotopic composition of evaporating water ($\delta_E$) was estimated using the Craig and Gordon (1965) model, expressed by Gonfiantini (1986) as follows:

$$\delta_E = \frac{(\delta_L - \varepsilon^+) / \alpha^+ - h\delta_A - \varepsilon_K}{1 - h + 10^{-3}\varepsilon_K}, \quad (3)$$

where $h$ is the relative humidity at the lake surface, $\delta_A$ is the local isotopic composition of the atmospheric moisture (‰), $\alpha^+$ is the equilibrium isotopic fractionation, $\varepsilon^+ = (\alpha^+ - 1) \cdot 1000$ is the equilibrium isotopic separation (‰), and $\varepsilon_K = C_K (1 - h)$ is the kinetic isotopic separation (‰), with $C_K$ being the ratio of molecular diffusivities between heavy and light molecules (Gibson et al., 2016).

In this study, $C_K$ values were considered to be representative of fully turbulent wind conditions and a rough surface for both oxygen ($C_K = 14.2\%$) and hydrogen ($C_K = 12.5\%$), based on experimental data (Horita et al., 2008). For calculating equilibrium fractionation factors, experimental values of Horita and Wesolowski (1994) were used:

$$\alpha^+ (^{18}\text{O}) = \exp(-7.685/10^3 + 6.7123/T - 1666.4/T^2 + 350410/T^3), \quad (4)$$

$$\alpha^+ (^{2}\text{H}) = \exp(1158.8 \times T^3/10^{12} - 1620.1 \times T^2/10^9 + 794.84 \times T/10^6 - 161.04/10^3 + 2999200/T^3), \quad (5)$$

where $T$ is temperature (K). The isotopic composition of atmospheric moisture ($\delta_A$, ‰) was calculated assuming equilibrium isotopic exchange between precipitation and vapor:

$$\delta_A = \frac{\delta_P - \varepsilon^+}{1 + 10^{-3}\varepsilon^+}, \quad (6)$$

where $\delta_P$ (‰) is the mean annual isotopic composition of precipitation. Assuming well-mixed conditions in the lake, the combination of Eq. (3) and Eq. (2) yields the following:
$$V \frac{d\delta_L}{dt} + \delta_L \frac{dV}{dt} = P \delta_P + I_G \delta_G - Q \delta_L \\ - \frac{E}{1 - h + 10^{-3} \varepsilon_K} \left( \frac{\delta_L - \varepsilon^+}{\alpha^+} - h \delta_A - \varepsilon_K \right). \quad (7)$$

A steady state was assumed, such that $dV/dt = 0$. Equation (7) can therefore be simplified to the follows:

$$V \frac{d\delta_L}{dt} = P \delta_P + I_G \delta_G - (P + I_G - E) \delta_L \\ - \frac{E}{1 - h + 10^{-3} \varepsilon_K} \left( \frac{\delta_L - \varepsilon^+}{\alpha^+} - h \delta_A - \varepsilon_K \right). \quad (8)$$

Resolving this calculation therefore allows isotopic composition of the lake water at time $t + dt$ to be determined, expressed as a function of its value at the previous time step, $t$, and two established parameters, $A$ ($\%c m^3 yr^{-1}$) and $B$ ($m^3 yr^{-1}$):

$$\delta_L^{t+dt} = \frac{A}{B} + \left( \delta_L^t - \frac{A}{B} \right) \exp \left( -\frac{B}{V} dt \right) \quad (9)$$

with

$$A = P \delta_P + I_G \delta_G - \frac{E}{1 - h + 10^{-3} \varepsilon_K} \left( -h \delta_A - \varepsilon_K - \varepsilon^+ / \alpha^+ \right) \quad (10)$$

$$B = P + I_G - E \left( 1 - \frac{1}{\alpha^+ (1 - h + 10^{-3} \varepsilon_K)} \right). \quad (11)$$

The monthly-mean isotopic composition of precipitation ($\delta_P$) was assessed in the four regions from the Global Network of Isotopes in Precipitation (GNIP) and Program for Groundwater Knowledge Acquisition (PACES) datasets. Future $\delta_P$ trends are uncertain; however, they have been shown to be mainly dependent on temperature evolution and local factors (Stumpp et al., 2014), and a recent study in Siberia showed that a long-term increase in precipitation $\delta^{18}O$ is close to the detection limit of the tracers (<1‰ per 50 years) (Butzin et al., 2014). Monthly current means were therefore used in the current simulations. The mean value of groundwater isotopic composition ($\delta_{Gi}$) was determined from the mean groundwater isotopic composition measured in wells, located in the same region and presenting no enrichment due to evaporation. The mean isotopic values used for groundwater are presented in Table 2.

The uncertainties associated with the Craig and Gordon (1965) model in the estimated isotopic composition of evaporating moisture ($\delta_E$) can be substantial, especially if relative humidity is greater than 0.8 (Kumar and Nachiappan, 1999). Moreover, a sensitivity analysis of $^{18}O$ isotopic balance of a small lake in Austria (Yehdegho et al., 1997) indicates that for flow-through, groundwater-dominated systems with limited evaporation, the isotopic composition of the lake water and the inflow water are the parameters critical to the overall uncertainty. Horita et al. (2008) recommended using

<table><caption>**Table 2.** Mean isotopic composition of groundwater obtained for the four regions, in degrees versus VSMOW.</caption><thead><tr><th>Region</th><th>δ<sup>18</sup>O</th><th>δ<sup>2</sup>H</th></tr></thead><tbody><tr><td>AT</td><td>-14.00</td><td>-101.3</td></tr><tr><td>OUT</td><td>-11.56</td><td>-81.6</td></tr><tr><td>LAU</td><td>-11.71</td><td>-80.9</td></tr><tr><td>SAG</td><td>-14.06</td><td>-103.1</td></tr></tbody></table>

time-averaged values of the parameters in the calculation of $\delta_E$ for the given period of interest. Moreover atmospheric parameters should be preferably evaporation-flux-weighted whereas liquid fluxes to a lake should be amount-weighted (Gibson, 2002; Gibson et al., 2016). Therefore, on an annual time step, $\delta_P$ is monthly precipitation-flux-weighted, except when it is used to estimate $\delta_A$; in this case, $\delta_P$ is monthly evaporation-flux-weighted. On a monthly timescale, monthly values are used for each parameter of the model, and evaporation is considered to be null during the ice-covered period. Moreover, in winter, when monthly-mean temperature is below zero, precipitation is assumed to be zero in the model. Then, when monthly temperature becomes equal to or higher than zero, accumulated precipitation and amount-weighted $\delta_P$ are added to the calculation during the melt period. Moreover, sensitivity tests on this model have been performed by Arnoux et al. (2017a, b) and show that it is mostly sensitive to $E, h$, and $\delta_G$.

## 2.3 Evolution scenarios

### 2.3.1 Climate models

Climatic parameters used in this study (evaporation, humidity, temperature, and precipitation) come from climate models. RCMs allow the downscaling of large-scale information from GCMs to scales closer to the watershed scale, leading to a better representation of surface forcings. In the present study, the fifth version of the Canadian RCM (CRCM5) was chosen, which has a 0.44‰ horizontal grid resolution (approx. 50 km; Sushima et al., 2010; Martynov et al., 2013; Šeparović et al., 2013). The CRCM5 is a grid-point model, based on a two-time-level, semi-Lagrangian, (quasi) fully implicit time discretization scheme (Alexandru and Sushima, 2016). The model includes a terrain-following vertical coordinate based on hydrostatic pressure (Laprise, 1991; Alexandru and Sushima, 2015), and an horizontal discretization on a rotated latitude-longitude, Arakawa C grid (Arakawa and Lamb, 1977; Alexandru and Sushima, 2015). Following CRCM4, changes that have been introduced into CRCM5 include, for example, evolution in the planetary boundary layer parameterization to suppress both turbulent vertical fluxes under very stable conditions and the interactively coupled one-dimensional lake model (Flake; Mironov et al., 2010; Martynov et al., 2012; Šeparović et al., 2013).
CRCM5 uses the Canadian Land-Surface Scheme (CLASS, version 3.5; Verseghy, 1991; Alexandru and Sushama, 2015). This model is described in detail in Martynov et al. (2013) and Šeparović et al. (2013).

The CRCMs were driven by the second-generation Canadian Earth System Model (CanESM2, improved from CanESM1; Arora et al., 2011), developed by the Canadian Center for Climatic Modeling and Analysis (CCCma). As explained in Šeparović et al. (2013), it consists of a fourth-generation atmospheric general circulation model CanAM4, coupled with (i) the physical ocean component OGCM4 developed from the NCAR CSM Ocean Model (NCOM; Gent et al., 1998), (ii) the Canadian Model of Ocean Carbon (CMOC; Christian et al., 2010), and (iii) Canadian Terrestrial Ecosystem Model (CTEM; Arora and Boer, 2010). The CanAM4 is a spectral model employing T63 triangular truncation with physical tendencies calculated on a 2.81 linear grid and 35 vertical levels (Arora et al., 2011; Šeparović et al., 2013).

### 2.3.2 Climate data

Four greenhouse gas concentration scenarios (representative concentration pathways) have been adopted by the IPCC in its Fifth Assessment Report (AR5) in 2014: RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5. The scenarios selected for the present study are RCP 4.5 and RCP 8.5, for which predicted climate data are available until 2100 for the study regions. The RCP 4.5 scenario considers that long-term global emissions of greenhouse gases, land use, and land cover stabilize radiative forcing at $4.5 \text{ W m}^{-2}$ (approximately 650 ppm CO$_2$-equivalent) by the year 2100, without ever exceeding that value. The RCP 8.5 scenario corresponds to the highest greenhouse gas emissions pathway scenario, with gas emissions and CO$_2$ concentrations increasing considerably over time, thus leading to a radiative forcing of $8.5 \text{ W m}^{-2}$ by the end of the century (approximately 1370 ppm CO$_2$ equivalent). The defining characteristics of these scenarios are enumerated in Moss et al. (2010).

In order to connect these RCP forecasts to our study and to visualize trends, yearly mean data are presented in Fig. 2. Based on previous literature on recharge changes (see Sect. 2.3.3), a reference period (2010–2040) is compared to a future period (2041–2071). It is noted that both evaporation and temperature display increases between the reference and future periods for both scenarios, although it is more pronounced for RCP 8.5. Moreover, precipitation and relative humidity do not show clear trends. However, it seems that precipitation variability will increase overall for both scenarios, although this is more pronounced for RCP 8.5. Moreover, the southern regions (i.e., OUT and LAU) have higher temperatures than the northern regions (i.e., AT and SAG), and precipitation is higher in LAU than in the other three regions. On a monthly timescale, surface temperatures in LAU show an increasing monthly trend, whereas evaporation in-

creases mainly during summer and stays relatively constant the rest of the year (data not shown). Meanwhile, precipitation does not show any clear trend. However, as temperatures increases in winter, melt periods likely will shift more frequently occur earlier in the year.

### 2.3.3 Recharge evolution

The mean annual recharge for each lake basin was obtained by dividing the lake drainage area by the calculated mean annual groundwater inflow to the lake (Meinikmann et al., 2013). In this study, recharge evolution is thus expressed in terms of changes in groundwater inflow to the lakes.

In the first step, recharge is assumed to be constant for the 2006–2014 period. Over this period, recharge is adjusted to fit the calculated lake isotopic compositions to those measured. In the second step, the results of Rivard et al. (2014) was chosen for the simulation of recharge scenarios, since this study focusses on the Annapolis Valley (Nova Scotia, Canada), not far from southern Quebec and with a similar latitude, geology, and climate. Therefore, the future recharge dynamics determined for the Annapolis Valley are assumed to be similar to those of the present study sites. Rivard et al. (2014) found that all scenarios predict an annual recharge to the aquifer within the range of +14 to +45 % higher than at present by 2041–2071. They also predict, on a seasonal basis, that recharge will undergo (i) a marked decrease in summer (from 4 to 33 %) and (ii) a spectacular increase in winter (more than 200 %), due to an earlier melt period starting date.

The following section focused firstly on monthly lake isotopic composition evolution (Sect. 3.1) and secondly on yearly lake isotopic composition evolution (Sect. 3.2). Monthly and yearly values are compared for the two standard periods (i.e., for reference (2010–2040) and future (2041–2071) periods).

For the first part of the study, Lake Lacasse, located in the LAU region, has been chosen, since it was subject to continuous monitoring (Arnoux et al., 2017b). Its groundwater inflow and variability has therefore already been well constrained throughout the year 2015–2016 (Fig. 3b). For this lake, the model was run from 2006 to 2071, and four different recharge evolution scenarios were applied to the 2041–2071 period, following the predictions of Rivard et al. (2014) for scenarios S1 and S2, as described below.

- NC: no change in recharge (groundwater inflow follows the pattern described in Fig. 3, obtained from Arnoux et al., 2017b);

- S0: a recharge decrease of 33 % during the summer period (from June to October);

- S1: a 200 % increase in recharge during the melt period (from January to March) and a 4 % decrease in the summer period;
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Climate data for the reference (Ref; 2010–2040) and future (Fut; 2041–2071) periods, obtained from CRCM5-CanESM2, with RCP 4.5 (left) and RCP 8.5 (right) scenarios for the four different study areas. The variables are (a) surface air temperature, (b) surface water evaporation (obtained from surface heat flux), (c) surface relative humidity (obtained from surface specific humidity), and (d) precipitation (Martynov et al., 2013; Šeparović et al., 2013); box whiskers describe median, first, and third quartiles as well as maximum and minimum values.</figcaption></figure>

- S2: a 200% increase in recharge during the melt period and a 33% decrease during the summer period.

For the second part, three annual recharge evolution scenarios were tested, following the predictions of Rivard et al. (2014): no change (NC), a 14% increase (Low), and a 45% increase (High) in mean annual recharge.

## 2.4 Population growth

Variations in the quantity and/or quality of groundwater feeding lakes can obviously impact the geochemistry, and thus the water quality of lakes, especially for lakes displaying a high G index (the percentage of groundwater comprising the total lake inflow; Arnoux et al., 2017a). Moreover, in rural areas of Quebec, lake and groundwater quality is likely to be influenced by changes in population density. The population of Quebec is aging, and many seasonal residences (e.g., cottages) around lakes in rural areas are ex-

pected to become year-round residences. Furthermore, these residences are not connected to wastewater treatment plants; rather, owners have their own private wells for drinking water and private septic tanks with subsurface seepage beds for wastewater. The predicted population changes are summarized in Table 3. Population is mainly expected to increase in the southern regions (OUT and LAU), with a mean increase of 24 and 28%, respectively (ISQ, 2014; Table 3). Scenarios of population growth are compared with scenarios of recharge evolution for each lake to assess their future quality evolution.

# 3 Results and discussion

## 3.1 Monthly evolution of lake isotopic composition

Figure 3 shows the measured (see Arnoux et al., 2017b, and the Supplement for more details about measured values) and
<table><caption>**Table 3.** Predicted population growth in the different study regions in 2036 relative to 2011 numbers, according to three different scenarios (ISQ, 2014).</caption><thead><tr><th rowspan="2">Region</th><th colspan="3">Scenarios</th></tr><tr><th>Reference (%)</th><th>Low (%)</th><th>High (%)</th></tr></thead><tbody><tr><td>OUT</td><td>24</td><td>13</td><td>36</td></tr><tr><td>AT</td><td>5</td><td>0</td><td>10</td></tr><tr><td>LAU</td><td>28</td><td>21</td><td>34</td></tr><tr><td>SAG</td><td>0</td><td>-4</td><td>4</td></tr></tbody></table>

modeled isotopic compositions of Lake Lacasse. It can be observed that the modeled values are more variable than the measured ones, undoubtedly due to the higher evaporation rate in the climatic model (459 mm) than that measured during the field monitoring period (204 mm). It is also shown that the model attributes greater weight to the contribution of the depleted snow value than in reality. This is probably due to the snow column (which is close to 0‰ C during the snow melt) being less dense than the lake surface water (which has a mean temperature of close to 4‰ C) and therefore bypassing the lake, flowing rapidly out of the lake outlet. In such a case, the snow does not influence the lake isotopic composition as much as the model predicts. Since similar results are obtained for δ²H values, only the δ¹⁸O results from the model will be presented in the following sections.

Lake Lacasse has a mean G index (i.e. the percentage of groundwater that makes up total lake inflows) of 69% during the reference period. Results for monthly simulations, with RCP 4.5 climate data, are illustrated in Fig. 4. Lake isotopic compositions are not significantly different between the reference and future periods if no change is applied to the recharge pattern (Fig. 4). Under scenarios S1 and S2, it can be observed that future δ¹⁸O is nearly 100% different from reference conditions during the first 2 months of the year (Fig. 4). It is at least 75% different for the month of March, but this month shows important variation during the future period. Throughout the rest of the year, ranges of variation are not completely different, but increasing or decreasing trends can be observed, depending on the season.

Indeed, Fig. 5 shows the monthly differences between mean lake δ¹⁸O in the reference period (which is the same for all scenarios) and mean lake δ¹⁸O in the future period, for the four recharge evolution scenarios.

During the year:

- regarding the reference period, the highest variation is observed in March for S1 (-1‰), S2 (-1‰), and NC (-0.4‰), after the melt period. For S0, the greatest change regarding the reference period is observed in September and October (+0.4‰), after the evaporation period;

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a) Isotopic composition of Lake Lacasse between June 2015 and May 2016, measured and modeled from stable isotopic mass balance model using climate data from climate model CRCM5 – CanESM2 and scenario RCP 4.5; (b) the pattern of groundwater inflow ($I_G$) to Lake Lacasse (Arnoux et al., 2017b).</figcaption></figure>

- regarding the NC future period, the greatest difference between winter recharge is in February (-0.6 and -0.5‰ for S1 and S2, respectively). This suggests that future changes in lake isotopic composition associated with recharge may be highest in February.

During the summer:

- regarding the reference period, the highest variation is in August for NC (+0.2‰), while it is in September and October for S0 (+0.4‰ for both months) and S2 (+0.2 and +0.3‰ in September and October, respectively). S1 does not show any variation;

- regarding the NC future period, the greatest change occurs in October for S0 (+0.3‰) and for S2 (+0.2‰), and in September for S1 (-0.1‰).

Results of scenario S2, characterized by the greatest changes in recharge in both summer and winter, highlight that the impact of decreasing recharge during summer attenuates the substantial impact of increasing recharge during winter. Indeed, during winter, S1 shows more depleted values than S2 (-0.5 versus -0.4‰ in January, and -0.8 and -0.7‰ with respect to the reference period for S1 and S2, respectively). Therefore, the more recharge decreases in the summer, the more lake isotopic composition increases in the summer, due to increased evaporation. Meanwhile, the more recharge increases in the winter, the more lake isotopic composition is depleted in the winter. If both phenomena occur
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Monthly Lake Lacasse isotopic composition, calculated using RCP 4.5 climatic data, for different periods and various recharge patterns: no change (NC), -33% in the summer (from June to October; S0), +200% during the melt period (from January to March) and -4% in the summer (S1), and +200% during the melt period and -33% in the summer (S2); (b) close-up of the winter months; (c) close-up of the summer months; box whiskers describe median, first, and third quartiles as well as maximum and minimum values.</figcaption></figure>

in a given year, the mean annual lake isotopic composition evolution will therefore not be expected to shift much, since their opposing impacts on lake isotopic composition will cancel each other out. As such, S1 is the scenario showing the highest variation in annual mean, of -3‰, compared with -2‰ for S2 and +2‰ for S0.

Based on these observations, it appears that isotopic signatures measured at the end of February and in September or October will provide information on the greatest changes during the winter and summer periods, respectively. The greatest changes in lake isotopic composition are likely to be at the end of the melt period.

Moreover, simulation results show that RCP 4.5 and 8.5 models provide similar results for Lake Lacasse isotopic composition evolution. Figure 6 shows the comparison of lake $\delta^{18}$O composition for both RCP climate scenarios, from 2010 to 2071, assuming the NC recharge scenario. In Fig. 6,

it can be observed that there is a small trend toward $\delta^{18}$O enrichment due to a higher evaporation rate, which is more pronounced for the RCP 8.5 than for the RCP 4.5 scenario. However, on a yearly timescale, the impact of evaporation increase in the summer seems to be attenuated by a precipitation increase throughout the rest of the year, likely implying that these climate changes result in a nearly non-measurable impact on lake isotopic composition evolution.

Finally, all these results show that extreme caution is required when interpreting trends in lake isotopic composition, and that their interpretation requires (i) a minimum background knowledge – at least 1 year of data – of lake isotopic composition evolution in relation to its hydrological balance and (ii) an accurate evaluation of weather data variability in the year of monitoring, with respect to their annual means for the study lake. A long-term change in recharge will definitely impact lake isotopic composition, but the lake is also
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Differences between mean Lake Lacasse $\delta^{18}$O in the reference period and mean Lake Lacasse $\delta^{18}$O in the future period, for the RCP 4.5 climate scenario and four scenarios of recharge evolution: no change (NC), -33% in the summer (from June to October; S0), +200% during the melt period (from January to March) and -4% in the summer (S1), and +200% during the melt period and -33% in the summer (S2).</figcaption></figure>

sensitive to changes in other water budget parameters. It may therefore still be difficult to definitively isolate the effect of recharge over long time periods. As such, it is also important to consider evolution in the yearly mean lake isotopic composition.

## 3.2 Annual evolution of lake geochemistry

### 3.2.1 Isotopic signature evolution

The model was run for the 10 study lakes, including Lake Lacasse (Table 1 for main lake characteristics). Figure 7 illustrates differences in $\delta^{18}$O in the reference period compared to the future period for lakes which have a range of G indices (see Arnoux et al., 2017a for more details about lakes measured values). It can be observed that, if the recharge is set as constant from 2010 to 2071 (NC recharge scenario), there is no significant difference between the reference and future period (Fig. 7), although evaporation shows a significant increase with time. The lack of a trend is probably mitigated by concurrent shifts in precipitation (Fig. 2). Without considering changes in groundwater inflow, it appears that lake isotopic composition will be impacted at least as much by changes in precipitation as by changes in evaporation.

Figure 7 illustrates that the range of lake isotopic composition variation depends significantly on climate conditions, lake volumes, and their associated G indices. It can be observed that lakes with a low G index and a small volume have higher potential variability in isotopic composition regarding climatic variations than those with a high G index and high volume. For example, for two lakes with a similar mean G index, such as Lake Ludovic (SAG; G index = 51%) and Lake Lacroix (OUT; G index = 53%), the former is

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between monthly results in $\delta^{18}$O for both scenarios RCP 4.5 and 8.5 for the 2010–2071 period.</figcaption></figure>

expected to have a greater spread in isotopic compositions than the latter, even though the SAG region will likely undergo less evaporation increase compared with the OUT region (Fig. 2). This difference is due to the lower volume of Lake Ludovic ($V = 400\,000\,\text{m}^3$), compared with Lake Lacroix ($V = 1\,080\,000\,\text{m}^3$; Table 1). In addition, when lakes have a high G index, the groundwater flux tends to buffer lake isotopic variations, and so they tend to be less sensitive to changes in climate data. The dominant control on lake isotopic variability therefore appears to be the G index. Another example is Lake Lanthier, which has a smaller volume ($V = 125\,000\,\text{m}^3$) and a higher G index (G index = 94%), and therefore shows a limited range of isotopic variation compared with Lake Lacroix, although both are located in the OUT region (Fig. 7).

If a changing recharge scenario is applied, a decreasing trend in lake isotopic composition is clearly observed (Fig. 7). However, it is also shown that lakes are sensitive to large changes in annual recharge (+45%), but the differences are not significant if a smaller change (+14%) occurs. Moreover, as the percentage of recharge change applied in the model is the same for all lakes, it can be observed that the trend intensity will depend on four main parameters: lake catchment size (which controls the intensity of the flux change), the region (which underlies climate condition), lake volume (which impacts the range of variation), and the G index. However, a relationship is only found with the latter.

Figure 8 illustrates variations in mean lake $\delta^{18}$O versus G index in both reference and future periods. As shown, lake isotopic composition is more sensitive to changes in recharge for G indices ranging from 50 to 80%, with a maximum of
<figure><img src="image_10.png" /></figure>

**Figure 7.** Reference period (Ref; 2010–2040) lake $\delta^{18}$O composition and that corresponding to three different future period (2041–2071) recharge scenarios: no change (NC), +14 % (Low), and +45 % (High), for RCP 4.5 (top) and RCP 8.5 (bottom) scenarios. The values in brackets correspond to the mean G index (percentage of groundwater flow in the total inflow) for each lake calculated for the reference period; left panels show OUT and LAU regions, middle panels AT, and right panels SAG; box whiskers describe median, first, and third quartiles as well as maximum and minimum values.

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

**Figure 8.** Differences between mean lake $\delta^{18}$O in the reference period (2010–2040) and future period (2041–2071), for the higher recharge change scenario, versus lake G indices. RCP 4.5 (black dots) and 8.5 (grey dots) scenarios are represented.

sensitivity observed for a G index of around 65 %. It can also be observed that RCP 8.5 predicts a more depleted isotopic composition than does RCP 4.5. This implies that for the same recharge scenario, variations in precipitation and melt period (duration and time in the year) may impact the lake

isotopic evolution more than precipitation. Finally, the polynomial relationship between the two variables in Fig. 8 highlights that the G index drives the response of lake isotopic composition to changes in recharge.

### 3.2.2 Lake quality evolution

As has been shown previously, isotopic composition can be sensitive to future recharge changes. This section is now a discussion about how lake water quality could be impacted by these future recharge changes depending on the location of the lake. In the study areas, one of the principal concerns about lake water quality, today and in the future, is to avoid blue-green algae blooms in limiting P loads to lakes. This study does not take into account several parameters that may impact blue-green algae blooms in lakes, such as the lake water biogeochemistry, chemical threshold processes, thermal and oxygen stratification, and warming of the water column. The purpose here is to show in which cases lakes could be under risk of too-high P loading, and therefore at risk of a decrease in water quality, depending on their catchment evolution as a function of recharge and population evolution. In
such cases it will be important to prevent P loads and therefore changes on the lake catchment.

Turning to the predictions of population growth summarized in Table 3, population is predicted to increase mainly in the southern regions, OUT and LAU, with a mean increase of 24 and 28% by 2036, respectively (ISQ, 2014). Assuming an identical per capita P load, total P load in groundwater originating from wastewater should increase by the same percentage.

Domestic sewage is the main contribution of anthropogenic sources to the total P load for most of Canadian lakes (Dillon and Evans, 1993; Paterson et al., 2006). The total P load from sewage systems is a function of (i) the population and (ii) the annual P consumption per capita (Paterson et al., 2006). Assuming an effluent concentration of 9 mg L⁻¹ (considering reductions in the phosphate content of detergents) and a daily water usage of 200 L capita⁻¹ day⁻¹ (as was done by Paterson et al., 2006), the P contribution is estimated to be 0.66 kg capita⁻¹ yr⁻¹. Investigated lakes in the OUT and LAU regions collect sewage from 4 (Lake Lachigan), 53 (Lake Lanthier), 117 (Lake Lacroix), and 17 houses (Lake Lacasse) within their catchments. If two inhabitants per house are assumed, P loading to groundwater will be increased from 1 to 39 kg yr⁻¹ in the studied lakes in these areas.

The impact of this P load increase on lakes can then roughly be estimated based on the ratio of change in annual P load versus change in annual recharge, as illustrated in Fig. 9. For an increase in recharge, if ΔP/ΔR < 1, the change in recharge over the catchment, and thus the evolution of the groundwater inflow to the lakes, will be greater than the P variation. In such a case, the lake water quality may not be impacted by this P variation. However, if ΔP/ΔR > 1, the lake water quality will be impacted, and precautions should be taken to minimize the risk of blue-green algae blooms and consequent eutrophication. For the study regions (Fig. 9), if recharge increases 14% by 2036, as estimated by Rivard et al. (2014), lakes in the LAU and OUT areas will experience a decrease in their water quality. However, if the recharge change is closer to +45% (Rivard et al., 2014), lake water quality should not be worse than today, providing all other things remain equal and assuming the population growth forecasts are accurate.

## 4 Conclusions

The main objectives of this study were to determine how future trends in groundwater recharge can affect lake geochemistry, and to assess whether stable isotopes might be an effective tool for identifying lakes that are susceptible to change, or are undergoing changes, in their water budget and quality.

Firstly, climate predictions from both RCP 4.5 and 8.5 scenarios and their impacts on future lake isotopic composition have been considered. By 2050, temperature and evaporation

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Population growth predictions versus changes in recharge. The shaded area represents the scenarios for which lakes may be under risk of too-high P loading, and therefore at risk of a decrease in water quality. Dots represent lakes in the four study areas for three recharge scenarios.</figcaption></figure>

are expected to increase, and precipitation to exhibit a slightly increasing trend, all trends being more intense under the RCP 8.5 scenario. On a monthly time step, it has been highlighted that future lake isotopic signatures will be more depleted with respect to the reference period, mainly in March and February, because of an earlier melt period. In the summer, lake isotopic composition will be more enriched, mainly in August, due to the higher evaporation rate expected. However, future variations with respect to the reference period are smaller in the summer than in the winter. Scenario RCP 8.5 induces more intense monthly variations, but no significant difference in future lake isotopic signatures is observed on a yearly time step between the two scenarios. This means that enrichment caused by increased evaporation compensates for depletion induced by precipitation variation. It is therefore unclear whether lakes will be impacted more by increased evaporation or precipitation changes. Caution is therefore recommended in the interpretation of isotopic trends in lakes where background knowledge – for at least 1 year – of their isotopic composition evolution with respect to weather data and their hydrologic balance is lacking.

It has then been demonstrated that future lake isotopic composition will also depend on recharge fluctuations, in addition to climate conditions. On a monthly basis, the highest impact of recharge evolution on future lake isotopic composition will be in February. Moreover, if recharge decreases during the summer, the main difference will be observed at the end of the summer, after the evaporation period and before
recharge stops decreasing, in September or October. Therefore, to clearly identify future changes in recharge through the lake isotopic signature evolution, sampling only at the end of February and in September or October will provide information on the greatest changes for the winter and summer periods, respectively.

On an annual time step, modeled evolutions of lake isotopic composition can clearly be sensitive to both +45 and +14 % changes in recharge – less so, nevertheless, to the latter. The intensity of the future trend of lake isotopic composition will depend on four main parameters: lake catchment (which controls the intensity of the flux change), the region (which drives climate conditions), lake volume (which impacts the range of variation), and the G index (which is the dominant control on water balance conditions). Based on these model simulations, stable isotopes appear to be especially useful for detecting changes in recharge to lakes with a G index of between 50 and 80 %.

It is important to keep in mind that if both a winter increase and summer decrease in recharge occur during the same year, the trend in mean annual lake isotopic composition will be nullified, because seasonal variation is impacted in opposing directions, canceling out the signal at the yearly time step. Consequently, if no clear annual trend is observed, it does not mean that recharge is not changing. Nevertheless, mean annual lake isotopic compositions will be observed to be impacted by recharge evolution only if it evolves in the same way throughout the year for the most part (i.e., consistently decreasing or increasing). In light of these results, it is a monthly time step is strongly suggested in such investigations, since seasonal recharge fluctuations can be canceled out in the yearly signal. Moreover, it is important to note that runoff has been considered negligible for the study lakes but can be important for other lakes and, in these cases, this model could underestimate the effect of spring melt on future lake isotopic composition.

It is also shown that changes in water quality in groundwater-connected lakes depend substantially on lake location and on the intensity of recharge change. For the studied lakes, in the case of a +14 % recharge increase by 2036, lakes in LAU and OUT regions may experience altered water quality (driven by phosphorous loading), but no change is expected in the case of a +45 % recharge intensification. If the percentage of recharge increase is at least equal to the percentage of population growth around the lake, lake quality should not become degraded, but if not, recharge evolution should be considered in lake management. Lake water quality in the SAG and AT areas may not decrease when considering population growth predictions. However, this study does not take into account several parameters that can impact blue-green algae blooms in lakes, such as the lake water residence time, chemical threshold processes, and the warming of the water column (Planas and Paquet, 2016).

Finally, even if small groundwater-fed lakes will be sensitive to climate, and especially to recharge and anthropogenic

changes, it is still difficult to predict how their geochemistry will be impacted, as it is very reactive to each slight variation in water balance parameters. However, more indicators are now available to predict lake geochemistry evolution, mainly depending on their location and their G index. To go further, a recharge model adapted to lake catchments and coupled with melt dynamics, closely dependent on climate forecasts, could provide more details on lake geochemical evolution, for more sustainable lake management.

*Data availability.* Data used are available in Arnoux et al. (2017a, b) and in the Supplement. Detailed data required to reproduce the modelling results are available upon request.

**The Supplement related to this article is available online at https://doi.org/10.5194/hess-21-5875-2017-supplement.**

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* We thank all participants of this project which is part of a PhD project funded by University Paris-Saclay/Paris-Sud and by grants to Florent Barbecot (CRSNG and FQRNT). We thank all participants of the field campaigns, Viorel Horoi and students for their help on the field, the GRIL team for their logistical support, and people living around the study lakes for their hospitality and support. We thank also the team of Laxmi Sushima for their help with climate data. Comments from the editor and two anonymous referees have helped to improve this paper and we gratefully acknowledge them.

Edited by: Bettina Schaefl

Reviewed by: two anonymous referees

**References**

Alexandru, A. and Sushima, L.: Impact of land-use and land-cover changes on CRCM5 climate projections over North America for the twenty-first century, Clim. Dynam., 47, 1197, https://doi.org/10.1007/s00382-015-2896-3, 2016.

Arakawa, A. and Lamb, W. R.: Computational design of the basic dynamical processes of the UCLA general circulation model General circulation models of the atmosphere, Academic Press Inc, New York (A78-10662 01-47), 173–265, 1977.

Arnoux, M., Barbecot, F., Gibert-Brunet, E., Gibson, J., Rosa, E., and Noret, A., and Monvoisin, G.: Geochemical and isotopic mass balances of kettle lakes in southern Quebec (Canada) as tools to document variations in groundwater quantity and quality, Environ. Earth Sci., 76, 106, https://doi.org/10.1007/s12665-017-6410-6, 2017a.

Arnoux, M., Gibert-Brunet, E., Barbecot, F., Gillon, S., Gibson, J., and Noret, A.: Interactions between groundwater and season-
ally ice-covered lakes: Using water stable isotopes and radon-222 multilayer mass balance models, Hydrol. Process., 31, 2566–2581, 2017b.

Arora, V. K. and Boer, G. J.: Uncertainties in the 20th century carbon budget associated with land use change, Glob. Change Biol., 16, 3327–3348, https://doi.org/10.1111/j.1365-2486.2010.02202.x, 2010.

Arora, V. K., Scinocca, J. F., Boer, G. J., Christian, J. R., Denman, K. L., Flato, G. M., Kharin, V. V., Lee, W. G., and Merryfield, W. J.: Carbon emission limits required to satisfy future representative concentration pathways of greenhouse gases, Geophys. Res. Lett., 38, https://doi.org/10.1029/2010GL046270, 2011.

Benn, D. I. and Evans, D. J. A.: Sediment-landform associations in: Glaciers and Glaciation, 2nd edn. Paperback, Hodder Arnold, London, 421–533, 2010.

Butzin, M., Werner, M., Masson-Delmotte, V., Risi, C., Frankenberg, C., Gribanov, K., Jouzel, J., and Zakharov, V. I.: Variations of oxygen-18 in West Siberian precipitation during the last 50 years, Atmos. Chem. Phys., 14, 5853–5869, https://doi.org/10.5194/acp-14-5853-2014, 2014.

Christian, J. R., Arora, V. K., Boer, G. J., Curry, C. L., Zahariev, K., Denman, K. L., Flato, G. M., Lee, W. G., Merryfield, W. J., Roulet, N. T., and Scinocca, J. F.: The global carbon cycle in the Canadian Earth system model (CanESM1): Preindustrial control simulation, J. Geophys. Res.-Biogeosci., 115, https://doi.org/10.1029/2008JG000920, 2010.

Craig, H. and Gordon, L. I.: Deuterium and oxygen-18 in the ocean and marine atmosphere, in: Stable Isotopes in Oceanographic Studies and Paleotemperatures, edited by: Tongiorgi, E., Spoleto, Italy, 9–130, 1965.

Crosbie, R. S., Scanlon, B. R., Mpelasoka, F. S., Reedy, R. C., Gates, J. B., and Zhang, L.: Potential climate change effects on groundwater recharge in the High Plains Aquifer, USA, Water Resour. Res., 49, 3936–3951, https://doi.org/10.1002/wcr.20292, 2013.

Croteau, A., Nastev, M., and Lefebvre, R.: Groundwater Recharge Assessment in the Chateauguay River Watershed, Can. Water Resour. J., 35, 451–468, https://doi.org/10.4296/cwrj3504451, 2010.

Dillon, P. J. and Evans, H. E.: A comparison of phosphorus retention in lakes determined from mass balance and sediment core calculations, Water Res., 27, 659–668, https://doi.org/10.1016/0043-1354(93)09176-I, 1993.

Dragoni, W. and Sukhija, B. S.: Climate change and groundwater: A short review, Geological Society Special Publication, 1–12, 2008.

Gent, P. R., Bryan, F. O., Danabasoglu, G., Doney, S. C., Holland, W. R., Large, W. G., and McWilliams, J. C.: The NCAR climate system model global ocean component, J. Climate, 11, 1287–1306, 1998.

Gibson, J. J.: Short-term evaporation and water budget comparisons in shallow Arctic lakes using non-steady isotope mass balance, J. Hydrol., 264, 242–261, https://doi.org/10.1016/S0022-1694(02)00091-4, 2002.

Gibson, J. J., Birks, S. J., and Yi, Y.: Stable isotope mass balance of lakes: a contemporary perspective, Quaternary Sci. Rev., 131, 316–328, https://doi.org/10.1016/j.quascirev.2015.04.013, 2016.

Gonfiantini, R.: Environmental isotopes in lake studies. In: Fritz P, Fontes JCh (eds) Handbook of environmental isotope geochem-

istry, The terrestrial environment, B, vol 2, Elsevier, Amsterdam, 113–168, 1986.

Healy, R. W.: Estimating groundwater recharge, Estimating Groundwater Recharge, 1–245, 2011.

Horita, J. and Wesolowski, D.: Liquid-vapour fractionation of oxygen and hydrogen isotopes of water from the freezing to the critical temperature, Geochim. Cosmochim. Acta, 58, 3425–3437, 1994.

Horita, J., Rozanski, K., and Cohen, S.: Isotope effects in the evaporation of water: a status report of the Craig-Gordon model, Isotop. Environ. Health Stud., 44, 23–49, https://doi.org/10.1080/10256010801887174, 2008.

IPCC: the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, available at: https://www.ipcc.ch/report/ar5/, 2014.

ISQ: Institut de la statistique du Québec, available at: http://www.stat.gouv.qc.ca (last access: 2016), 2014.

Isokangas, E., Rozanski, K., Rossi, P. M., Ronkanen, A.-K., and Kløve, B.: Quantifying groundwater dependence of a sub-polar lake cluster in Finland using an isotope mass balance approach, Hydrol. Earth Syst. Sci., 19, 1247–1262, https://doi.org/10.5194/hess-19-1247-2015, 2015.

Jyrkama, M. I. and Sykes, J. F.: The impact of climate change on spatially varying groundwater recharge in the grand river watershed (Ontario), J. Hydrol., 338, 237–250, https://doi.org/10.1016/j.jhydrol.2007.02.036, 2007.

Kumar, B. and Nachiappan, R. P.: On the sensitivity of Craig and Gordon model for the estimation of the isotopic composition of lake evaporates, Water Resour. Res., 35, 1689–1691, https://doi.org/10.1029/1999WR900011, 1999.

Kuniansky, E. L.: Geohydrology and simulation of groundwater flow in the “400-foot”, “600-foot”, and adjacent aquifers, Baton Rouge area, Louisiana, Louisiana Department of Transportation and Development Technical Report, 49, 1989.

Kurylyk, B. L. and MacQuarrie, K. T. B.: The uncertainty associated with estimating future groundwater recharge: A summary of recent research and an example from a small unconfined aquifer in a northern humid-continental climate, J. Hydrol., 492, 244–253, https://doi.org/10.1016/j.jhydrol.2013.03.043, 2013.

Laprise, R.: The Euler equations of motion with hydrostatic pressure as an independent variable, Mon. Weather Rev., 120, 197–207, 1991.

Lemieux, J.-M., Hassaoui, J., Molson, J., Therrien, R., Therrien, P., Chouteau, M., and Ouellet, M.: Simulating the impact of climate change on the groundwater resources of the Magdalen Islands, Québec, Canada, J. Hydrol.: Reg. Stud., 3, 400–423, https://doi.org/10.1016/j.ejrh.2015.02.011, 2015.

Levison, J., Larocque, M., Fournier, V., Gagné, S., Pellerin, S., and Ouellet, M. A.: Dynamics of a headwater system and peatland under current conditions and with climate change, Hydrol. Process., 28, 4808–4822, https://doi.org/10.1002/hyp.9978, 2014.

Martynov, A., Sushama, L., Laprise, R., Winger, K., and Dugas, B.: Interactive lakes in the Canadian Regional Climate Model, version 5: The role of lakes in the regional climate of North America, Tellus A, 64, 16226, https://doi.org/10.3402/tellusa.v64i0.16226, 2012.

Martynov, A., Laprise, R., Sushama, L., Winger, K., Šeparović, L., and Dugas, B.: Reanalysis-driven climate simulation over CORDEX North America domain using the Canadian Re-
gional Climate Model, version 5: Model performance evaluation, Clim. Dynam., 41, 2973–3005, https://doi.org/10.1007/s00382-013-1778-9, 2013.

Meinikmann, K., Lewandowski, J., and Nützmann, G.: Laccustrine groundwater discharge: Combined determination of volumes and spatial patterns, J. Hydrol., 502, 202–211, https://doi.org/10.1016/j.jhydrol.2013.08.021, 2013.

Mironov, D., Heise, E., Kourzeneva, E., Ritter, B., Schneider, N., and Terzhevik, A.: Implementation of the lake parameterisation scheme FLake into the numerical weather prediction model COSMO, Boreal Environ. Res., 15, 218–230, 2010.

Moss, R. H., Edmonds, J. A., Hibbard, K. A., Manning, M. R., Rose, S. K., Van Vuuren, D. P., Carter, T. R., Emori, S., Kainuma, M., Kram, T., Meehl, G. A., Mitchell, J. F. B., Nakicenovic, N., Rishi, K., Smith, S. J., Stouffer, R. J., Thomson, A. M., Weyant, J. P., and Wilbanks, T. J.: The next generation of scenarios for climate change research and assessment, Nature, 463, 747–756, https://doi.org/10.1038/nature08823, 2010.

Paterson, A. M., Dillon, P. J., Hutchinson, N. J., Futter, M. N., Clark, B. J., Mills, R. B., Reid, R. A., and Scheider, W. A.: A review of the components, coefficients and technical assumptions of Ontario's Lakeshore Capacity Model, Lake and Reservoir Management, 22, 7–18, 2006.

Planas, D. and Paquet, S.: Importance of climate change-physical forcing on the increase of cyanobacterial blooms in a small, stratified lake, J. Limnol., 75, 201–214, https://doi.org/10.4081/jlimnol.2016.1371, 2016.

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

Rivard, C., Paniconi, C., Vigneault, H., and Chaumont, D.: A watershed-scale study of climate change impacts on groundwater recharge (Annapolis Valley, Nova Scotia, Canada), Hydrol. Sci. J., 59, 1437–1456, https://doi.org/10.1080/02626667.2014.887203, 2014.

Rutledge, A. T.: Computer programs for describing the recession of ground-water discharge and for estimating mean ground-water recharge and discharge from streamflow data – update, US Geol Surv Water-Resour Invest Rep, 98–4148:43, 1998.

Rutledge, A. T.: Update on the use of the RORA program for recharge estimation, Ground Water, 45, 374–382, 2007.

Scanlon, B. R., Healy, R. W., and Cook, P. G.: Choosing appropriate techniques for quantifying groundwater recharge, Hydrogeol. J., 10, 18–39, https://doi.org/10.1007/s10040-001-0176-2, 2002.

Schindler, D. W.: Evolution of phosphorus limitation in lakes, Science, 195, 260–262, 1977.

Šeparović, L., Alexandru, A., Laprise, R., Martynov, A., Sushama, L., Winger, K., Tete, K., and Valin, M.: Present climate and climate change over North America as simulated by the fifth-generation Canadian regional climate model, Clim. Dynam., 41, 3167–3201, https://doi.org/10.1007/s00382-013-1737-5, 2013.

Stumpp, C., Klaus, J., and Stichler, W.: Analysis of long-term stable isotopic composition in German precipitation, J. Hydrol., 517, 351–361, https://doi.org/10.1016/j.jhydrol.2014.05.034, 2014.

Sushama, L., Khaliq, N., and Laprise, R.: Dry spell characteristics over Canada in a changing climate as simulated by the Canadian RCM, Global Planet. Change, 74, 1–14, https://doi.org/10.1016/j.gloplacha.2010.07.004, 2010.

Tiedeman, C. R., Goode, D. J., and Hsieh, P. A.: Numerical simulation of ground-water flow through glacial deposits and crystalline bedrock in the Mirror Lake area, Grafton County, New Hampshire, US Geological Survey Professional Paper, 1572, 1997.

Turner, K. W., Wolfe, B. B., and Edwards, T. W. D.: Characterizing the role of hydrological processes on lake water balances in the Old Crow Flats, Yukon Territory, Canada, using water isotope tracers, J. Hydrol., 386, 103–117, https://doi.org/10.1016/j.jhydrol.2010.03.012, 2010.

versegy, D. L.: Class – A Canadian land surface scheme for GCMS. I. Soil model, Int. J. Climatol., 11, 111–133, https://doi.org/10.1002/joc.3370110202, 1991.

Wang, H. J., Liang, X. M., Jiang, P. H., Wang, J., Wu, S. K., and Wang, H. Z.: TP ratio and planktivorous fish do not affect nutrient-chlorophyll relationships in shallow lakes, Freshw. Biol., 53, 935–944, 2008.

Winter, T. C.: Relation of streams, lakes, and wetlands to ground-water flow systems, Hydrogeol. J., 7, 28–45, 1999.

Yehdegho, B., Rozanski, K., Zojer, H., and Stichler, W.: Interaction of dredging lakes with the adjacent groundwater field: An isotope study, J. Hydrol., 192, 247–270, https://doi.org/10.1016/S0022-1694(96)03102-2, 1997.