Quantifying floodwater impacts on a lake water budget via volume-dependent transient stable isotope mass balance

Isotope mass balance models have undergone significant developments in the last decade, demonstrating their utility for assessing the spatial and temporal variability in hydrological processes and revealing significant value for baseline assessment in remote and/or flood-affected settings where direct measurement of surface water fluxes to lakes (i.e. stream gauging) are difficult to perform. In this study, we demonstrate that isotopic mass balance modelling can be used to provide evidence of the relative importance of direct floodwater inputs and temporary subsurface storage of floodwater at ungauged lake systems. A volume-dependent transient isotopic mass balance model was developed for an artificial lake (named lake A) in southern Quebec (Canada). This lake typically receives substantial floodwater inputs during the spring freshet period as an ephemeral hydraulic connection with a 150 000 km² large watershed is established. First-order water flux estimates to lake A allow for impacts of floodwater inputs to be highlighted within the annual water budget. The isotopic mass balance model has revealed that groundwater and surface water inputs account for 60 %–71 % and 39 %–28 % of the total annual water inputs to lake A, respectively, which demonstrates an inherent dependence of the lake on groundwater. However, when considering the potential temporary subsurface storage of floodwater, the partitioning between groundwater and surface water inputs tends to equalize, and the lake A water budget is found to be more resilient to groundwater quantity and quality changes. Our findings suggest not only that floodwater fluxes to lake A have an impact on its dynamics during springtime but also significantly influence its long-term water balance and help to inform, understand, and predict future water quality variations. From a global perspective, this knowledge is useful for establishing regional-scale management strategies for maintaining water quality at flood-affected lakes, for predicting the response of artificial recharge systems in such settings, and for mitigating impacts due to land use and climate changes.

1 Introduction

Lakes are complex ecosystems which play a valuable economic, social, and environmental role within watersheds (Kløve et al., 2011). In fact, lacustrine ecosystems can provide a number of ecosystem services, such as biodiversity, water supply, recreation and tourism, fisheries, and sequestration of nutrients (Schallenberg et al., 2013). The actual benefits that can be provided by lakes depend on the water quality, and poor resilience to water quality changes can lead to benefit losses (Mueller et al., 2016). Globally, the quantity and quality of groundwater and surface water resources are known to be affected by land use (Baudron et al., 2013; Cunha et al., 2016; Lerner and Harris, 2009; Scanlon et al., 2005) and climate changes (Delpla et al., 2009). As both surface water and groundwater contribute to lake water balances (Rosenberry et al., 2015), changes that affect the surface water/groundwater apportionment can potentially modify or threaten lake water quality (Jeppesen et al., 2014). Understanding hydrological processes in lakes can help to depict the vulnerability and/or resilience of a lake to pollution (Rosen, 2015) and to invasive species (Walsh et al., 2016) and, thus, secure water quantity and quality over time for drinking water production purposes (Herczeg et al., 2003). In Quebec (Canada), there are an important number of municipal wells that receive contributions from surface water resources (i.e. lakes or rivers) and are, thus, performing unintentional (Patenaude et al., 2020) or intentional (Masse-Dufresne et al., 2019, 2021) bank filtration.

Over the past few decades, significant developments have been made in the application of isotope mass balance models for assessing the spatial and temporal variability in hydrological processes in lakes, most notably the quantification of groundwater and evaporative fluxes (Herczeg et al., 2003; Bocanegra et al., 2013; Gibson et al., 2016; Arnoux et al., 2017b). In remote environments, such as in northern Canada, the application of isotopic methods is particularly convenient as direct measurements of surface water and groundwater fluxes is time-consuming, expensive, and difficult (Welch et al., 2018). Isotopic mass balance models can notably be applied to ungauged lake systems to efficiently characterize the impacts of floods on water apportionment (Haig et al., 2020). While isotopic frameworks were successfully used to assess the relative importance of floodwater inputs to lakes (Turner et al., 2010; Brock et al., 2007), no attempt was made to evaluate the timing of the floodwater inputs and to differentiate between the role of (i) direct floodwater inputs and (ii) temporary subsurface storage of floodwater on a lake's annual water budget. In this study, direct inputs refer to the floodwater that enters a lake via the surface (e.g. by inundating and/or flowing through a stream), while temporary subsurface storage of floodwater encompasses the floodwater-like inputs that reach the lake via the subsurface (e.g. through floodplain recharge or bank storage).

To gain information on the timing of hydrological processes, one may use a transient and short time step isotopic mass balance. A previous study by Zimmermann (1979) used a transient isotope balance to estimate groundwater inflow and outflow, evaporation, and residence times for two young artificial groundwater lakes near Heidelberg, Germany, although these lakes had no surface water connections, and volumetric changes were considered negligible. Zimmermann (1979) showed that the lakes were actively exchanging with groundwater, which controlled the long-term rate of isotopic enrichment to isotopic steady state, but the lakes also responded to seasonal cycling in the magnitude of water balance processes. While these findings are informative, Zimmermann (1979) did not attempt to build a predictive isotope mass balance model but rather used a best-fit approach to obtain a solitary long-term estimate of water balance partitioning for each lake. Petermann et al. (2018) also constrained groundwater connectivity for an artificial lake near Leipzig, Germany, with no surface inlet or outlet. By comparing groundwater inflow rates obtained via stable isotope and radon mass balances on a monthly time step, Petermann et al. (2018) highlighted the need to consider seasonal variability when conducting lake water budget studies. Our approach builds on that of Zimmermann (1979) and Petermann et al. (2018), developing a predictive model of both atmospheric and water balance controls on isotopic enrichment and accounting for volumetric changes on a daily time step.

The main objective of this study is to provide evidence of the relative importance of direct floodwater inputs and temporary subsurface storage of floodwater at ungauged lake systems using an isotopic mass balance model. To do so, we first aim to establish an isotopic framework based on the local water cycle, to verify the applicability of isotopic mass balance in the present setting, as contrasting isotopic signatures are required between various water reservoirs and fluxes, including floodwater inputs. Second, we quantify the water budget according to two reference scenarios (A and B) to grasp the impact of site-specific uncertainties on the computed results. Then, we analyse the temporal variability in the groundwater inputs and the sensitivity of the lake to floodwater-driven pollution. Finally, we demonstrate the implications of floodwater-like subsurface inputs on the water balance partition.

The water balance is computed via a volume-dependent transient isotopic mass balance model, which is applied to predict the daily isotopic response of an artificial lake in Canada that is ephemerally connected to a 150 000 km² watershed during spring freshet. During these flood events, the surficial water fluxes entering the study lake are not constrained in a gaugeable river or canal but occur over a 1 km wide surficial flood area. Our study period spans a flood with an average recurrence interval of 100 years and is therefore an example of the response of the system to a major hydrological event.

2 Study site

2.1 Geological and hydrological settings

The study site is located in the area of Greater Montreal and borders the lake Deux Montagnes (further referred to as lake DM), which corresponds to a widening of the Ottawa River at the confluence with the St Lawrence River in Quebec (Canada; Fig. 1). The Ottawa River is the second-largest river in eastern Canada, draining a watershed of approximately 150 000 km² (MDDELCC, 2015). The water level of lake DM is partly controlled by flow regulation structures (e.g. hydroelectric dams) upstream on the Ottawa River. Lake DM water levels also show seasonal fluctuations in response to precipitation and snowpack melting over the Ottawa River watershed. High water levels at lake DM are typically observed during springtime (April–May) and, less prominently, during autumn (November–December), while lowest water levels normally occur at the end of the summer (September; Centre d'expertise hydrique du Québec, 2020).

Figure 1 (a–c) Location of the study site, relative to the Ottawa River watershed, lake Deux Montagnes (DM), and the area of Greater Montreal. (d) Location of lake A and lake B relative to lake DM and a schematic representation of the hydrogeological context. The grey dashed lines illustrate the approximative extent of the paleo valley. LA-S1 and LB-S1 are surface water sampling points at lake A and lake B, respectively. LA-P1 to LA-P4 correspond to vertical profile sampling locations at lake A. The maps were created from openly available data and used in accordance with the Open Government Licence – Canada, or the Open Data Policy, M-13-13 of the United States Census Bureau. Detailed source information is provided in Appendix A.

Lake A (2.79×10⁵ m²) and lake B (7.6×10⁴ m²) are two small artificial lakes created from sand-dredging activities and are located at approximately 1 km from the shore of lake DM. The dredging is still ongoing at lake A, while it ceased a few decades ago at lake B. Both lakes are approximately 20 m deep (Masse-Dufresne et al., 2019) and were excavated within alluvial sands which were deposited in a paleo valley extending in the NE–SW direction and carved into the Champlain Sea clays (Ageos, 2010). Lithostratigraphic data (i.e. well logs) suggest that the paleo valley is approximately 600 m wide and has a maximum depth of 25 m. Between lake DM and lake A, a thin layer (few centimetres to roughly 2 m) of alluvial sands are deposited on top the clayey sediments (Fig. S1 in the Supplement; Ageos, 2010).

Lake A is connected to a small stream (S1) with a mean and maximum annual discharge of 0.32 and 1.19 m³ s⁻¹, respectively (Ageos, 2010). Maximum discharge typically occurs during the month of April as S1 drains snowmelt water from a small watershed (14.4 km²; Centre d'expertise hydrique du Québec, 2019), whereas low flow is recorded for the rest of the hydrological year. For the springtime in 2017, the surface water flows from S1 are deemed negligible compared to the floodwater inputs and are thus not considered in this study.

In total, two channelized outlet streams (S2 and S3) allow water to exit lake A and flow towards lake DM. The direction of the surface water fluxes at S2 can be reversed if water level at lake DM exceeds both a topographic threshold at 22.12 m a.s.l. (above sea level; determined from a topographic land survey along S2) and the water level at lake A (Ageos, 2010). Flow reversal also occurs in S3, but the elevation of the topographic threshold is unknown.

Lake A and lake B both contribute to the supply of a bank filtration system which is composed of eight wells and is designed to supply drinking water for up to 18 000 people (Ageos, 2010). Typically, two to three wells are operated on a daily basis at a total pumping rate ranging from 4000 m³ d⁻¹ (in wintertime) to 7500 m³ d⁻¹ (in summertime; Masse-Dufresne et al., 2019). Although the operation of the bank filtration system does not form a complete hydraulic barrier between the two artificial lakes, it does lead to a lowering of the lake B water level to below that of lake A (Ageos, 2010).

2.2 Hydrodynamics of the major flood event

In 2017, a major flood event occurred in the peri-urban region of Montreal and was caused by the combination of intense precipitation and snowpack melting over the Ottawa River watershed (Teufel et al., 2019). Rapid water level rise at lake DM occurred in late February, early April, and early May at rates of approximately 0.11, 0.19, and 0.16 m d⁻¹, respectively. A historical maximum water level (i.e. 24.77 m a.s.l.) was reached on 8 May 2017, corresponding to a net water level rise of >2.7 m compared to early February (Fig. 2). High water levels at lake DM resulted in the inundation of the area between lake A and lake DM (Fig. 1d), and the surface water fluxes were not constrained in S2 and S3 but occurred over a 1 km wide area.

Figure 2 Daily mean water levels at lake A, lake DM, and observation well VP from 9 February 2017 to 25 January 2018. The grey shaded area corresponds to the floodwater control period. Q_max and Q_min indicate the timing of the adjusted maximum and minimum output from the lake.

The water level in lake A was equivalent to lake DM during the flood peak (on 8 May 2017) and daily mean water levels at lake A and lake DM show good correlation (R²=0.98 and p value < 0.01) for the observed period (27 April to 17 May 2017). Daily mean water levels at observation well VP and lake DM also follow a similar pattern from late February to late July 2017 (R²=0.93 and p value < 0.01). Considering the above and a visible hydraulic connection between the lake DM and lake A, the data indicate that the daily water level variations at observation well VP were controlled by lake DM from late February to late July 2017. The lake A water level was also presumably controlled by lake DM until late July 2017, but technical issues prevented confirmation (i.e. the logger in lake A broke on 17 May 2017).

Then, from August to late October 2017, the water level in lake DM was below the topographical threshold, and there is no similarity between the evolution of the water level at lake DM and observation well VP (R²=0.11 and p value > 0.01). It is, thus, possible to infer that the lake A water level was not controlled by lake DM from August to late October 2017. This is also supported by the manual measurement of the lake A water level in September.

The water level of lake DM exceeded the topographic threshold again from November 2017 to January 2018, but the daily mean water levels at lake DM and observation well VP show a moderate correlation (R²=0.63 and p value < 0.01). The manual measurements also indicate a discrepancy between lake DM and lake A water levels in December 2017 and January 2018. The weaker correlation between the water level measurements suggests that lake DM was not controlling the dynamics of the lake A water level. It is, thus, likely that lake A received little to no surface water inputs from lake DM from November 2017 to January 2018. In this context, surface water inflows from lake DM during autumn and winter are considered negligible in this study and not included in the developed stable isotope mass balance model (Sect. 4.2).

2.3 Conceptual model of lake A water balance

Based on the geological and hydrological setting of the study site (Sect. 2.1) and flood-specific considerations (Sect. 2.2), we established a conceptual model of the lake A water balance, as described below.

Considering that lake A is sitting in alluvial sands (i.e. a highly permeable material), it is assumed that groundwater inputs (I_G) and outputs (Q_G) contribute to the water budget. Although it is difficult to interpret the location of I_G, it appears evident that Q_G occur along the NE bank of lake A. In fact, there are subsurface fluxes across the sandy bank that contribute to the bank filtration system or discharge into lake B, as its water level has been lower since the initiation of the bank filtration system (Masse-Dufresne et al., 2019). Given the regional groundwater flow in the NE to SW, Q_G can also presumably occur along the SW bank of lake A. Besides, it is likely that little to no subsurface fluxes exist in the area between lake A and lake DM, where clayey sediments are found.

Figure 3 Schematic representation of the hydrological processes at lake A during (a) groundwater control and (b) floodwater control periods. Inputs include precipitation (P), surface water (I_S), and groundwater (I_G), while outputs include evaporation (E), surface water outflow (Q_S), and groundwater outflow (Q_G). The area between lake DM and lake A is flooded in (b), and I_S from lake DM contribute to the water balance of lake A.

For the study period, it is conceptualized that the direction of the surface water fluxes in S2 and S3 is from lake A to lake DM, except from 27 February to 8 May 2017. During this period (hereafter referred to as the floodwater control period), the water level of lake DM exceeds the topographic threshold, and lake A receives surface water inflow (I_S) from lake DM. Also, it is likely that the high water level in lake A imposed a hydraulic gradient at the lake–aquifer interface, which allowed for Q_G from the lake and inhibited I_G. Then, as lake A and lake DM water levels started to decrease (from 8 May 2017), it is assumed that water exits lake A as surface water outputs (Q_S) towards lake DM or as Q_G. Although lake DM's water level again exceeded the topographic threshold from November 2017 to January 2018, the weaker correlation between the water levels suggests that the lake A water level was not controlled by lake DM. In this context, we conceptualized that lake A water level variations are mainly controlled by groundwater flows (I_G and Q_G). Surface water inputs (I_S) are set to zero during this period (see Sect. 2.2).

To summarize, for the year 2017, the lake A water budget can be conceptualized with two distinct hydrological periods, i.e. (a) the groundwater control period and (b) the floodwater control period (Fig. 3). While the groundwater control period concerns most of the hydrological year, the floodwater control period only applies from 23 February to 8 May 2017. During the groundwater control period (Fig. 3a), it is assumed that groundwater inflows (I_G) and precipitation (P) constitute the total water inputs to lake A, while surface water inflows (I_S) are negligible. During this period, the outputs are occurring through evaporative fluxes (E), surface water outflows (Q_S), and groundwater outflows (Q_G). In contrast, it is assumed that I_S and P represent the total water inputs to lake A during the floodwater control period (Fig. 3b). High water levels at lake A impose a hydraulic gradient at the lake–aquifer interface, which allows for Q_G and inhibits I_G.

3 Methods

3.1 Field measurements

Pressure–temperature loggers (Diver^®; TD-Diver and CTD-Diver; Van Essen Instruments B.V., Delft, the Netherlands) were used to measure surface water levels at lake A and groundwater levels at observation well VP on a 15 min time step. Water levels were recorded from 27 April 2017 (after the ice cover melted) to 17 May 2017 at lake A and from 29 March 2017 to 25 January 2018 (except between 19 July and 6 August 2017) at observation well VP. All the level loggers' clocks were synchronized with the computer's clock when launching automatic measurements. This procedure was done via the Diver-Office 2018.2 software. Manual measurements of the water level were regularly performed to calibrate (relative to a reference datum) and validate the automatic water level measurements. A level logger was also used to measure on-site atmospheric pressure and perform barometric compensation on water level measurements. Also note that water levels in lake A were not continuously recorded after 17 May 2017 due to a logger failure, but manual water level measurements (in September 2017, December 2017, and January 2018) depict the general evolution of lake A's water level.

Mean daily water levels at lake DM were retrieved with permission from the Centre d'expertise hydrique du Quebec database (Centre d'expertise hydrique du Québec, 2020). Meteorological data were measured at land-based meteorological stations near the study site and obtained from the Environment and Climate Change Canada database (available online at https://www.weatherstats.ca/, last access: 30 September 2019). Daily air temperature, relative humidity, wind speed, dew point, and atmospheric pressure were measured at the Montréal–Mirabel International Airport station (45.68^∘ N, −74.04^∘ E; 18 km from the study site). Daily precipitation and solar radiation were measured at Sainte-Anne-de-Bellevue station (45.43^∘ N, −73.93^∘ E; 10 km from the study site) and Montréal–Trudeau International Airport station (45.47^∘ N, −73.75^∘ E; 17 km from the study site), respectively.

3.2 Water sampling and analytical techniques

Physico-chemical parameter measurements and water sampling were performed at lake A at approximately 0.3 m below the surface and 1 m from the lake shoreline (at LA-S1) on a weekly to monthly basis from 9 February 2017 to 25 January 2018. Physico-chemical parameters (including temperature, electrical conductivity, pH, and redox potential) were measured using a multiparameter probe (YSI Pro Plus 6051030 and Pro Series pH/ORP/ISE and Conductivity Field Cable 6051030-1; YSI Incorporated, Yellow Springs, OH, USA). Additionally, vertical profile measurements and depth-resolved water sampling were conducted on 9 February 2017, 17 August 2017, and 25 January 2018 (at LA-P1 to LA-P4). Lake A water sampling was performed in the northern part of the lake for logistical reasons and due to ease of accessibility. As horizontal homogeneity has been previously demonstrated by Pazouki et al. (2016), the water samples were deemed representative of the whole waterbody.

Floodwater was sampled at two locations (near S2 and S3) on 19 April 2017 and at lake DM on 10 May 2017. Water samples were also collected at the surface and at depth within lake B and at observation well Z16, which is upstream of lake B and, thus, representative of the regional groundwater contributing to the latter (Ageos, 2016).

Water samples were analysed for major ions, alkalinity, and stable isotopic compositions of water (δ¹⁸O and δ²H). Water was filtered in the field using 0.45 µm hydrophilic polyvinylidene fluoride (PVDF) membranes (Millex-HV; Millipore, Burlington, MA, USA) prior to sampling for major ions and alkalinity. From December to March, cold weather prevented field filtration, so this procedure was performed in the laboratory on the same day. All samples were collected in 50 mL polypropylene containers and cooled during transport to the laboratory. The samples were then kept refrigerated at 4 ^∘C until analysis, except for stable isotopes, which were stored at room temperature. Major ions were analysed within 48 h via ionic chromatography (ICS 5000 AS-DP Dionex; Thermo Fisher Scientific, Saint-Laurent, QC, Canada) at Polytechnique Montréal (Montreal, Quebec). The limit of detection was ≤0.2 mg L⁻¹ for all major ions. Bicarbonate concentrations were derived from alkalinity, which was measured manually in the laboratory according to the Gran method (Gran, 1952) at Polytechnique Montréal (Montreal, Quebec). On samples with measured alkalinity (n=12), the ionic balance errors were all below 8 %. The mean and median ionic balance errors were 1 %. Stable isotopes of oxygen and hydrogen were measured with a water isotope analyser with off-axis integrated cavity output spectroscopy (LGR-T-LWIA-45-EP; Los Gatos Research, San Jose, CA, USA) at Geotop-UQAM (Montreal, Quebec). In total, 1 mL of water was pipetted in a 2 mL vial and closed with a septum cap. Each sample was injected (1 µL) and measured 10 times. The first two injections of each sample were rejected to limit memory effects. A total of three internal reference waters (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=\mathrm{0.23}\pm \mathrm{0.06}$ ‰, $-\mathrm{13.74}\pm \mathrm{0.07}$ ‰ and $-\mathrm{20.35}\pm \mathrm{0.10}$ ‰; ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{1.28}\pm \mathrm{0.27}$ ‰, $-\mathrm{98.89}\pm \mathrm{1.12}$ ‰, and $-\mathrm{155.66}\pm \mathrm{0.69}$ ‰) were used to normalize the results on the Vienna Standard Mean Ocean Water–Standard Light Antarctic Precipitation (​VSMOW–SLAP) scale. A fourth reference water (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{4.31}\pm \mathrm{0.08}$ ‰; ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{25.19}\pm \mathrm{0.83}$ ‰) was analysed as an unknown to assess the exactness of the normalization. The overall analytical uncertainty (1σ) is better than ±0.1 ‰ for δ¹⁸O and ±1.0 ‰ for δ²H. This uncertainty is based on the long-term measurement of the fourth reference water and does not include the homogeneity nor the representativity of the sample.

3.3 Stable isotope mass balance

Stable isotope mass balances for lakes can either be performed based on (i) a well-mixed single layer model or (ii) a depth resolved multi-layered model. Arnoux et al. (2017c) performed a comparison of both methods and reported that well-mixed and depth-resolved multi-layered models yielded similar results and showed that groundwater inputs and outputs play an important role in lake water budgets. Arnoux et al. (2017c) further highlighted that the multi-layer model additionally allowed for the determination of groundwater flow with depth but required a temporally and depth-resolved sampling in order to ensure a thorough understanding of the stability/mixing of the different layers. Such time-consuming sampling and monitoring efforts are, however, often unrealistic in remote and/or flood-affected contexts. Additionally, Gibson et al. (2017) showed that the timing of the lake water sampling may introduce greater bias in a well-mixed isotopic mass balance model than the uncertainty related to the lake stratification. For these reasons, we opted to develop a well-mixed model in the context of this study. Note that, despite the biases underlying well-mixed models, this approach remains adequate for characterizing the relative importance of hydrological processes and is particularly useful for giving first-order estimates of water fluxes in ungauged basins.

The water and stable isotope mass balance of a well-mixed lake can be described, respectively as Eqs. (1) and (2) as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}=I-E-Q\text{(2)}& {\displaystyle }& {\displaystyle }V{\displaystyle \frac{\mathrm{d}{\mathit{\delta }}_{\mathrm{L}}}{\mathrm{d}t}}+{\mathit{\delta }}_{\mathrm{L}}{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}t}}=I{\mathit{\delta }}_{\mathrm{I}}-E{\mathit{\delta }}_{\mathrm{E}}-Q{\mathit{\delta }}_Q,\end{array}$$

where V is the lake volume, t is time, I is the instantaneous inflow, E is evaporation, and Q is the instantaneous outflow. I corresponds to the sum of surface water inflow (I_S), groundwater inflow (I_G), and precipitation (P). Similarly, Q is the sum of surface water outflow (Q_S) and groundwater outflow (Q_G). δ_L, δ_I, δ_E, and δ_Q are the isotopic compositions of the lake, I, E, and Q, respectively. In the context of this study, the balance equations can be simplified based on the conceptual model. During the groundwater control period, I_S=0 and, thus, $I=I_{\mathrm{G}}+P$ and ${\mathit{\delta }}_{\mathrm{I}}=({\mathit{\delta }}_{\mathrm{G}}{I}_{\mathrm{G}}+{\mathit{\delta }}_{\mathrm{P}}{I}_{\mathrm{P}})/I$. In contrast, I_G=0 during the floodwater control period, $I=I_{\mathrm{S}}+P$ and ${\mathit{\delta }}_{\mathrm{I}}=({\mathit{\delta }}_{I_{\mathrm{s}}}{I}_{\mathrm{S}}+{\mathit{\delta }}_{\mathrm{P}}{I}_{\mathrm{P}})/I$. Note that δ_G and ${\mathit{\delta }}_{I_{\mathrm{s}}}$ are the isotopic signatures of groundwater and surface water inputs, respectively.

The application of Eqs. (1) and (2) for both δ¹⁸O and δ²H is valid during the ice-free period and also assumes the constant density of water (Gibson, 2002). In this study, the potential impacts of the ice cover formation and melting are neglected, as the ice volume is likely to represent only a small fraction (<2 %) of the entire water body. Moreover, considering the ice water isotopic separation factor, i.e. 3.1 ‰ for δ¹⁸O and 19.3 ‰ for δ²H (O'Neil, 1968) and assuming well-mixed conditions, the lake water isotopic variation would be comprised within the analytical uncertainty. Also, floodwater inputs from lake DM were expected to be much more important and occurring simultaneously with ice melt during the freshet period.

Thus, a volume-dependent model is applied, as described in Gibson (2002). The change in the isotopic composition of the lake (δ_L) with f (i.e. the remaining fraction of lake water) can be expressed as Eq. (3) as follows:

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}_{\mathrm{L}}\left(f\right)={\mathit{\delta }}_{\mathrm{S}}-\left({\mathit{\delta }}_{\mathrm{S}}-{\mathit{\delta }}_{\mathrm{0}}\right)f^{\left[\frac{-(\mathrm{1}+mX)}{\mathrm{1}-X-Y}\right]},\end{array}$$

where $X=E/I$ is the fraction of lake water lost by evaporation, $Y=Q/I$ is the fraction of lake water lost to liquid outflows, m is the temporal enrichment slope (see Appendix B), δ₀ is the isotopic composition of the lake at the beginning of the time step, and δ_S is the steady-state isotopic composition the lake would attain if f tends to 0 (see Appendix B).

A step-wise approach is used to solve Eq. (3) on a daily time step. At each time step, recalculation of $f=V/V_{\mathrm{0}}$ is needed, where V is the residual volume at the end of the time step and V₀ the original volume at the beginning of the time step (or V^t−dt). Hence, Eq. (3) is based on the water level difference between 2 d. The water flux parameters (E, I, and Q) and isotopic signatures (δ_E, δ_A, δ_I, and δ_Q) are thus evaluated on a daily time step.

The flushing time (t_f) is defined as the ratio of the volume of water in a system to the rate of renewal (Monsen et al., 2002). In this study, t_f by groundwater inputs is considered and is expressed as follows:

$$\begin{array}{}\text{(4)}& t_{\mathrm{f}}=V/I_{\mathrm{G}}.\end{array}$$

3.4 Daily volume changes at lake A and water fluxes

The initial lake volume (4.7×10⁶ m³) was estimated from the observed lake surface area (2.79×10⁵ m²) and the maximal depth (20 m) and assuming bank slopes of 25^∘. Assuming bank slopes of 20 or 30^∘, a typical range for saturated sands (Holtz and Kovacs, 1981), would result in an estimated initial lake volume of 4.84×10⁶ m³ (+3 %) and 4.32×10⁶ m³ (−8 %). Lake A volume variations are estimated from daily water level changes and assuming a constant lake area. As water level measurements are only available for a short period at lake A, water levels at lake DM and observation well VP are used as proxies. Water levels at observation well VP were used as a proxy from 24 August to 30 October 2017, while water levels at lake DM were assumed representative of lake A for the rest of the study period (i.e. from 9 February to 23 August 2017 and from 31 October 2017 to 25 January 2018). This approximation is deemed acceptable because the simulation of δ_L depends on the remaining fraction of lake water f (not the absolute water level), and daily variations in the water levels at lake A, lake DM, and observation well VP were shown to be similar (see Sect. 2.2).

Evaporative fluxes (E) are calculated using the standardized Penman-48 evaporation equation, as described in Valiantzas (2006), as follows:

$$\begin{array}{}\text{(5)}& E_{\mathrm{Penman}\text{-}\mathrm{48}}={\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}\cdot {\displaystyle \frac{R_{\mathrm{n}}}{\mathit{\lambda }}}+{\displaystyle \frac{\mathit{\gamma }}{\mathrm{\Delta }+\mathit{\gamma }}}\cdot {\displaystyle \frac{\mathrm{6.43}f\left(u\right)D}{\mathit{\lambda }}},\end{array}$$

where R_n is the net solar radiation (megajoules per square metre per day), Δ is the slope of the saturation vapour pressure curve (kilopascals per degree Celsius; hereafter kPa ^∘C⁻¹), γ is the psychrometric coefficient (kPa ^∘C⁻¹), λ is the latent heat of vaporization (megajoules per kilogram), f(u) is the wind function (see Appendix B), and D is the vapour pressure deficit. For comparative purposes, estimation of the daily evaporative fluxes was also conducted with the Linacre-OW equation (Linacre, 1977) and the simplified Penman-48 equation (Valiantzas, 2006). These methods yielded similar evaporation estimates from April to August but underestimated total evaporation by 24 % to 33 % compared to the standardized Penman-48 equation. The discrepancy between the models is restricted to late summer and autumn (see Appendix C; Fig. C1) and is attributed to the difference between the air and water surface temperature, which was estimated based on the equilibrium method as described by de Bruin (1982; see Appendix D). Note that E and P are set to zero during the ice cover period (i.e. from 1 January to 31 March, based on meteorological data and field observations).

For well-mixed conditions, the ${\mathit{\delta }}_{Q_{\mathrm{s}}}$ and ${\mathit{\delta }}_{Q_{\mathrm{g}}}$ are assumed to be equal to δ_L. Hence, no separation of these two fluxes is attempted, and they are merged into one variable, i.e. the outflow (Q). The direction and intensity of the water flux at the lake–aquifer interface can be conceptually described by Darcy's law, which states that Q=KAi, where K is the hydraulic conductivity, A is the cross-sectional area through which the water flows, and i is the hydraulic gradient. Given the significant depth of lake A (i.e. 20 m) in comparison to the maximum water level change during the flooding event (i.e. 2.7 m), the variations in A and K are expected to have a minor impact on Q. Hence, the change in outflows from the lake is expected to be mainly controlled by i changes and, consequently, to be roughly proportional to the change in lake water level. Considering the above, it was assumed that the daily outflow flux from lake A varied linearly according to the lake water level; the minimum and maximum outflow (Q_min and Q_max) correspond to the minimum and maximum water level, respectively. The outflow range (i.e. minimum and maximum values) was adjusted to obtain the best fit between the calculated and observed δ_L.

Total daily inflow (sum of daily P, I_S, and I_G) into lake A compensates for the adjusted daily outflow and daily lake volume difference. The precipitation (P) is evaluated from the available meteorological data (see Sect. 3.1), while direct measurement of I_S and I_G was not possible in this hydrogeological context (see Sect. 2.1). Consequently, further assumptions are needed to apportion these contributions. Considering the proposed conceptual model of the groundwater–surface water interactions (see Sect. 2.2), I_S is set to zero, while I_G is contributing to the lake during the groundwater control period. On the other hand, during the floodwater control period (i.e. from 23 February to 8 May 2017), it is assumed that the rising water level at lake A results in a hydraulic gradient forcing the lake water to infiltrate into the aquifer, inhibiting I_G.

4 Results

From 23 February to 8 May 2017, the net water fluxes are mainly positive, and an overall volume increase is observed at lake A. The maximum volume change of lake A was 7.6×10⁵ m³, which represents 16 % of the lake's initial volume. The maximum net water flux was 1.2×10⁵ m³ d⁻¹, corresponding to a water level rise of 0.43 m (on 5 April 2017 only). From 9 May to mid-August 2017, lake A volume was decreasing, and the daily net water fluxes were mainly negative. In early August 2017, lake A regained its initial volume. Then, in autumn and winter, the volume of lake A was oscillating, and the net water fluxes were ranging from $-\mathrm{6.4}\times {\mathrm{10}}^{\mathrm{4}}$ to 5.3×10⁴ m³ d⁻¹. At the end of the study period (i.e. on 25 January 2018), a net volume difference of 1.5×10⁵ m³ remained at lake A compared to 9 February 2017.

However, the evolution of lake A volume and the net water fluxes are not representative of the surface water–groundwater interactions. Indeed, gross water fluxes are likely to exceed net water fluxes at natural and dredged lakes sitting in permeable sediments (Zimmermann, 1979; Arnoux et al., 2017a; Jones et al., 2016). In the context of this study, we conceptualized two main hydrological periods during which the lake water can either drain towards lake DM or exit the lake as groundwater output. To balance out these outputs, the inflows to lake A must therefore be greater than the net water fluxes.

For that reason, the development of a volume-dependent transient stable isotope mass balance was required to correctly depict the importance of the floodwater inputs on the water mass balance of the lake.

4.1 Isotopic and geochemical framework

The isotopic composition of precipitation (δ_P), lake A, and floodwater are depicted in Fig. 4. The local meteoric water line (LMWL) was defined using an ordinary least squares regression (Hughes and Crawford, 2012) using isotope data in precipitation from the Saint-Bruno station RIGR (Research Infrastructure on Groundwater Recharge) database (n=27; from December 2015 to June 2017).

Figure 4 Isotopic composition of precipitation, lake A water, and floodwater from March 2017 to January 2018. Hollow and solid blue circles correspond to samples collected at ≤2 and >2 m depth, respectively. Analytical precision is 0.15 ‰ and 1 ‰ at 1σ for δ¹⁸O and δ²H. Precipitation data are retrieved from the research infrastructure on groundwater recharge database (Barbecot et al., 2019).

For the study period, the isotopic composition of bulk precipitation was available on a biweekly to monthly time step (n=15) and ranged from −19.19 ‰ to −6.85 ‰ for δ¹⁸O and −144 ‰ to −38 ‰ for δ²H. Interpolation was used to simulate the δ_P on a daily time step for the isotope mass balance model computation.

Isotopic compositions of lake A water samples (n=39) are linearly correlated (see solid blue line) and all plotted below the LMWL, which confirms that lake A is influenced by evaporation. Linear regression of lake A water samples defines the local evaporation line (LEL), which is ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{5.68}(\pm \mathrm{0.27})\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{12.80}(\pm \mathrm{2.83})$ (R²=0.92). Some samples from the surface of lake A are plotted below the LEL, likely indicating snowmelt water inputs as noted in previous studies of Canadian lakes (Wolfe et al., 2007).

The isotopic composition of the floodwater samples (n=3) is indeed more depleted than lake A waters (i.e. δ¹⁸O from −11.85 ‰ to −11.18 ‰ and δ²H from −81 ‰ to −78 ‰) and is most likely reflecting the significant contribution from heavy isotope-depleted snowmelt waters. The floodwater samples are also linearly correlated and are plotted along a line (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{5.33}{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{18.82}$) for which the slope is similar to lake A LEL, suggesting that the sampled floodwater evaporated under the same conditions as lake A water samples. For simplification purposes, the isotopic composition of the surface water inflow (${\mathit{\delta }}_{I_{\mathrm{s}}}$) was set to the intersection between the floodwater LEL and the LMWL (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{12.00}$ ‰ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{83}$ ‰). The long-term (1997–2008) average and minimum and maximum isotopic signatures of Ottawa River water at Carillon (∼34 km upstream from lake DM; see Fig. 1b) for the month of April are −11.19 ‰, −12.01 ‰ and −10.23 ‰ for δ¹⁸O and −81 ‰, −85 ‰ and −77 ‰ for δ²H, respectively (Rosa et al., 2016). The mean and minimum values compare well with the observed isotopic signatures at lake DM during springtime in 2017.

The isotopic composition of groundwater (δ_G) can be determined from direct groundwater samples or indirectly from the amount-weighted mean δ_P. However, in highly seasonal climates, there is a widespread cold season bias to groundwater recharge (Jasechko et al., 2017), and estimating δ_G via groundwater samples or an amount-weighted mean δ_P may be misleading. In fact, it has been argued that the LMWL–LEL intersection better represents the isotopic composition of the inflowing water to a lake and is, thus, commonly used to depict the δ_G in isotopic mass balance applications (Gibson et al., 1993; Wolfe et al., 2007; Edwards et al., 2004). Concerning the study site, the estimated δ_G is −11.26 ‰ for δ¹⁸O and −77 ‰ for δ²H (i.e. the Saint-Bruno LMWL and lake A LEL intersection). The latter compares well with the mean isotopic signature of groundwater at Saint-Télesphore station (−11.1 ‰ for δ¹⁸O and −78.5 ‰ for δ²H; Barbecot et al., 2018) and is more depleted than the long-term amount-weighted mean δ_P at Ottawa (−10.9 ‰ for δ¹⁸O and −75 ‰ for δ²H; IAEA/WMO, 2018). Note that the locations of Saint-Télesphore station and Ottawa are depicted in Fig. 1b.

The geochemical facies of lake A and lake DM samples are illustrated in Fig. 5 by the means of a Piper diagram. Mean values for lake B and regional groundwater (GW) geochemical facies are also plotted for comparison purposes. Both lake A and floodwater were found to be Ca-HCO₃ types, which is typical for precipitation- and snowmelt-dominated waters (Clark, 2015). The geochemistry of lake A is relatively constant throughout the year and reveals a depth-wise homogeneity. The geochemistry of lake B is distinct from lake A and appears to be influenced by regional groundwater characterized by a Na-Cl water type.

Figure 5 Geochemical facies of lake A (n=23) and floodwater (n=1). Mean values for lake B (n=42) and regional groundwater (GW) (n=11) geochemical facies are also plotted. Lake A and floodwater are characterized by Ca-HCO₃ water types, while lake B and regional GW correspond to Na-Cl water types. Note that regional GW was sampled upstream of lake B.

4.2 Evaluation of the water budget

4.2.1 Volume dependent isotopic mass balance model

As described in Sect. 3.3, the isotopic mass balance model was solved iteratively by recalculating δ_L on a daily time step. This model was developed assuming (1) well-mixed conditions and (2) that the outflow flux changes are roughly proportional to the lake's water level changes. We adjusted the minimum and maximum outflow fluxes (Q_min and Q_max) so that they corresponded to the minimum and maximum water levels (see Fig. 3).

In total, three sampling campaigns (i.e. on 9 February 2017, 17 August 2017, and 25 January 2018) were conducted at lake A in order to collect water samples for isotopic analyses from the epilimnion, metalimnion, and hypolimnion (Fig. 6; Appendix E; Fig. E1) to account for the vertical stratification of the isotopic signature (Gibson et al., 2017). The vertical isotopic profiles were volume-weighted according to the representative layer for each discrete measurement in order to obtain the observed δ_L for each campaign (Table 1). The depth-averaged isotopic composition of the lake on 9 February 2017 (i.e. ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{10.15}$ ‰ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{70}$ ‰) was used as the initial modelled δ_L.

Figure 6 Observed and modelled depth-average isotopic composition of the lake (δ_L) for δ¹⁸O (a, b) and δ²H (c, d) from 9 February 2017 to 25 January 2018 for scenarios A and B. The hollow and solid blue circles correspond to lake A water samples collected at ≤2 and >2 m, respectively. The modelled δ_L is fitted against the three depth-averaged δ_L and an additional sample collected at ≤2 m depth on 9–10 May 2017 (scenario A) and 27 April 2017 (scenario B). These samples are marked by the hollow red squares. The grey shaded area corresponds to the floodwater control period. The error bars correspond to the standard error in the samples for each campaign. The dashed lines represent the modelled δ_L when considering that 25 % to 100 % of the outputs from the lake during the floodwater control period were temporally stored in the aquifer and discharged to the lake as floodwater-like inputs (${\mathit{\delta }}_{I_{\mathrm{s}}}$).

Table 1 Observed depth averaged (or mean) and standard deviation (SD) of isotopic composition of lake A for the sampling campaigns in February 2017, August 2017, and January 2018 and all samples. The isotopic composition of the samples collected at the surface of lake A on 9–10 May and 27 April 2017 are also listed. The asterisks (^*) indicate that a mean value was calculated (instead of a depth-averaged value).

While depth-averaged δ_L was not available during the floodwater control period (i.e. late February to early May), water samples from the surface of lake A provide relevant evidence to better constrain the model. It is likely that lake A was fully mixed during the floodwater control period, and that the water samples collected at the surface of lake A on 27 April or 9–10 May 2017 are representative of the whole water body. Indeed, the observed surface water temperature was <5 ^∘C until early May (see Fig. C1) and suggests a limited density gradient along the water column which does not allow for the development of thermal stratification. In this context, we opted to simulate two scenarios (A and B) for which the isotopic mass balance model is either constrained at ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{11.20}$ ‰ and${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{76}$ ‰ on 9–10 May 2017 or at ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{11.86}$ ‰ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{80.68}$ ‰ on 27 April 2017.

The results of the volume-dependent isotopic mass balance for δ¹⁸O and δ²H are illustrated in Fig. 6. The fitted Q_min and Q_max from lake A are 3.7×10⁴ and 8.0×10⁴ m³ d⁻¹ for scenario A and 1.0×10³ and 2.8×10⁵ m³ d⁻¹ for scenario B. These first-order water flux estimates represent equivalent water level variations ranging from 0.004 and 1.0 m d⁻¹. From 23 February to 8 May 2017 (see the grey shaded area), hydraulic conditions allowed for surface inputs (I_S) from lake DM to lake A at a mean rate of 6.61×10⁴ m³ d⁻¹, with a total floodwater volume of 4.82×10⁶ m³ for scenario A. The total floodwater volume was twice as important (9.96×10⁶ m³) for scenario B. Then, from 9 May 2017, we considered that these floodwater inputs stopped as the lake water level started to decrease. As a consequence, the model yielded a gradual enrichment of δ_L due to the combined contribution from I_G and E for both scenarios. From 9 May 2017 to 25 January 2018, the total I_G were 1.16×10⁷ and 1.48×10⁷ m³ for scenario A and B, respectively. Overall, the δ¹⁸O and δ²H models were better at reproducing the January 2018 and August 2017 observed δ_L, respectively. This is likely linked to the uncertainties and representativeness of the meteorological data, which is controlling the isotopic fractionation due to evaporation.

While the computed flows for scenario A are within a plausible range for the combination of surface and groundwater outflow processes (i.e. minimum and maximum equivalent water level variations of 0.13 and 0.29 m d⁻¹), scenario B yielded less realistic results (i.e. minimum and maximum equivalent water level variations of 0.004 and 1.0 m d⁻¹). As mentioned above, scenario B was constrained at ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{11.86}$ ‰ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{80.68}$ ‰ in late April (Fig. 6), based on a surface water sample which was taken during a temporarily decreasing water level period (Fig. 3) and is, thus, likely less representative of the overall lake's dynamics compared to scenario A. This demonstrates the limit of the approach and shows that it is important to correctly constrain the model during flood events in order to perform precise estimations of the water balance.

The water mass balance of lake A from 9 February 2017 to 23 January 2018 is summarized in Table 2 for both scenarios. The difference between the total inputs and total outputs correspond to the lake volume difference (1.48×10⁵ m³) between the start and the end of the model run. Groundwater inputs (I_G) and surface water inputs (I_S) account for 71 % and 28 % of the total water inputs to the lake for scenario A. While I_S are twice as important for scenario B, they only account for 39 % (+11 %) of the total inputs, and the I_G are 60 % (−11 %). It thus appears that the annual dynamic of lake A is dominated by groundwater inputs for both scenarios, despite the intensity of the flood event. For scenarios A and B, t_f, as defined in Eq. (4), is similar (i.e. 135 and 110 d). Precipitation is contributing 1 % of the total annual inputs and evaporation only accounts for 2 % of the total annual outputs. Although the establishment of a hydraulic connection between lake DM and lake A is a recurring yearly hydrological process, it is important to note that the magnitude and duration of the flooding event of 2017 was particularly important and, thus, had a greater impact on the dynamic of lake A in comparison to other years.

Table 2 Water mass balance of lake A for scenarios A and B. The difference between the total inputs and total outputs corresponds to the lake volume difference over the study period. The total inputs (I) correspond to the sum of precipitation (P), surface water inflow (I_S), and groundwater inflow (I_G). The total outputs (Q) correspond to the sum of evaporation (E) and surface water and groundwater outflow (Q). The mean flushing time (t_f) is the ratio of the lake volume to the mean groundwater inputs (I_G).

4.2.2 Sensitivity analysis

A one-at-a-time (OAT) sensitivity analysis was performed to grasp the relative impact of the input parameters' uncertainties on the model outputs. For each parameter, we tested two scenarios which delimit the uncertainty for each parameter. First, we tested the sensitivity of the model for V+3 % and V−8 % (i.e. estimated with slopes of 30 and 20^∘). Concerning ${\mathit{\delta }}_{I_{\mathrm{s}}}$ and δ_G, the model was tested for ±0.5 ‰ for δ¹⁸O and ±4 ‰ for δ²H, assuming they would both evolve along the LMWL (see Fig. 3). Then, we assessed the sensitivity of the model to δ_A by fixing the seasonality factor k at 0.5 and 0.9. Evaporation was computed with ±20 %, whereas the meteorological parameters (i.e. relative humidity (RH), T_air, U, P, and Rs) were tested for ±10 %. As E and δ_A are dependent on the water surface temperature, we also tested the sensitivity of the model when considering that T is equal to the daily mean air temperature (T_air). Finally, we tested for the uncertainties concerning the definition of the LMWL. For the reference scenario, the LMWL (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{8.13}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{14.78}$) was estimated using an ordinary least square regression (OLSR). For the sensitivity analysis, we estimated the LMWL via a precipitation amount weighted least square regression (PWLSR), which was developed by Hughes and Crawford (2012). Using the PWLSR method, the LMWL is defined as ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{8.28}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}+\mathrm{17.73}$, and ${\mathit{\delta }}_{I_{\mathrm{s}}}$ and δ_G are estimated at −12.39 ‰ and −11.74 ‰ for δ¹⁸O and at −85 ‰ and −79 ‰ for δ²H, respectively. Recalculation of ${\mathit{\delta }}_{I_{\mathrm{s}}}$ and δ_G was needed, as they were both assumed to be plotted on the LMWL (see Sect. 4.1).

The results of this sensitivity analysis are listed in Tables F1 and F2 (Appendix F) for scenarios A and B. Overall, the model was found to be highly sensitive to the uncertainties associated with ${\mathit{\delta }}_{I_{\mathrm{s}}}$, δ_G, and E, as the annual mean water fluxes (Q and I) varied up to −31 % and +46 % compared to the reference scenarios A and B. A negligible to slight change in the modelled δ_L was found when considering the uncertainties for V, δ_A, RH, T_air, U, P, and Rs. For these variables, the mean flux estimate (Q and I) changes ranged from −8 % to +4 % compared to the reference scenarios A and B. As expected, the value of ${\mathit{\delta }}_{I_{\mathrm{s}}}$ affects the modelled δ_L exclusively during the floodwater control period. Similarly, the values of δ_G and E particularly influence the modelled δ_L from late summer to early winter. This is due to the fact that Q and E are the dominant fluxes during this period. When considering that T is equal to T_air, despite the significantly different maximum and minimum values for Q, the mean Q was relatively similar to the reference scenarios, and only a small change for t_f (+3 % and +2 % compared to reference scenarios A and B) was found. Finally, the model is highly sensitive to the uncertainties associated with the LMWL, as a translation of the LMWL implies an enrichment or depletion of both the ${\mathit{\delta }}_{I_{\mathrm{s}}}$, δ_G at the same time. Such modifications result in mean flux estimate (Q and I) changes of up to −38 % and −43 % compared to reference scenarios A and B.

4.3 Importance of temporary floodwater storage on the water balance partition

The developed isotopic mass balance model yielded significant floodwater inputs during springtime to be a best fit to the observed δ_L. A first-order estimate of the total floodwater volume summed to 4.82×10⁶ m³ (for scenario A), which is nearly equal to the lake's initial volume (i.e. 4.70×10⁶ m³). Similar results were obtained by Falcone (2007), who studied the hydrological processes influencing the water balance of lakes in the Peace–Athabasca Delta, Alberta (Canada) using water isotope tracers. They reported that a springtime freshet (in 2003) did replenish the flooded lakes from 68 % to >100 % (88 % on average).

As mentioned in Sect. 2.3, it was conceptualized that the high surface water elevation of lake A during springtime resulted in hydraulic gradients that forced lake water to infiltrate into the aquifer and induce local recharge (see Fig. 3). An important volume of flood-derived water could, thus, be stored during the increasing water level period and eventually discharged back to the lake as its water level decreased. Hence, the groundwater inputs to lake A following the flooding event likely corresponded to flood-derived surface water originating from lake DM. Considering that these fluxes are characterized by a floodwater-like isotopic signature (${\mathit{\delta }}_{I_{\mathrm{s}}}$), rather than the isotopic signature of groundwater (δ_G), the temporal evolution of the modelled δ_L would be modified. Such a consideration is noteworthy for a better depiction of the importance of floodwater inputs in the water balance partition.

Assuming that 25 % to 100 % of the outputs (Q) from the lake during the floodwater control period were temporally stored in the aquifer and did eventually discharge back to the lake, the modelled δ_L diverges more or less from the reference scenarios A and B (Fig. 6; see the dashed lines). It is noteworthy that a better fit between the modelled δ_Land depth-averaged δ_L is obtained when considering that 25 % to 50 % of the outputs (Q) from the lake during the floodwater control period discharge back to the lake. In fact, it is likely that part of the potential stored floodwater could have effectively discharged back to the lake. For instance, part of the floodwater-like groundwater could have been abstracted by the pumping wells at the adjacent bank filtration site or could have been discharged to lake B. These results illustrate the importance of considering temporary subsurface floodwater storage when assessing water balances, especially as the magnitude and frequency of floods are likely to be more important in the future (Aissia et al., 2012).

4.4 Temporal variability in the water balance partition

The water balance presented in Table 2 provides an overview of the relative importance of the hydrological processes at lake A for the study period (i.e. February 2017 to January 2018). As the surface water inputs (as floodwater) only occurred during springtime at lake A, it is also important to decipher the temporal variability in the water fluxes. The dependence of a lake on groundwater can be quantified via the G index, which is the ratio of cumulative groundwater inputs to the cumulative total inputs (Isokangas et al., 2015). Figure 7 shows the temporal evolution of the G index from 9 February 2017 to 25 January 2018 for scenario A and the associated scenarios (A1 to A22) considered in the sensitivity analysis. Note that the G index is calculated at a daily time step, based on the cumulative water fluxes. It is used to understand the relative importance of groundwater inputs over the studied period and does not consider the initial state of the lake. In early February, the G index was 100 % because no surface water inputs (I_S) or precipitation (P) had yet contributed to the water balance. During the floodwater control period (see grey shaded area), the G index rapidly decreased and reached 12 % on 8 May 2017 (for the reference scenario A). A gradual increase in the G index is then computed for the rest of the study period. On 25 January 2018, the G index is 71 % and is likely more representative of annual conditions. Despite the sensitivity of the model to the input parameters, all scenarios yielded similar results. The G index ranged from 62 % to 75 % on an annual timescale for the different scenarios.

Figure 7 Temporal evolution of the G index from 9 February 2017 to 25 January 2018 for scenario A and the associated scenarios considered in the sensitivity analysis (i.e. A1 to A22). The grey shaded area corresponds to the floodwater control period. The dashed lines correspond to the G index when considering that 25 % to 100 % of the outputs from the lake during the floodwater control period were temporally stored in the aquifer and discharged to the lake as floodwater-like inputs (I_S, ${\mathit{\delta }}_{I_{\mathrm{s}}}$).

The impact of the potential temporary subsurface storage is also depicted in Fig. 7 (see dashed lines). As highlighted in Sect. 4.3, part of the potentially stored floodwater in the aquifer could have discharged back to the lake as floodwater-like inputs (${\mathit{\delta }}_{I_{\mathrm{s}}}$) after the flooding event. Considering these fluxes as surface water inputs (I_S), rather than groundwater inputs (I_G), would alter the temporal evolution of the G index. Assuming that 25 % to 100 % of the outputs from the lake during the floodwater control period did eventually discharge back to the lake, the floodwater inputs would contribute to the lake water balance until early June to early August (Fig. 7). Lake A would, thus, be dependent on flood-derived water during a 1- to 3-month period after the flooding event. On an annual timescale, the temporary subsurface storage could lower the G index to a minimum value of 47 %.

5 Discussion

5.1 Resilience of lakes to groundwater pollution

The resilience of a system has been defined as its capacity to cope with perturbations (i.e. internal and/or external changes) while maintaining its state (Cumming et al., 2005). In the case of a lake, perturbations can manifest as a change in the water quantity and quality contributing to the water balance. According to Arnoux et al. (2017b), the impact of a perturbation to a lake is not only dependent on the relative importance of water budget fluxes but also on the residence time of water in the lake. Thus, they proposed an interpretation framework which relates the response time of a lake to changes in groundwater quantity and/or quality, thereby linking the G index with t_f, which is the mean flushing time by groundwater fluxes (Fig. 8). They depict a general case, applicable to any pollution, regardless of reactivity or the fate of contaminants. Hence, care should be taken when interpreting the sensitivity to specific contaminants which are subject to attenuation processes, such as degradation and sorption.

Figure 8 Resilience of lakes to groundwater quantity and quality changes for lake A (this study) and kettle lakes (Arnoux et al., 2017b) in southern Quebec (Canada). The G index is the ratio of groundwater inputs to total inputs, and t_f is the mean flushing time by groundwater. This representation is adapted from Arnoux et al. (2017b).

In their study, Arnoux et al. (2017b) assessed the resilience of kettle lakes (n=20), located in southern Quebec (Canada), in similar morpho-climatic contexts to lake A. The surveyed lakes were found to be characterized by a wide range of conditions, from resilient (i.e. G index < 50 % and t_f>5 year) to highly sensitive to groundwater changes (i.e. G index > 50 % and t_f<1 year). This is related to the variability in the hydrogeological contexts, resulting in variations in the importance of groundwater contributions and the range of mean flushing times of the lakes (see the grey arrow in Fig. 8). The majority of the lakes (i.e. 50 %) were found to be characterized by intermediate conditions (G index > 50 % and $\mathrm{5}<t_{\mathrm{f}}<\mathrm{1}$ year) and, thus, were classified as being relatively resilient to both surface and groundwater changes.

Concerning lake A, studied scenarios (i.e. reference scenario A and the sensitivity analysis) yielded values for G index > 50 % and t_f<1 year, i.e. indicating that lake A is highly sensitive to groundwater changes but resilient to surface pollution. Nevertheless, it was shown that temporary floodwater storage and discharge to lakes are crucial to correctly represent the G index by accounting for the origin of water fluxes (Fig. 7; Sect. 4.4). While floodwater storage lowers the G index, the t_f slightly increases (see orange arrow in Fig. 8). Therefore, the studied lake receives a reduced groundwater contribution relative to the initial estimated apportionment when not accounting for floodwater storage but is still characterized by a rapid flushing time. This implies that flood-affected lakes are more likely to be characterized by an intermediate condition and, thus, are relatively resilient to groundwater quantity and quality changes. The geochemical data (Sect. 4.2) are in agreement with this interpretation. Indeed, a low mineralization and Ca-HCO₃ water type at lake A is consistent with the significant floodwater contributions (to the lake and aquifer). In comparison, the neighbouring lake (i.e. lake B) does not undergo yearly recurrent flooding and was shown to be more mineralized with a Na-Cl water type, likely originating from the road salt contamination of regional groundwater (Pazouki et al., 2016). Biehler et al. (2020) similarly reported hydrological controls on the geochemistry of a shallow aquifer in an hyporheic zone, where the river stage influenced the mixing ratio between river water and the deeper aquifer.

Considering the above, it is possible to speculate about the potential future impacts of climate change on lake A. Globally, future meteorological scenarios are predicting changes in precipitation and climate extremes, including floods and droughts (Salinger, 2005). In Quebec (Canada), river stages are expected to increase across various watersheds in response to future climate scenarios (Roy et al., 2001; Dibike and Coulibaly, 2005; Minville et al., 2008). These hydrological responses could result in floods of longer duration and higher intensity (Aissia et al., 2012) and more pronounced droughts (Wheaton et al., 2007). Such changes could directly affect the quality of lake A. If flooding becomes more prevalent, enhanced floodwater input to lake A would likely occur. In this case, the surface water inputs from floods would buffer the sensitivity of lake A to groundwater quality changes originating from its watershed. On the other hand, if floods become less important and/or less frequent, we can expect that the water quality of lake A would be more dependent on regional groundwater quality. In such a case, the geochemistry of lake A could potentially shift towards that of lake B, and an increase in the salinity and the concentration of Na⁺, Ca²⁺, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}+}$, and Cl⁻ would be expected for lake A.

5.2 Implications for water management

Water budget assessments at natural lakes can serve as a tool for quantifying local human impacts (i.e. land use changes and climate changes) on water resources (Arnoux et al., 2017b). Based on the results of this study, it becomes apparent that water budget assessments at artificial lakes (such as lake A) can also be used to track human impacts on water resources. Recurring water budget assessments at a specific lake over time will serve to document changes in groundwater and surface water apportionment and can help to detect changes in local groundwater availability, and to anticipate impacts on local water supply utilities. As the response time of a lake to changes is controlled by its flushing time, the evolution of the G index will manifest at various rates. Indeed, lakes with different t_f would reflect changes at different timescales. For instance, lakes with t_f>5 year would be expected to respond to decadal changes, while lakes with t_f<5 year would track annual or interannual variability. By analogy, we might postulate that it would be informative to study lakes with rapid response times (i.e. t_f<1 year), as they will act as precursors of the evolution of nearby surface water bodies characterized by longer flushing times.

As demonstrated, isotopic approaches may be efficiently employed to solve water budget unknowns as the method can be performed at low cost and requires limited sampling and monitoring efforts for flood-affected environments which may be difficult or dangerous to monitor using traditional approaches. To enhance the effectiveness of our approach, the sampling strategy may be improved. First, surface water sampling for isotopic analyses is recommended during turnover periods (i.e. springtime and autumn) and should be combined with depth-resolved measurements of physico-chemical parameters to confirm the vertical homogeneity or stratification. Second, for long-duration flood events, monitoring of potential evolution in floodwater isotopic signatures could help to improve the accuracy and realism of the model. Groundwater-level monitoring and groundwater sampling in the vicinity of the lake could also help to strengthen the conceptual model by providing data to interpret the direction of groundwater fluxes and the variability in isotopic composition through time.

6 Conclusions

In this study, a volume-dependent transient isotopic mass balance model was developed and applied to a flood-affected lake in an ungauged basin in southern Quebec (Canada). This allowed for a better understanding of the resilience of a flood-affected lake to changes in the surface–groundwater water balance partition, for an understanding of the role of floodwater, and for predicting the resilience to groundwater quantity and quality changes for a local water supply. Given the contrasting isotopic signature of the floodwater, the isotopic mass balance model was effectively applied at the study site. We anticipate that the isotopic framework is likely to be transferable to other lake systems subject to periodic flooding, including lowland lakes fed by mountain floodwaters, river deltas, wadis, or nival (snowmelt-dominated) regimes, the latter of which dominates the high-latitude and high-altitude cold regions, including much of the Canadian landmass.

The isotopic mass balance model revealed that groundwater inputs dominated the annual water budget. To test the sensitivity, representativeness, and resilience of the model, several model scenarios were evaluated to account for uncertainty in important input variables. Despite sensitivity to some variables, all model scenarios considered in the sensitivity analysis converged on the results that lake A is mainly dependent on groundwater inputs and has a rapid (<1 year) flushing time by groundwater, suggesting that lake A would be highly sensitive to groundwater quantity and quality changes.

When taking into account potential subsurface storage, a better fit could be obtained between the modelled and depth-averaged isotopic signature of the lake, suggesting that the contribution of floodwater-like subsurface inputs is important to consider when assessing for water balance at flood-affected lakes. In fact, the increased contribution of surface water (from subsurface storage) resulted in a lower contribution from groundwater and, consequently, in an increased resilience to groundwater changes. This finding provides a basis for postulating the impact of climate change on the water quality of lake A. If the importance of floods increases, more floodwater inputs to lake A can be expected during springtime, causing increased recharge. In this case, the surface water inputs from floods would increase the resilience of flood-affected lakes to groundwater quantity and quality changes at the watershed scale. On the other hand, if floods become less severe and/or less frequent, we can expect that the water quality of flood-affected lakes becomes more dependent on regional groundwater quality. From a global perspective, performing water balance assessments at lakes with rapid flushing time (<1 year) can help to predict the evolution of other surface water bodies with longer flushing times in their vicinity and, therefore, is useful for establishing regional-scale management strategies for maintaining lake water quality.

Appendix A

Table A1 Detailed source information and download links for the openly available geospatial data in Fig. 1. All data are openly available and are used in accordance with the Open Government Licence – Canada or the Open Data Policy, M-13-13, of the United States Census Bureau.

Figure A1 Geological cross section along the pumping wells showing the buried valley carved into the Champlain Sea clays and filled with alluvial gravels and sands.

Appendix B

Computation of isotope mass balance parameters

The parameter f(u), for the estimation of E (Eq. 5), is calculated according to the area-dependent expression described by McJannet et al. (2012) as follows:

$$\begin{array}{}\text{(B1)}& f\left(u\right)=(\mathrm{2.36}+\mathrm{1.67}u)A^{-\mathrm{0.05}},\end{array}$$

where u is the wind speed (metres per second) measured at 2 m above the ground, and A is the area (square metres) of the lake. Note that Eq. (6) was developed for land-based meteorological data.

The isotopic composition of the evaporating moisture (δ_E) is estimated based on the Craig and Gordon (1965) model and, as described by Gonfiantini (1986), is as follows:

$$\begin{array}{}\text{(B2)}& {\mathit{\delta }}_{\mathrm{E}}={\displaystyle \frac{\frac{\left({\mathit{\delta }}_{\mathrm{L}}-{\mathit{\epsilon }}^{+}\right)}{{\mathit{\alpha }}^{+}}-h{\mathit{\delta }}_{\mathrm{A}}-{\mathit{\epsilon }}_{\mathrm{K}}}{\mathrm{1}-h+{\mathrm{10}}^{-\mathrm{3}}{\mathit{\epsilon }}_{\mathrm{K}}}}\left(\mathrm{‰}\right),\end{array}$$

where h is the relative humidity normalized to water surface temperature (in decimal fraction), δ_A is the isotopic composition of atmospheric moisture (described later on), ε⁺ is the equilibrium isotopic separation, and ε_K is the kinetic isotopic separation, with ${\mathit{\epsilon }}^{+}=({\mathit{\alpha }}^{+}-\mathrm{1}){\mathrm{10}}^{\mathrm{3}}$ and ${\mathit{\epsilon }}_{\mathrm{K}}=\mathit{\theta }\cdot C_{\mathrm{K}}$ (1 h). α⁺ is the equilibrium isotopic fractionation, θ is a transport resistance parameter, and C_K is the ratio of molecular diffusivities of the heavy and light molecules. θ is expected to be close to 1 for small lakes (Gibson et al., 2015), and C_K is typically fixed at 14.2 ‰ and 12.5 ‰ for δ¹⁸O and δ²H, respectively, in lake studies as these values represent fully turbulent wind conditions (Horita et al., 2008). Experimental values for α⁺ were used (Horita and Wesolowski, 1994) as follows:

$$\begin{array}{}\text{(B3)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}^{+}\left({}^{\mathrm{18}}\mathrm{O}\right)& =\exp \left[-{\displaystyle \frac{\mathrm{7.685}}{{\mathrm{10}}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{6.7123}}{(T+\mathrm{273.15})}}\right.\\ & \left.-{\displaystyle \frac{\mathrm{16666.4}}{(T+\mathrm{273.15}{)}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{350410}}{(T+\mathrm{273.15}{)}^{\mathrm{3}}}}\right]\end{array}\text{(B4)}& {\displaystyle }\begin{array}{rl}{\mathit{\alpha }}^{+}\left({}^{\mathrm{2}}\mathrm{H}\right)& =\exp \left[\mathrm{1158.8}\left({\displaystyle \frac{(T+\mathrm{273.15}{)}^{\mathrm{3}}}{{\mathrm{10}}^{\mathrm{12}}}}\right)\right.\\ & +\mathrm{1620.1}\left({\displaystyle \frac{(T+\mathrm{273.15}{)}^{\mathrm{2}}}{{\mathrm{10}}^{\mathrm{9}}}}\right)\\ & +\mathrm{794.84}\left({\displaystyle \frac{(T+\mathrm{273.15})}{{\mathrm{10}}^{\mathrm{6}}}}\right)\\ & \left.-{\displaystyle \frac{\mathrm{161.04}}{{\mathrm{10}}^{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{2}\mathrm{999}\mathrm{200}}{(T+\mathrm{273.15}{)}^{\mathrm{3}}}}\right],\end{array}\end{array}$$

where T is the water surface temperature (degrees Celsius), which was estimated according to the equilibrium method as described by de Bruin (1982) (see Appendix D).

The parameters m and δ_s, for the computation of δ_L (Eq. 3), are calculated as follows (Gibson, 2002):

$$\begin{array}{}\text{(B5)}& {\displaystyle }& {\displaystyle }m={\displaystyle \frac{\left(h-{\mathrm{10}}^{-\mathrm{3}}\times \left({\mathit{\epsilon }}_{\mathrm{K}}+\frac{{\mathit{\epsilon }}^{+}}{{\mathit{\alpha }}^{+}}\right)\right)}{\left(\mathrm{1}-h+{\mathrm{10}}^{-\mathrm{3}}\times {\mathit{\epsilon }}_{\mathrm{K}}\right)}}\text{(B6)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{S}}={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{I}}+mX{\mathit{\delta }}^{*}}{\mathrm{1}+mX}},\end{array}$$

where δ^* is the limiting isotopic composition that the lake would approach as V→0 and is calculated as follows:

$$\begin{array}{}\text{(B7)}& {\mathit{\delta }}^{*}=\left(h{\mathit{\delta }}_{\mathrm{A}}+{\mathit{\epsilon }}_{\mathrm{K}}+{\displaystyle \frac{{\mathit{\epsilon }}^{+}}{{\mathit{\alpha }}^{+}}}\right)/\left(h-{\mathrm{10}}^{-\mathrm{3}}\times \left({\mathit{\epsilon }}_{\mathrm{K}}+{\displaystyle \frac{{\mathit{\epsilon }}^{+}}{{\mathit{\alpha }}^{+}}}\right)\right).\end{array}$$

The isotopic composition of atmospheric moisture (δ_A) is estimated using the partial equilibrium model of Gibson et al. (2015) as follows:

$$\begin{array}{}\text{(B8)}& {\mathit{\delta }}_{\mathrm{A}}={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{P}}-k{\mathit{\epsilon }}^{+}}{\mathrm{1}+{\mathrm{10}}^{-\mathrm{3}}\times k{\mathit{\epsilon }}^{+}}},\end{array}$$

where δ_P is the isotopic composition of precipitation, and k is a seasonality factor fixed at 0.5 in this study. The k value (ranging from 0.5 to 1) is selected to provide a best fit between the measured and modelled local evaporation line. In Eq. (12), δ_P and monthly exchange parameters (ε⁺, α⁺, and ε_K) are evaporation flux weighted based on daily evaporation records.

Appendix C

Comparison of the evaporative flux (E) estimations

Figure C1 Cumulative evaporative fluxes from lake A via the standardized Penman-48, simplified Penman-48 (Valiantzas, 2006), and Linacre-OW (Linacre, 1977) equations.

Appendix D

D1 Estimation of the water surface temperature based on the equilibrium method (de Bruin, 1982)

The water surface temperature (T) was estimated via the equilibrium method presented by de Bruin (1982) because no continuous measurements were available. This model is based on the assumption of a well-mixed surface body and was developed from standard land-based weather data. It was tested on two adjacent reservoirs in the Netherlands, with average depths of 5 and 15 m, respectively. Similar to de Bruin (1982), we used the 10 d mean values because we are interested in the annual variations in the water temperature. Moreover, the 10 d mean values were found to better simulate the observed water surface temperature. The difference between the observed and modelled water temperature is typically ≤1 ^∘C, except in July and December where discrepancies of up to 5 ^∘C were observed (Fig. D1). This is likely because lake A develops a thermal stratification over summertime and in wintertime. Potential uncertainties in isotopic mass balance models due to stratification in lakes up to 35 m were previously described and discussed by Gibson et al. (2017, 2019). They reported that sampling methods and lake stratification can lead to volume-dependent bias in the water balance partition. In this study, not accounting fully for thermal stratification will lead to overestimation of evaporation fluxes, and groundwater exchange will potentially be underestimated.

Figure D1 Temporal evolution of air temperature and observed and estimated water surface temperatures at lake A. Water surface temperature estimations were computed according to the equilibrium method described by de Bruin (1982).

Appendix E

E1 Isotopic composition along vertical profiles in lake A

Figure E1 Isotopic composition of lake A water samples against depth on 9 February 2017, 17 August 2017, and 25 January 2018. The hollow circles and solid circles represent samples collected at ≤2 and >2 m depth, respectively.

Appendix F

F1 Results of the sensitivity analysis for reference scenarios A and B

See Tables F1 and F2.

Table F1 Sensitivity analysis on the input parameters of the isotopic mass balance model. Q is the output flux from lake A, I the input flux, and t_f the mean flushing time by groundwater.

Table F2 Sensitivity analysis on the input parameters of the isotopic mass balance model for the reference scenario B. Q is the output flux from lake A, I the input flux, and t_f the mean flushing time by groundwater.