# Effects of seasonal and interannual variability in water isotopes (δ2H, δ18O) on estimates of water balance in a chain of seven prairie lakes

H A Haig, N M Hayes, G L Simpson, Y Yi, B Wissel, K R Hodder, P R Leavitt

► To cite this version:

H A Haig, N M Hayes, G L Simpson, Y Yi, B Wissel, et al.. Effects of seasonal and interannual variability in water isotopes (δ2H, δ18O) on estimates of water balance in a chain of seven prairie lakes. Journal of Hydrology X, 2020, 10, 10.1016/j.hydroa.2020.100069. hal-04835329

HAL Id: hal-04835329

https://hal.science/hal-04835329v1

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research papers

# Effects of seasonal and interannual variability in water isotopes (δ²H, δ¹⁸O) on estimates of water balance in a chain of seven prairie lakes

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

H.A. Haig<sup>a, *,</sup> , N.M. Hayes<sup>a, b</sup> , G.L. Simpson<sup>c</sup> , Y. Yi<sup>d, e</sup> , B. Wissel<sup>c</sup> , K.R. Hodder<sup>f</sup> , P.R. Leavitt<sup>a, c, g</sup>

ARTICLE INFO

Keywords:

Water isotopes

Hydrology

Scales of variability

High resolution

Long-term monitoring

Water balance

Prairie lakes

Qu'Appelle

ABSTRACT

Stable isotopes of hydrogen (δ²H) and oxygen (δ¹⁸O) provide important quantitative measures of lake hydrology and water balance, particularly in lakes where monitoring of fluxes is incomplete. However, little is known of the relative effects of seasonal variation in water isotopes on estimates of lake hydrology, particularly over decadal scales. To address this gap, we measured water isotopes bi-weekly May-September during 2003–2016 in seven riverine lakes within the 52,000 km² Qu'Appelle River drainage basin of the Canadian Prairies. Analyses revealed that within-year variation in δ¹⁸O values routinely exceeded that among years, reflecting rapid changes in water source, particularly in lakes with water residence times <1 year. Isotopic variation was greatest during spring following snowmelt, except in large deep lakes which exhibited limited differences among seasons or years. In contrast, large hydrological events (e.g., 1-in-140-year flood in 2011) homogenized isotopic values, even among riverine lakes separated by over 150 km, and exerted particularly strong legacy effects on large lakes. Overall, study lakes exhibited a strongly positive moisture balance (evaporation < inflow), despite regional precipitation deficits of 30 cm yr⁻¹, with greater reliance on rainfall (vs. snow) and possibly evaporation in downstream lakes within more humid regions. We conclude that seasonal samples of water isotopes are required to characterize the hydrology of shallow lakes, or those with unknown reliance on snowmelt waters, as well as to better quantify lake susceptibility to climate variability.

# 1. Introduction

Better understanding of the vulnerability of surface waters to climate change is necessary to assess how water resources may change in the future. With one-sixth of the Earth's population residing in regions where water supplies are regulated mainly by snowmelt (Barnett et al., 2005), warming winter temperatures that reduce snowfall and its accumulation in cold regions are of particular concern (Akinremi et al., 1999; Coles et al., 2017). In addition, endorheic basins in sub-humid regions may exhibit the greatest sensitivity to climate change because they already show pronounced variations in water supply and lake level (Wang et al., 2018). However, data on surface water availability and fluxes are limited in many areas across the globe, largely due to the predominance of ungauged basins (Kirchner, 2006; Wood et al., 2011).

For example, only 12% of terrestrial land in Canada has sufficient instrumentation to evaluate lake and river vulnerability to changes in climate (Coulibaly et al., 2013). In areas where instrumental records are insufficient, alternative methods are needed to assess the hydrology of surface water bodies. With recent advances in analysis and modeling, measurements of stable isotopes of water (δ¹⁸O, δ²H) provide a promising means to quantify water balance, better assess controls of hydrologic variation, and evaluate surface water vulnerability to climate change (Bam and Ireson, 2019; Remmer et al., 2018; Seibert and McDonnell, 2002; Tunaley et al., 2017).

In principle, δ¹⁸O and δ²H of water (hereafter water isotopes) can be used as natural tracers to infer differences in water balance (evaporation to inflow ratios; E:I) and source waters (rain, snow, groundwater) in diverse continental surface waters (Gibson et al., 2019, 2016b, 2016a;

<footnote><sup>a</sup> Limnology Laboratory, Department of Biology, University of Regina, Regina, Saskatchewan S4S 0A2, Canada</footnote>

<footnote><sup>b</sup> College of Biological Sciences, University of Minnesota, Minneapolis-St. Paul, MN 55108, USA</footnote>

<footnote><sup>c</sup> Institute of Environmental Change and Society, University of Regina, Regina, Saskatchewan S4S 0A2, Canada</footnote>

<footnote><sup>d</sup> Environmental Monitoring and Science Division, Alberta Environment and Parks, Edmonton, Alberta T5J 5C6, Canada</footnote>

<footnote><sup>e</sup> Department of Geography, University of Victoria, Victoria, BC V8W 3R4, Canada</footnote>

<footnote><sup>f</sup> Prairie Environmental Processes Laboratory, Department of Geography and Environmental Studies, University of Regina, Regina, Saskatchewan S4S 0A2, Canada</footnote>

<footnote><sup>g</sup> Institute for Global Food Security, Queen's University Belfast, Belfast, Co. Antrim BT9 5DL, United Kingdom</footnote>

<footnote>* Corresponding author.
E-mail address: haig.a.heather@gmail.com (H.A. Haig).</footnote>
Halder et al., 2013; Pham et al., 2009; Rosa et al., 2016). For example, seasonal variation in river-water isotopes has been used to quantify changes in source water through a single season and among years (Jasechko et al., 2017; Reckerth et al., 2017; Welch et al., 2018). These studies suggest that water isotopes are capable of detecting hydrological variation associated with source waters on the scale of weeks-to-months (Bam and Ireson, 2019; Gibson and Reid, 2014). The capability of water isotopes to integrate the relative contribution of different inflow sources has been highlighted as their advantage over conventional monitoring approaches, particularly in situations where key fluxes are not measured. Similarly, isotopic studies have been successfully used in small wetlands, to determine the relative importance of loss pathways such as infiltration and evaporation (Bam and Ireson, 2019). However, despite these successes, there have been few studies that use water isotopes to quantify short-term variability in hydrology, water source, or water balance (E:I ratio) of lake ecosystems (Cui et al., 2018).

Seasonal records of δ¹⁸O and δ²H with high temporal resolution are rare for lake ecosystems (Cui et al., 2018). For example, Gibson and Reid (2014) quantified sub-monthly changes in δ²H and δ¹⁸O isotope values during a 20-year period in five basins to demonstrate that lakes in the Canadian tundra exhibit up to 4‰ variation in δ¹⁸O values within a single ice-free season. In this case, isotopic shifts were attributed to changes in evaporative forcing rather than variation in inputs within this dry environment (Gibson and Reid, 2014). In contrast, Tyler et al. (2007) analyzed bi-weekly samples over a five-year period in Lake Lochnagar, Scotland, and suggested that the observed 2‰ range in lake water δ¹⁸O values mainly reflected changes in quantity and source of water entering the lake. Similarly, Jones et al. (2016) found that seasonal variability in lake water isotope values can be explained only when models include both evaporation and precipitation sources. Finally, Brooks et al. (2014) surveyed ca. 100 lakes in the continental USA twice during the ice-free period and determined that seasonal variation was insignificant relative to among-lake variation in isotope values. Taken together, these findings suggest that isotopic variability can exceed 5‰ within a year in some regions, due to a combination of evaporation (Bam and Ireson, 2019; Steinman et al., 2010), altered source waters (Pham et al., 2009), and changes in basin connectivity (Brooks et al., 2014). These findings demonstrate the potential of water isotopes to record variation in lake hydrology on a number of timescales (weeks to years), but also illustrate the need to quantify how water isotope values vary across diverse temporal and spatial scales within large lake districts.

The primary goal of this paper was to quantify how temporal variability in water isotopes affects their use as metrics of the water budgets of lakes in a cold sub-humid prairie region. Specifically, we used a 14-year record of bi-weekly (May-August) δ¹⁸O measurements from seven inter-connected lakes to quantify how variation in source waters (rain vs. snow), basin properties (residence time, volume), and landscape position affected estimates of water balance (E:I) within each lake. Comparison among basins that differed by over 100-fold in most hydrological features (Vogt et al., 2011) helped quantify whether regional lakes exhibited common or unique patterns of temporal variation. Finally, this study fills a gap in high-resolution lake water isotope research by providing several records for lakes of the Northern Great Plains. Given the importance of snowmelt to annual hydrological budgets (Pomeroy et al., 2007; Pham et al., 2009) and the pronounced precipitation deficit (>30 cm yr⁻¹) in the study region (Bonsal et al., 2017), we hypothesized that these riverine study lakes would exhibit high seasonal and interannual variation in δ¹⁸O values that greatly affect isotope-derived estimates of lake hydrology (c.f., Gibson and Reid, 2014).

## 2. Methods

### 2.1. Site description

Study lakes are situated in a cold, sub-humid, grassland region of

southern Saskatchewan (SK), Canada (50°00′-51°30′ N, 101°30′-107°10′ W). Basins span over two orders of magnitude in most morphological and hydrological parameters (Table 1). Five lakes (Diefenbaker, Buffalo Pound, Pasqua, Katepwa, Crooked) are part of a central chain of basins following the main stem of the Qu'Appelle River, whereas Last Mountain and Wascana lakes drain into the river via mid-reach tributaries (Fig. 1). Under natural conditions, source waters for the Qu'Appelle River originate from wetlands west of Buffalo Pound Lake and flow eastward to Crooked Lake ~ 400 km downstream (Fig. 1). Since the mid-1960s, water has been withdrawn from the South Saskatchewan River at the Lake Diefenbaker reservoir and used to supplement flow in the Qu'Appelle River. In addition, the hydrology of all Qu'Appelle Valley lakes is managed by the Saskatchewan Water Security Agency (SWSA), with control structures that maintain water level while ensuring adequate water availability downstream (Saskatchewan Water Security Agency, 2013, 2008). Control structures include a bifurcation system at the south end of sub-saline Last Mountain Lake (1–3 g total dissolved solids TDS L⁻¹) with flow direction determined by lake and river levels, as well as water needs within the entire drainage basin. In general, ~80% of annual runoff is associated with spring snowmelt, even though ~ 75% of annual precipitation is received as rain during summer and fall (Coles et al., 2017; Coles and McDonnell, 2018; Fang et al., 2007).

The gross Qu'Appelle River drainage basin encompasses ~ 52,000 km² (excluding Lake Diefenbaker) of agricultural cropland (75% of land cover), with the remaining area including grasslands (12%), surface waters (5%) and the urban centers of Moose Jaw and Regina (<5%) (Vogt et al., 2011). Regional climate is characterized as cool-summer humid continental (Köppen Dfb classification), with annual average temperatures of ~ 3 °C, median relative humidity of ~74%, short warm summers ($T_{july}$ ≈ 8 °C), and cold winters ($T_{jan}$ ≈ -14 °C). In general, mean summer cloud cover (% sky cover) increases, while precipitation deficit (rain minus potential evaporation; Thornthwaite, 1948) declines along a gradient from western headwaters to eastern lakes (Pham et al., 2009). Average summer wind speed at each lake varies from 9 to 15 km h⁻¹ (Vogt et al., 2011), resulting in polymictic conditions in most years, with the exception of Katepwa Lake, which often exhibits thermal stratification in summer.

Hydrological monitoring encompassed a 14-year period (2003–2016) and included daily lake-level records for six of the seven lakes from the Water Survey of Canada (WSC; http://wateroffice.ec.gc.ca/). As lake-level records were not available for Pasqua Lake, we substituted those from Echo Lake, a smaller interconnected basin located <500 m downstream of the Pasqua Lake outlet. Area-capacity curves for all lakes were obtained from SWSA to estimate changes in surface area and volume through time. Inflow for each lake was estimated using gauge and model data from SWSA (Supplementary Table 1), including predictions from the Water Resources Management Model (WRMM) for most years, and measured flows from WSC (http://wateroffice.ec.gc.ca/).

Meteorological conditions were variable during the study period (Supplementary Figs. 1 and 2) with significant seasonal spring flooding (Ahmari et al., 2016; Blais et al., 2016; Stadnyk et al., 2016) and periods of summer aridity (Bonsal and Shabbar, 2008; Wheater and Gober, 2013). Low summer rain resulted in a pronounced precipitation deficit (rain - potential pan evaporation) of over 63 cm in 2008–2009 (Supplementary Fig. 2b). In contrast, heavy snow during winter 2010–2011 and a rapid spring melt resulted in a nearly 1-in-140-year flood in April 2011 (Blais et al., 2016), with damages exceeding $800 million (Brimelow et al., 2014; Wheater and Gober, 2013). Since that time, the Qu'Appelle River drainage basin has experienced elevated summer precipitation (Supplementary Figs. 1 and 2a), thereby reducing the mean precipitation deficit to 32.8 cm (2011–2016) from the long-term average of 46.1 cm (2000–2014) (Supplementary Fig. 2b). Despite high variation among seasons, mean temperatures and relative humidity
<caption>**Table 1**</caption>

Lake characteristics and stable isotope descriptive statistics for the seven study sites. Volume and surface area estimates were made using area capacity curves from Saskatchewan Water Security Agency and lakes levels from the Water Survey of Canada at each site. Drainage basin area has been defined as both GDA (Gross Drainage Area), and SFDA (Sink Free Drainage Area). $\Delta\delta^{18}\text{O}(\text{‰})$, $\Delta\delta^{2}\text{H}(\text{‰})$, $\Delta\text{E}/\text{I}^{18}\text{O}(\text{‰})$ maximum values represent the maximum seasonal and inter-annual disparities relative to mean values as suggested by Cui et al. (2018). All other values are presented as mean ± SD, or by minimum and maximum, separated by commas.

<table><thead><tr><th></th><th>Diefenbaker</th><th>Buffalo Pound</th><th>Last Mountain</th><th>Wascana</th><th>Pasqua</th><th>Katepwa</th><th>Crooked</th></tr></thead><tbody><tr><td>Latitude, Longitude</td><td>51.0, -106.5</td><td>50.6, -105.4</td><td>51.0, -105.2</td><td>50.4, -104.6</td><td>50.8, -104.0</td><td>50.7, -103.6</td><td>50.6, -102.7</td></tr><tr><td>Surface Area (m<sup>2</sup> × 10<sup>6</sup>)</td><td>371 ± 14</td><td>30 ± 0.2</td><td>186 ± 9</td><td>2 ± 0.2</td><td>19 ± 1</td><td>16 ± 0.2</td><td>14 ± 2</td></tr><tr><td>Volume (m<sup>3</sup> × 10<sup>6</sup>)</td><td>7487 ± 350</td><td>93 ± 4</td><td>1863 ± 83</td><td>3 ± 0.5</td><td>117 ± 8</td><td>232 ± 3</td><td>114 ± 9</td></tr><tr><td>SFDA (m<sup>2</sup> × 10<sup>6</sup>)</td><td>8.2 × 10<sup>4</sup></td><td>1.7 × 10<sup>2</sup></td><td>1.2 × 10<sup>3</sup></td><td>27</td><td>1.6 × 10<sup>2</sup></td><td>2.5 × 10<sup>2</sup></td><td>1.6 × 10<sup>2</sup></td></tr><tr><td>GDA (m<sup>2</sup> × 10<sup>9</sup>)</td><td>1.5 × 10<sup>2</sup></td><td>32</td><td>15</td><td>2</td><td>37</td><td>39</td><td>44</td></tr><tr><td>Annual Inflow (m<sup>3</sup> × 10<sup>6</sup>)</td><td>7402 ± 3202</td><td>136 ± 38</td><td>229 ± 157</td><td>34 ± 45</td><td>401 ± 325</td><td>422 ± 363</td><td>476 ± 405</td></tr><tr><td>τ (years)</td><td>1.0 ± 0.3</td><td>0.6 ± 0.1</td><td>6.8 ± 2.6</td><td>0.2 ± 0.4</td><td>0.5 ± 0.3</td><td>0.9 ± 0.5</td><td>0.4 ± 0.3</td></tr><tr><td>d-e × cess (‰)</td><td>1.1 ± 1.4</td><td>-5.0 ± 2.1</td><td>-12.6 ± 2.8</td><td>-5.4 ± 4.6</td><td>-9.0 ± 3.4</td><td>-9.2 ± 2.9</td><td>-8.6 ± 3.4</td></tr><tr><td>Mean δ<sup>18</sup>O (‰)</td><td>-16.9</td><td>-14.2</td><td>-9.6</td><td>-12.4</td><td>-11.6</td><td>-11.6</td><td>-11.4</td></tr><tr><td>Range δ<sup>18</sup>O (‰)</td><td>-18.3, -15.8</td><td>-17.8, -11.9</td><td>-11.4, -8.1</td><td>-19.9, -6.8</td><td>-17.7, -9.1</td><td>-15.4, -10.2</td><td>-13.9, -5.3</td></tr><tr><td>Δδ<sup>18</sup>O (‰) Ma ×, Seasonal</td><td>1.1</td><td>1.6</td><td>0.84</td><td>8.0</td><td>4.6</td><td>3.19</td><td>3.7</td></tr><tr><td>Δδ<sup>18</sup>O (‰) Ma ×, Interannual</td><td>0.94</td><td>1.6</td><td>1.5</td><td>2.4</td><td>2</td><td>2.7</td><td>2.4</td></tr><tr><td>Mean δ<sup>2</sup>H (‰)</td><td>-134</td><td>-119</td><td>-89</td><td>-105</td><td>-102</td><td>-102</td><td>-100</td></tr><tr><td>Range δ<sup>2</sup>H (‰)</td><td>-142, -128</td><td>-141, -104</td><td>-100, -82</td><td>-159, -68</td><td>-141, -89</td><td>-125, -93</td><td>-115, -43</td></tr><tr><td>Δδ<sup>2</sup>H (‰) Ma ×, Seasonal</td><td>7.1</td><td>13.0</td><td>3.1</td><td>56.9</td><td>32.9</td><td>20.3</td><td>38.5</td></tr><tr><td>Δδ<sup>2</sup>H (‰) Ma ×, Interannual</td><td>5.4</td><td>10.3</td><td>8.8</td><td>-17</td><td>13.1</td><td>15.8</td><td>18.3</td></tr><tr><td>Mean E/I <sup>18</sup>O</td><td>0.03</td><td>0.1</td><td>0.24</td><td>0.10</td><td>0.07</td><td>0.07</td><td>0.06</td></tr><tr><td>ΔE:I<sup>18</sup>O Ma ×, Seasonal</td><td>0.03</td><td>0.04</td><td>0.05</td><td>0.13</td><td>0.12</td><td>0.05</td><td>0.09</td></tr><tr><td>ΔE:I<sup>18</sup>O Ma ×, Interannual</td><td>0.02</td><td>0.04</td><td>0.12</td><td>0.07</td><td>0.10</td><td>0.07</td><td>0.07</td></tr></tbody></table>

have remained relatively constant since 2003 (Supplementary Fig. 1). Meteorological variability across the Qu'Appelle River drainage basin generally followed a longitudinal gradient during the study period, with downstream sites in the east receiving greater annual precipitation (mean 420.9 mm) than those in the west (mean 367.4 mm).

Lake water isotope values were compared to those from the local meteoric water line (LMWL) derived from samples collected in Saskatoon, SK (1990–2013) (IAEA/WMO, 2019) and calculated using the precipitation-amount-weighted least squares regression technique (Hughes and Crawford, 2012). In addition, Saskatoon precipitation data (Supplementary Fig. 3) were compared to three years of precipitation samples collected directly at Wascana Lake by the University of Regina (2013–2016). No significant difference was found between isotopic values of local precipitation at Saskatoon and Regina, therefore the more robust decadal data from Saskatoon were used for all subsequent calculations. Briefly, $\delta^{18}\text{O}$ of local precipitation increases from mean values in April of -17.1‰ to -12.6‰ in July before decreasing in August to -14.3‰ (Supplementary Fig. 3). Stable isotope values of winter precipitation (November–April) are more depleted than those of summer rain, with a mean $\delta^{18}\text{O}$ value of -21.1‰ reflecting seasonal changes in atmospheric temperature and the air-mass from which source water is derived. Summer precipitation mainly originates from southern air masses derived from the Gulf of Mexico or Gulf of California, whereas winter precipitation reflects zonal atmospheric flow from the northern Pacific Ocean (Liu and Stewart, 2003).

## 2.2. Isotope analysis

Depth-integrated water samples were collected bi-weekly from seven study lakes between 01 May and 31 August from 2003 to 2016. Samples were filtered sequentially through a GF/C (nominal pore size 1.2 µm) and membrane filter (0.45-µm pore) to remove particulate matter, and stored in airtight amber glass jars at 4 °C to prevent sample evaporation. Samples were analyzed for $\delta^{2}\text{H}$ and $\delta^{18}\text{O}$ following Haig et al. (2020), using a Picarro L2120-I cavity ring-down spectrometer (CRDS), at the Institute of Environmental Change and Society (IECS), University of Regina, Regina, SK, Canada. Isotopic values were standardized to internal and international standards of Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation 2 (SLAP2). All isotopic values were post-processed in a Microsoft Access relational-database application, Laboratory Information Management System (LIMS) for lasers (Coplen and Wassenaar, 2015) and reported in

$\delta$ notation in per mil units (‰) with analytical uncertainty of 0.1‰ for $\delta^{18}\text{O}$ and 0.5‰ for $\delta^{2}\text{H}$.

## 2.3. Isotope mass balances

Comparisons between lentic water isotopes and the LMWL were used to determine the relative importance of water inflow and evaporation to each lake system following Haig et al. (2020). Relationships among parameters in isotope models below were visualized using a flow-chart procedure (Supplementary Fig. 4). Evaporation was also estimated from d-excess values, where d-excess = $\delta^{2}\text{H} - 8\delta^{18}\text{O}$, and values < 10 indicate samples which have undergone evaporation (Dansgaard, 1964). A theoretical local evaporative line (LEL) was established to visualize the projected evolution of lake water isotopes in this region, moving from source waters towards a state of evaporative enrichment. Here the LEL was estimated for each lake individually using a linear regression of isotopic values of weighted-average precipitation (weighted by precipitation amount, $\delta_{P}$), the isotopic composition of lake waters in a closed-basin in which evaporation equals inflow ($\delta_{SS}$), and the theoretical maximum (limiting) isotopic enrichment ($\delta^{*}$). $\delta_{P}$ was calculated using the data from Saskatoon, while $\delta_{SS}$ and $\delta^{*}$ were calculated for each study lake and for each data point using the following formulae:

$$ \delta_{SS}(\delta_{18}\text{O} or \delta_{2}\text{H}) = [\delta_{P} + (m \times \delta^{*})] \times (1 + m)^{-1} \quad (1) $$

$$ \delta_{(\delta_{18}\text{O} or \delta_{2}\text{H})}^{*} = [h\delta_{A} + \varepsilon_{k} + \varepsilon^{*}/\alpha^{*}] \times [h - 10^{-3} \times (\varepsilon_{k} + \varepsilon^{*}/\alpha^{*})]^{-1} \quad (2) $$

$$ m_{(\delta_{18}\text{O} or \delta_{2}\text{H})} = [h - 10^{-3} \times (\varepsilon_{k} + \varepsilon^{*}/\alpha^{*})] \times [1 - h + 10^{-3} \times \varepsilon_{k}]^{-1} \quad (3) $$

(Allison and Leaney, 1982). Here $m$ represents the enrichment slope, $h$ is the relative humidity in decimal notation, $\delta_{A}$ is the isotopic composition of atmospheric moisture (assumed to be in equilibrium with local precipitation), $\varepsilon_{k}$ and $\varepsilon^{*}$ are kinetic and equilibrium isotopic separations, respectively, and $\alpha^{*}$ is the liquid-vapour equilibrium isotopic fractionation factor, which is related to $\varepsilon^{*}$ as $\alpha^{*} = 1 + \varepsilon^{*}/1000$ (Gonfiantini, 1986).

In all equations, both temperature and relative humidity values were flux-weighted, using estimates of potential evapotranspiration (Thornthwaite, 1948) to account for the seasonality of evaporation (Gibson, 2002; Gibson et al., 2016a). Meteorological data for all calculations were obtained from the nearest monitoring station operated by Environment and Climate Change Canada (https://climate.weather.gc).
<figure><img src="image_3.png" /><figcaption>Fig. 1. Map showing location of the Qu'Appelle watershed and study sites, Saskatchewan, Canada. Dominant flow direction is noted by arrows along streams. A bifurcation control structure is located at the south end of Last Mountain Lake with flow direction determined by lake and river levels as well as water needs within the system. Seven study lakes are noted, and gauging stations are indicated as triangles. Tan shading in the inset indicates the Qu'Appelle Valley gross drainage area within central Canada and the USA.</figcaption></figure>

ca). This weighting allows for more a realistic estimate of LEL slope in regions where lakes are ice covered for part of the year. Due to the inherently high temporal variability of climate in our study region, meteorological variables were averaged over an antecedent period equal to the estimated residence time of each lake (Table 1). The use of such a lake-specific time frame uniquely allowed us to relate current lake isotopic values to the appropriate antecedent hydrological and meteorological conditions.

All lakes are assumed to be in instantaneous isotopic and hydrologic steady state, allowing the water balance to be presented as:

$$I = Q + E \qquad (4)$$

$$I\delta_I = (Q \times \delta_Q) + (E \times \delta_E), \qquad (5)$$

where I, Q, and E are lake inflow, outflow, and evaporation (m³) with associated delta values (δI, δQ, and δE) in ‰, respectively summed over the hydrological year. The assumption of steady state is common in isotopic studies despite knowledge that changes in isotopic values or lake volumes are expected during the ice-free season (Gao et al., 2018;

Turner et al., 2014). Recommendations of Yi et al. (2008), and Cui et al. (2018) suggest that samples taken late in the ice-free period minimize deviations from steady state and are most reliable when the objective of the investigation is to obtain estimates of the long-term hydrological status of a system. However, here we sought to use all sample dates to allow quantification of the evolution of isotopic values through the ice-free season and to more objectively assess the validity of steady-state assumptions. In steady state, isotopic values of outflow are the same as that of the lake water (δQ ≈ δL). Therefore, to quantify the water balance for each sample, Eq. (5) was rearranged and Q = 1-E (from Eq. (4)) incorporated to obtain an estimate of the ratio of evaporation to inflow (x):

$$x = \frac{E}{I} = \frac{(\delta_I - \delta_L)}{(\delta_E - \delta_L)}. \qquad (6)$$

Two methods were used to calculate the isotopic composition of inflow, δI, each with a different assumption specific to lake position within the hydrologic landscape (Haig, 2019; Haig et al., 2020). First, for headwater lakes (Diefenbaker, Last Mountain, and Wascana), we
used the coupled isotope tracer (CITM) method of Yi et al., (2008) to calculate δI. This method was also applied to sites (Buffalo Pound) where information on input waters was inadequate to quantify upstream sources. Briefly, this approach uses the intercept between the LMWL and a regression line from sample-specific observations of δE and δL to estimate δI and is conservative with respect to both mass and isotopes (Haig et al., 2020). This method was also used for all lake samples collected after 2014 due to absence of inflow data from the WRMM model for those two years.

Second, for lakes with inflow sources that have previously undergone evaporation (all other basins), we assumed that inflow was not solely of meteoric origin, but rather was also influenced by surface water flows (Haig et al., 2020). Here we sought to directly estimate the isotopic values of input water based on upstream sources to avoid overestimation of the importance of evaporation in the calculation of mass balances. In these basins, we calculated lake δI by assuming no change in lake volume over the sampling period and that this simplified mass balance could be used to estimate inflow (I) from precipitation (P), runoff (R), and upstream flow (J), all in m³, as;

$$I = P + R + J. \qquad (7)$$

Expanding Eq. (7) to consider the isotopic signatures of these sources, yielded the following model:

$$\delta_I = \frac{(P \times \delta_P) + (R \times \delta_R) + (J \times \delta_J)}{P + R + J} \qquad (8)$$

where δI, δP, δR, and δJ are the isotopic composition (‰) of inflow, precipitation, runoff, and upstream source waters, respectively. To estimate runoff (m³), we assumed δR ≈ δP (Gibson and Reid, 2014), an estimation which is appropriate given the location of this study, but may not be appropriate for other lake regions. Water flux was calculated by rearranging equation (7) and using the water balance from Eq. (6);

$$R = \frac{E}{x} - P - J. \qquad (9)$$

Substituting Eqs. (9) and (6) into Eq. (8) allowed calculation of δI for lakes receiving evaporated inflow as:

$$\delta_I = \frac{-J\delta_J\delta_L + E\delta_P\delta_E - E\delta_L\delta_P + J\delta_L\delta_P}{E\delta_E - E\delta_L - J\delta_J + J\delta_P} \qquad (10)$$

where δJ is specific to a given lake. Lake-specific δJ was estimated from mean stream isotopic values or mean upstream lake values, depending on the study basin. Mean values were used in this calculation because upstream contributions were not regularly sampled.

To complete the calculation for the isotopic inputs in Eq. (9), and for calculation of δI using the coupled isotope tracer method, δE was estimated using theoretical models (Craig et al., 1965; Gibson and Edwards, 2002, Supplementary Fig. 4):

$$\delta_E = \frac{(\delta_L - \varepsilon^*) / \alpha^* - h\delta_A - \varepsilon_k}{1 - h + 10^{-3}\varepsilon_k} \qquad (11)$$

where ε* is the isotopic equilibrium separation, α* is the liquid-vapour isotope fractionation factor, and εK is the isotopic kinetic separation. Estimates of ε* were calculated independently for hydrogen and oxygen using known relationships from literature (Horita and Wesolowski, 1994; Majoube, 1971);

$$\begin{aligned} \varepsilon^*_{(\delta^H)} &= 10^3 \ln \alpha \\ &= 1158.8 + (T^3 \times 10^{-9}) - 1620.1 \times (T^2 \times 10^{-6}) + 794.84 \times (T^2 \times 10^{-3}) - 161.04 + 2.9992 \times (10^9 \times T^{-3}) \end{aligned} \qquad (12)$$

$$\begin{aligned} \varepsilon^*_{(\delta^{18}O)} &= 10^3 \ln \alpha \\ &= -7.685 + [6.7123 \times (10^{-3} \times T^{-1})] - [1.6664 \times (10^6 \times T^{-2})] + 0.3504 \times (10^9 \times T^{-3}) \end{aligned} \qquad (13)$$

with temperature (T) in kelvin. The isotopic composition of atmospheric moisture (δA) was estimated using δP from the open-water season (April–November) in the following equation:

$$\delta_A = \frac{\delta_P - \varepsilon^*}{\alpha^*} \qquad (14)$$

Finally, the kinetic enrichment factor (εK), was calculated as;

$$\varepsilon_K = C_K(1-h)(1-h) \qquad (15)$$

where C_K represents the ratio of molecular diffusivities of the heavy and light isotopes of hydrogen (12.5) and oxygen (14.2) (Gonfiantini, 1986), and humidity (h) is represented in decimal notation (Supplementary Fig. 1). Values estimated using Eqs. (12)–(15) were then used to determine the isotopic values of evaporated water (δE) using Eq. (11) (Supplementary Fig. 4).

To fully evaluate the influence of isotopic seasonality on estimates of lake hydrology, a sample-specific antecedent period was determined based on the residence time of each lake. This sample-specific period was used to estimate the meteorological parameters, as well as calculate inflow volumes for each of the equations listed above. By treating each sampling excursion as a unique event, the effect of seasonality on inflow sources (δI) and water mass balance (E/I) could be systematically evaluated.

## 2.4. Statistical analysis

This study used a bivariate copula additive model (Marra and Radice, 2017) to simultaneously model the trends within, and correlation between, δI¹⁸O and δI¹⁸O time series. This approach is based on generalized additive models (GAMs; Wood, 2011; Wood et al., 2016) that are an effective approach for estimation of non-linear seasonal or interannual trends in environmental time series (Augustin et al., 2013; Ferguson et al., 2008; Orr et al., 2015; Pedersen et al., 2019). In GAMs, penalized splines are used to derive the shapes of intra- and interannual trends from the observed time series using a single linear predictor comprising one or more smooth functions, fO, to estimate the change in the mean of the response over time. More recently, GAMs have been extended to allow the estimation of additional parameters of the conditional distribution of the response, such as the variance; these models are known as GAMs for location, scale, shape (GAMLSS; Stasinopoulos et al., 2017). However, GAMLSS are univariate and, as such, do not allow estimation of how the correlation or dependence between two response variables may change through time. Copula-based models (Joe, 2015; Nelsen, 2006) enable the estimation of dependence structures between two or more response variables (Genest and Favre, 2007). Recent developments in statistical theory have combined the properties of GAMLSS and copula-based models for the simultaneous estimation of marginal and joint distributions using smooth functions, resulting in the bivariate copula additive model used herein (Marra and Radice, 2017). In this paper, we restricted our choice of copula function to the Gaussian copula (θ), as it can be interpreted as the correlation between two time series. Details of the application of bivariate copula additive models to δI¹⁸O and δI¹⁸O time series are provided in the Supplementary Information.

# 3. Results

## 3.1. Lake-specific variability in water isotopes

Lake water isotope (δL) values varied substantially among years and study lakes, with the lowest variability recorded in the largest lakes
(Figs. 2 and 3, Table 1). For example, Lake Diefenbaker ($\delta^{18}$O range 2.5‰) and peripheral sub-saline Last Mountain Lake ($\delta^{18}$O range 3.3‰) both exhibited a narrow range in $\delta_L$ values. Wider $\delta^{18}$O seasonal ranges were observed for Buffalo Pound (5.9‰), Pasqua (8.6‰), Katepwa (5.3‰) and Crooked Lakes (8.6‰), with $\delta^{18}$O values generally becoming more enriched as a function of distance downstream (Fig. 2, Table 1). The most dynamic range of lakewater $\delta^{18}$O was observed in rapidly-flushing Wascana Lake, where values ranged from -6.8‰ to -19.9‰ over the ice-free period (Fig. 2). At most sites, the variation in lakewater isotope values lay parallel to the LMWL, except for that in Last Mountain Lake where isotopic variation was positioned between the LMWL and the LEL, suggesting a system more influenced by evaporation (Fig. 2). The longitudinal increase in mean $\delta_L$ values from west to east in lakes joined by the Qu'Appelle River (Diefenbaker, Buffalo Pound, Pasqua, Katepwa, Crooked) suggested a greater importance of rain to downstream lakes and, possibly, a progressive effect of evaporation in sequential basins (Table 1, Supplementary Fig. 5).

## 3.2. Landscape patterns of seasonal isotope variability

The bivariate copula additive model quantified the effects of season (as DoY), year, and lake identity on variability in $\delta_L^{18}$O and $\delta_I^{18}$O (Fig. 3, Supplementary Figs. 5 and 6). Some sites exhibited approximately linear increase in $\delta_L^{18}$O values during the ice-free period (e.g., Buffalo Pound), while isotopic values in other lakes increased to a plateau in late-summer (e.g., Wascana). The extent of seasonal variation in both $\delta_L$ and $\delta_I$ was generally lower in lakes with the largest water volumes (Diefenbaker, Last Mountain, Katepwa) and greater in shallower lakes

with lower volumes (Wascana, Buffalo Pound) (Table 1). In lakes with residence times routinely < 1 year (all except Last Mountain and Diefenbaker), seasonal variability in $\delta_L$ was greater than the among-year variation in mean annual values (Fig. 3, Supplementary Fig. 5). In contrast, lakes with residence times routinely > 1 year were more likely to exhibit greater interannual than seasonal variability in lakewater isotope values. Increases in $\delta_L$ and $\delta_I$ values were observed in most lakes as summer progressed; however, lakewater isotope values also declined mid-summer (ca. DoY 175) in Lake Diefenbaker (Fig. 3), likely reflecting the mid-summer influx of meltwaters from the Rocky Mountains over 500 km to the west (Supplementary Fig. 7). Similarly, late-summer changes in $\delta_I$ values were recorded occasionally in several other lakes (Fig. 3), suggesting a shift in source waters (e.g., changing influence of rain).

Over all lakes and years, the relationship between seasonal variation in lake $\delta^{18}$O values ($\delta_L$; red line) and that of inflow ($\delta_I$; blue line) was often strong (Fig. 3), as indicated by the highly positive $\theta$ values in many lake-year pairs (Fig. 4). Specifically, relationships between $\delta_L$ and $\delta_I$ were strongest in Wascana Lake (mean $\mu_0$ = 0.67), a basin with a short water residence time (annual mean < 3 months, Table 1) and weakest in Last Mountain Lake ($\mu_0$ = -0.04, residence time ~ 7 year). The temporal pattern in $\theta$ values for Last Mountain Lake declined steadily from strong positive to strong negative relationship, leading to a mean value that was close to zero despite years of very strong correspondence between $\delta_L$ and $\delta_I$. In contrast, for lakes with sub-decadal residence times, agreement between $\delta_L$ and $\delta_I$ was greatest during years of high river discharge ($\mu_{0.2010}$ = 0.71, $\mu_{0.2011}$ = 0.68), whereas $\delta_L$ appeared to be more strongly influenced by evaporative fractionation in years with slower water

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Raw water isotope data for Qu'Appelle lakes, including a) Diefenbaker, b) Buffalo Pound, c) Last Mountain, d) Wascana, e) Pasqua, f) Katepwa, and g) Crooked lakes, as well as h) mean values (coloured by lake) and i) precipitation samples from Saskatoon, Canada (rain red, snow blue). Local meteoric water line (LMWL, solid) and lake-specific local evaporative line (LEL, dashed) intercept at the amount-weighted mean isotope value for precipitation ($\delta_p$). For graphical purposes, annual LELs were averaged by lake to present a common LEL during the study period (h, i). Positioning of panels represents the flow path within the catchment, with lakes along the center line representing basins within the central chain and off-axes lakes positioned mid-reach at their point of outflow confluence with the Qu'Appelle River.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Seasonal trends by lake and year for isotopes from lakewater (δL; red) and inflow (δI; blue) as δ¹⁸O ‰. Lines represent the GAM model predictions to the original data (points). Scales are independent in this figure to allow for visualization of lake specific trends over time. See Methods Section 2.2 for calculation details.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Modeled annual theta (θ) values describing the relationship between δL and δI, where 1.0 = high positive correspondence, -1.0 = strong inverse correspondence, and 0.0 = no correspondence. Study lakes colours as Fig. 2.</figcaption></figure>
renewal, leading to divergence of δ<sub>L</sub> and δ<sub>I</sub> (Figs. 3 and 4, Supplementary Fig. 7) and weak relationships (μ<sub>0</sub> 2008 = 0.25) between parameters.

### 3.3. Estimates of lake hydrology

Estimates of water balance (E:I) derived from lakewater isotopes suggested that the hydrologic character of all study lakes was regulated mainly by influx of surface water (Fig. 5a). In general, only sub-saline Last Mountain Lake exhibited elevated E:I ratios consistent with the effects of evaporative concentration. Shallow or small lakes exhibited higher E:I values when calculated from samples taken in late summer than from those collected in spring (Fig. 5b), whereas there were only minor effects of collection season on inferred water balance in large or deep lakes. For example, E:I of Wascana Lake ranged from < 0.05 in spring, to values above 0.2 by fall, whereas both Lake Diefenbaker and Last Mountain Lake showed minimal seasonality in E:I estimates, despite persistent differences in E:I ratios among basins. Differences in E:I values taken in spring and late summer (Fig. 5) were greatest in 2003 (μ = 0.07,

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

**Fig. 5.** Modeled evaporation to inflow ratios based in water isotopes where site is denoted by shape and year by colour. a) Comparison of E:I calculated from May (X axis) and August (Y axis) water samples. Lakes plotting above the solid 1:1 line suggest that analysis of late summer samples suggested a greater contribution of evaporation to water balance than does analysis of spring samples. Observations that plot below the 1:1 line suggesting a higher proportion of late summer inflow in the water balance. Dashed lines represent maximum expected variation due to sampling errors and other non-hydrological factors (Wolfe et al., 2007). b) Violin plots illustrate kernel probability density, i.e. the width of the shape area represents the proportion of the data located at a given E:I value.

σ = 0.12) an arid year associated with elevated water deficit, and least in 2012 (μ = 0.03, σ = 0.02) following a major spring flood event (Supplementary Fig. 2). Overall, variability in E:I in large lakes was less than inter-annual variability estimated for isotopic mass balances in other cold regions (ca. ± 20%; Wolfe et al., 2007; Gibson et al., 2019), while small and shallow systems often exceeded this empirical threshold (Fig. 5a, Table 1).

### 3.4. Seasonal patterns in isotope variability

Analysis of bi-weekly records of lakewater δ<sup>18</sup>O demonstrated that sampling at low temporal resolution did not capture the inter-annual variability in timing, source, and volume of inflow waters (Fig. 3, Supplementary Figs. 6 and 7). In the smallest lakes, isotopic variability among years was greatest early in May and June (Supplementary Fig. 8), likely reflecting large differences in snowpack accumulation, spring melt conditions (e.g., rate of soil thaw), and local spring meteorology (Supplementary Figs. 1 and 2). For example, linear regression analysis of August δ<sup>18</sup>O values with those observed in the following May (Fig. 3) demonstrated only non-significant relationships (r<sup>2</sup> < 0.20, p range > 0.10) for lakes with residence times < 1 year, consistent with a strong effect of winter conditions. In contrast, there was little change from the final sample taken in August to the first sample taken the next May in large or deep systems, with significant (r<sup>2</sup> > 0.6, p < 0.001) predictive relationships between δ<sub>L</sub> August and the following May in Last Mountain, Diefenbaker, and Katepwa lakes. Finally, substantial meteorological events, such as flooding in April 2011, harmonized δ<sup>18</sup>O among lakes, with more depleted values being recorded after inundation in the central lakes linked by the Qu'Appelle River, especially in the eastern half of the Qu'Appelle drainage basin (Fig. 3).

## 4. Discussion

Quantification of water isotopes at bi-weekly resolution demonstrated that variation in the isotopic composition of lake waters occurred because of both changes in source waters and evaporative forcing (Jones et al., 2016; Tyler et al., 2007), but that reliance on a specific form of precipitation (rain or snow), or sensitivity to evaporation, varied widely among lakes within a single river basin. Similarly, all lakes were highly reliant on surface inflow (low E:I) (Figs. 3 and 4, Supplementary Fig. 7), unlike previous surveys of large closed-basin lakes (Pham et al., 2009) and very small wetlands (Bam and Ireson, 2019). As well, small and shallow basins exhibited marked isotopic variation among seasons and years (Fig. 3, Table 1, Supplementary Fig. 5, 6), suggesting that not all fundamental hydrologic properties were constrained by single isotopic measurements. Landscape patterns of lake hydrology were also impacted by extreme events (Fig. 3, Supplementary Fig. 5–7), such as vernal flooding (Blais et al., 2016), a phenomena which will be more common in the future in this region (Asong et al., 2016; Khaliq et al., 2015; Masud et al., 2017). Applications of isotopic approaches demonstrated that the importance of seasonal changes in water source varies with both lake size and landscape position, and suggested that effects of future climate change will be site-specific in lakes of the Northern Great Plains.

### 4.1. Isotopic variability among lakes

Longitudinal variation in δ<sub>L</sub> among lakes (Fig. 2a, b, e, f, g) appeared to arise from a combination of basin-specific variation in the importance of precipitation and, secondarily, evaporative forcing during summer. For example, the high degree of similarity of lakewater δ<sup>18</sup>O values and modelled δ<sub>I</sub> (Fig. 3), low estimates of E:I (median < 0.1; Fig. 5b), and seasonal evolution of lakewater isotope values parallel to the LMWL (Fig. 2), suggest that temporal variability in δ<sub>L</sub> was regulated mainly by site-specific differences in delivery of local precipitation to lakes (Fig. 2i). However, unlike other Qu'Appelle lakes, sub-saline Last
Mountain Lake exhibited the effects of evaporation including generally enriched δL values (Fig. 3), with median E:I ratios two-fold greater than those in other Qu'Appelle lakes (Fig. 5a) and isotopic evolution during summer along the LEL (Fig. 2c). Both these patterns are consistent with those seen in closed-basin lakes in cold prairie regions (e.g., Pham et al., 2009). In previous studies where lake water is subject to significant evaporative forcing, either within the lake itself (Pham et al., 2009) or before inflow enters the basin (Gibson and Reid, 2014, 2010; Wolfe et al., 2007), isotope values usually evolve along the LEL and variation among lakes mainly reflects differences in source water (rain, snow) to individual basins (Yu et al., 2002). Instead, the relatively consistent position of mean annual δL values of individual lakes along the LMWL in this study (Fig. 2h) suggests that Qu'Appelle source waters changed progressively from headwaters to downstream sites, as noted in other lake chains (Yu et al., 2002). In particular, eastern study lakes appeared to exhibit a greater reliance on summer precipitation (Fig. 2i), consistent with the greater cloud cover (Vogt et al., 2011) and precipitation records at those locations over the past 30 years (Akinremi et al., 1999; McCullough et al., 2012). Taken together, these findings suggest that downstream lakes may be particularly sensitive to recently projected changes in the amount, intensity, and seasonality of precipitation during the 21st century (Asong et al., 2016; Bonsal et al., 2017), as well as the more-commonly-expressed concerns about rising temperatures and variation in evaporative forcing (Dibike et al., 2017; Gan and Tanzeeba, 2012; PaiMazumder et al., 2013).

Evaporative enrichment of lakewater isotope ratios may have contributed to longitudinal gradients of δL seen across the east-west axis of the 52,000 km² Qu'Appelle River drainage basin (Fig. 2a, b, e, f, g). For example, the observation that mean d-excess (Cui et al., 2017) decreased from ~ 1‰ to values low as -16.5‰ (Table 1), while individual E:I values ranged up to 0.34, suggested that evaporation also plays a role in shaping the hydrology of downstream sites (Fig. 5, Table 1) despite elevated regional precipitation (McCullough et al., 2012). Given these observations, we suggest that further modelling research may be needed to refine predictions of the relative importance of changes in evaporation and precipitation as controls of soil and surface water availability (Asong et al., 2016; Dibike et al., 2017; Masud et al., 2017; Wang et al., 2018).

## 4.2. Decadal-scale patterns of water source

Analysis of water isotopes since 2003 suggests that the reliance on snowmelt versus rain varies seasonally for many study lakes (Fig. 3, Supplementary Figs. 4 and 5). Prior isotopic research suggests that closed-basin prairie lakes in Canada sampled in mid-summer rely mainly on winter precipitation to persist (Fang et al., 2007; Pham et al., 2009; Pomeroy et al., 2007) as annual precipitation deficits routinely exceed 30 cm yr⁻¹ (Supplementary Fig. 2b). In contrast, the hydrologically-open (E:I < 0.1) and riverine nature of the Qu'Appelle system was characterized by high seasonal variability in δL during spring in the smallest lakes (Fig. 3; Supplementary Fig 8), likely reflecting the high variation in quantity and importance of vernal runoff to total lake volume (Pomeroy et al., 2007). In most years, the transport of winter precipitation into lakes is determined by a complex interaction between precipitation in the antecedent fall and winter, soil moisture in fall and spring, drainage basin size and aspect, and vernal meteorology (temperature, wind) which regulate the rate of snow melt, infiltration, and runoff (Coles et al., 2017; Coles and McDonnell, 2018; Shook and Pomeroy, 2010).

Despite the high variation in spring runoff among years, the observation that δI and δL in headwater lakes (Diefenbaker and Buffalo Pound) consistently plot below the amount-weighted values of precipitation (Fig. 2; δP, δ¹⁸O = -14.97‰, δ²H = -117.3‰), whereas eastern lakes (Katepwa, Crooked) most often plot above δP, demonstrated that there were persistent differences in the relative importance of snow and rain to individual basin hydrology. Such variation was even evident

within a single longitudinal position. For example, δL of Wascana Lake water (Fig. 2d, 3) varied from depleted values characteristic of mid-winter to those representative of only rain (Supplementary Fig. 3), even within a single year (Fig. 3), whereas Last Mountain Lake consistently plotted above local precipitation values, a position that would traditionally be interpreted as lower reliance on snowmelt inflow (Remmer et al., 2018). However, Last Mountain Lake was also isotopically-enriched (Fig. 2) and sub-saline (1–3 g TDS L⁻¹) relative to other Qu'Appelle lakes (Vogt et al., 2011, 2018), resulting in isotopic values evolving along a lake-specific trajectory (Fig. 2c). The intersection of this basin-specific evaporative line and the LMWL at -15.7‰ suggested that Last Mountain Lake receives considerable inputs from local groundwater (δGroundwater = -16.4‰). Together these patterns indicate that local hydrologic features (e.g., effective drainage area, land use, lake morphometry, etc.) interacted with regional gradients of climate (drier in west, wetter in east) to regulate the importance of seasonal precipitation to hydrological budgets of individual lakes.

The importance of winter precipitation to lake hydrology was especially evident during 2011, when the Qu'Appelle River experienced a 1-in-140 year inundation (Blais et al., 2016; Stadnyk et al., 2016). During that flood, the relationship between δL and δI values was strengthened (Fig. 4) and δI converged to a common value, particularly in downstream Pasqua, Katepwa and Crooked lakes (Fig. 3, Supplementary Fig. 6). Subsequently, δI shifted towards values characteristic of local precipitation (δP) (Supplementary Fig. 6), a pattern consistent with meteorological records showing above-average precipitation in both 2010 and 2011 (Supplementary Fig. 2). Such elevated precipitation throughout the annual cycle can result in water-saturated soils in fall, more profound soil freezing in winter, enhanced overland flow in spring, and higher representation of winter precipitation in the overall hydrological budget of lakes in cold regions (Coles et al., 2017; Coles and McDonnell, 2018; Shook et al., 2015).

## 4.3. Seasonal variability in inflow

Variations in δI and δL values among seasons were greater than those among years in small and shallow Qu'Appelle lakes (Table 1; Fig. 3; Supplementary Figs. 5 and 6), but contrasted observations from other lake surveys (Brooks et al., 2014; Gibson et al., 2018; Pham et al., 2009). For example, while large basins such as Diefenbaker and Last Mountain showed a limited deviation from their mean isotopic composition (Δδ¹⁸O < 1.5‰), smaller lakes were much more variable (Wascana Δδ¹⁸O = 8.0‰) (Table 1). Such a dependence of hydrological variability on lake size is not observed in all isotope studies, as large lakes elsewhere (e.g., Lake Chad, Qinghai Lake) exhibit over 3‰ variation around mean isotopic values (Bouchez et al., 2016; Cui et al., 2018; Wu et al., 2015). Thus, while we infer that large or deep lakes may exhibit more constrained isotopic values due to runoff being a smaller proportion of total lake volume (Brooks et al., 2014; Gibson et al., 2016b; Gibson and Edwards, 2002), further research is needed to develop a hierarchical understanding of the importance of factors other than residence time (groundwater, infiltration, evapo-transpiration, etc.) which may influence seasonal variability of water isotopes (Bam and Ireson, 2019).

Observed ranges in δL values (Supplementary Fig. 5) were similar to those recorded from more northern latitudes, where seasonal variation in δL can exceed 4‰ due to evaporative forcing (Gibson and Reid, 2014) or high connectivity of lakes with rivers (Brooks et al., 2014; MacDonald et al., 2017; Remmer et al., 2018; Wolfe et al., 2007). In general, small lakes in such cold regions often shift from snowmelt- to rainfall-dominated inflow during the ice-free season, reflecting the limited storage capacity of small water bodies, low residence time, and high seasonality of overland flow (MacKinnon et al., 2016; Tondu et al., 2013; Turner et al., 2014). Thus while Last Mountain Lake exhibited a sudden change in isotope values associated with the vernal flood of 2011, its δL values have generally varied little within individual years (Figs. 2 and 3; Supplementary Fig. 5, 6) consistent with its 7-year
residence time based on instrumental estimates of surface inflow (Table 1). Instead, the pronounced co-variance of $δ_I$ and $δ_L$ values in small basins (Figs. 3 and 4), combined with their greater seasonal variation, suggests that effects of short-term changes in precipitation will have more profound effects on the small lakes that predominate inland waters (Downing et al., 2006; Hayashi et al., 2016). This effect may be particularly pronounced in lakes of the Canadian Prairies and other regions of low topographic relief, where catchment areas are large relative to lake surface area, and effective drainage areas vary greatly with snowpack depth (Pomeroy et al., 2007; Shook et al., 2015, 2013).

Comparison of variability among months (Supplementary Fig. 8) suggests that samples collected during spring may have limited ability to predict the hydrology in shallow lakes later in the ice-free period (Fig. 5b). In particular, spring values were highly variable for shallow Wascana and Buffalo Pound lakes and showed substantial enrichment during summer (Fig. 3) resulting in underestimations of E:I based on vernal sampling (Fig. 5b). Although this observation is well known from rivers (Cui et al., 2017; Rosa et al., 2016; Smith et al., 2015), many lake surveys are based on single annual samples (Yuan et al., 2011) and rarely consider possible effects of seasonality or changes in reproducibility among years (Gibson and Reid, 2010; Pham et al., 2009; Tyler et al., 2007). The occurrence of high seasonal variability, and the dependence of that variability on lake size, may reduce the effectiveness of comparative studies based on different seasons or between lakes with greatly differing residence times. To minimize bias and to better capture the full range of water balance, we suggest that lakes with unknown residence be sampled for water isotopes at least twice, once immediately after peak spring inflow and again late in summer when inflow is at a minimum.

## 5. Conclusions

This paper explored a 14-year record of bi-weekly water isotope analyses during the ice-free period of seven hydrologically-diverse lakes to quantify how $δ^2H$ and $δ^{18}O$ values vary among seasons, years, and sites in a cold sub-humid catchment in central Canada. Seasonal patterns of water isotopes demonstrated that the importance of individual water sources changed substantially during the ice-free period and as a function of landscape position, with a longitudinal gradient of increasing importance of summer precipitation from west to eastern basins. While the mean water balance could be captured by sub-annual (low resolution) sampling of large basins with water residence times > 1 year, lakes with shorter residence required seasonal sampling to constrain both sources of water and overall water balance. Despite substantial precipitation deficits (Bonsal et al., 2017), we found that water balances in all lakes were dependent mainly on water inflow, during both spring (runoff maximum) and summer (rain). Long-term records also revealed the profound effect of extreme hydrological events (spring flood) which both homogenized regional $δ_I$ and $δ_L$ values and introduced long-lived legacy effects in large basins such as Last Mountain Lake. These findings will be useful to water managers as they illustrate that the importance of future changes in precipitation regimes (Ahmari et al., 2016; Blais et al., 2016; Brimelow et al., 2014) will be predictable on the basis of lake morphometry and geographic position within the hydrological landscape. Further, because E:I values were in the lowest 10th percentile of lakes in the Prairies (Pham et al., 2009) or Boreal Plains (Wolfe et al., 2007; Gibson et al., 2019), we infer that Qu'Appelle basin lakes, which supply water to one third of the provincial population, will be particularly sensitive to future changes in conveyance, precipitation, and runoff (Lapp et al., 2013; Sauchyn et al., 2016).

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

We thank members of the Limnology Laboratory for assistance with data collection since 2003. We also thank Curtis Hallborg and Trent Wurtz of the Saskatchewan Water Security Agency for data on surface flow in the Qu'Appelle River drainage basin, and area capacity-curves. We thank two reviewers for helpful comments on an earlier version of this paper. This work was supported by the NSERC Canada Discovery Grants program, Canada Research Chairs, Canada Foundation for Innovation, the Province of Saskatchewan, and the University of Regina. We acknowledge that the study lakes are located on Treaty 4 and Treaty 6 territory and thank the First Nations of Saskatchewan for sharing and protecting the water resources. This is a contribution of the Qu'Appelle Valley long-term ecological research program (QU-LTER).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.hydroa.2020.100069.

## References

* Ahmari, H., Blais, E.-L., Greshuk, J., 2016. The 2014 flood event in the Assiniboine River Basin: causes, assessment and damage. Can. Water Resour. J. Rev. Can. Ressour. Hydriq. 41 (1-2), 85-93. https://doi.org/10.1080/07011784.2015.1070695.

* Akinremi, O., McGinn, S., Cutforth, H., 1999. Precipitation trends on the Canadian Prairies. J. Clim. 12, 2996-3003. https://doi.org/10.1175/1520-0442(1999)012<2996:PTOTCP>2.0.CO;2.

* Allison, G.B., Leaney, F.W., 1982. Estimation of isotopic exchange parameters, using constant-feed pans. J. Hydrol. 55 (1-4), 151-161. https://doi.org/10.1016/0022-1694(82)90126-3.

* Asong, Z.E., Khalig, M.N., Wheater, H.S., 2016. Projected changes in precipitation and temperature over the Canadian Prairie Provinces using the Generalized Linear Model statistical downscaling approach. J. Hydrol. 539, 429-446. https://doi.org/10.1016/j.jhydrol.2016.05.044.

* Augustin, N.H., Trenkel, V.M., Wood, S.N., Lorance, P., 2013. Space-time modelling of blue ling for fisheries stock management: Space-time modelling of blue ling. Environmetrics 24 (2), 109-119. https://doi.org/10.1002/env.2196.

* Bam, E.K.P., Ireson, A.M., 2019. Quantifying the wetland water balance: A new isotope-based approach that includes precipitation and infiltration. J. Hydrol. 570, 185-200. https://doi.org/10.1016/j.jhydrol.2018.12.032.

* Barnett, T.P., Adam, J.C., Lettenmaier, D.P., 2005. Potential impacts of a warming climate on water availability in snow-dominated regions. Nature 438 (7066), 303-309. https://doi.org/10.1038/nature04141.

* Blais, E.-L., Greshuk, J., Stadnyk, T., 2016. The 2011 flood event in the Assiniboine River Basin: Causes, assessment and damages. Can. Water Resour. J. Rev. Can. Ressour. Hydriq. 41 (1-2), 74-84. https://doi.org/10.1080/07011784.2015.1046139.

* Bonsal, B., Shabbar, A., 2008. Impacts of large-scale circulation variability on low streamflows over canada: a review. Can. Water Resour. J. 33 (2), 137-154. https://doi.org/10.4296/cwrj3302137.

* Bonsal, B.R., Cuell, C., Wheaton, E., Sauchyn, D.J., Barrow, E., 2017. An assessment of historical and projected future hydro-climatic variability and extremes over southern watersheds in the Canadian Prairies: Canadian Prairie hydro-climatic variability and extremes. Int. J. Climatol. 37 (10), 3934-3948. https://doi.org/10.1002/joc.4967.

* Bouchez, C., Goncalves, J., Deschamps, P., Vallet-Coulomb, C., Hamelin, B., Doumnang, J.-C., Sylvestre, F., 2016. Hydrological, chemical, and isotopic budgets of Lake Chad: A quantitative assessment of evaporation, transpiration and infiltration fluxes. Hydrol. Earth Syst. Sci. 20, 1599-1619. https://doi.org/10.5194/hess-20-1599-2016.

* Brimelow, J., Stewart, R., Hanesiak, J., Kochtubajda, B., Szeto, K., Bonsal, B., 2014. Characterization and assessment of the devastating natural hazards across the Canadian Prairie Provinces from 2009 to 2011. Nat. Hazards 73 (2), 761-785. https://doi.org/10.1007/s11069-014-1107-6.

* Brooks, J.R., Gibson, J.J., Birks, S.J., Weber, M.H., Rodecap, K.D., Stoddard, J.L., 2014. Stable isotope estimates of evaporation: Inflow and water residence time for lakes across the United States as a tool for national lake water quality assessments. Limnol. Oceanogr. 59 (6), 2150-2165. https://doi.org/10.4319/lo.2014.59.6.2150.

* Coles, A.E., McConkey, B.G., McDonnell, J.J., 2017. Climate change impacts on hillslope runoff on the northern Great Plains, 1962-2013. J. Hydrol. 550, 538-548. https://doi.org/10.1016/j.jhydrol.2017.05.023.

* Coles, A.E., McDonnell, J.J., 2018. Fill and spill drives runoff connectivity over frozen ground. J. Hydrol. 558, 115-128. https://doi.org/10.1016/j.jhydrol.2018.01.016.

* Coplen, T.B., Wassenaar, L.I., 2015. LIMS for lasers 2015 for achieving long-term accuracy and precision of $^2H$, $^{17}O$, and $^{18}O$ of waters using laser absorption spectrometry. Rapid Comm. Mass. Spec. 29, 2122-2130. https://doi.org/10.1002/ rcm.7372.

* Coulibaly, P., Samuel, J., Pietroniro, A., Harvey, D., 2013. Evaluation of Canadian National Hydrometric Network density based on WMO 2008 standards. Can. Water Resour. J. 38 (2), 159-167. https://doi.org/10.1080/07011784.2013.787181.
Craig, H., Gordon, L., Tongiorgi, E., 1965. Deuterium and oxygen 18 variations in the
ocean and the marine atmosphere; in: Stable Isotopes in Oceanographic Studies and
Paleotemperatures. Laboratorio di Geologia Nucleare, pp. 9-130.

Cui, J., Tian, L., Biggs, T.W., Wen, R., 2017. Deuterium-excess determination of evaporation to inflow ratios of an alpine lake: Implications for water balance and modeling. Hydrol. Process. 31 (5), 1034–1046. https://doi.org/10.1002/hyp.11085.

Cui, J., Tian, L., Gibson, J.J., 2018. When to conduct an isotopic survey for lake water balance evaluation in highly seasonal climates. Hydrol. Process. 32 (3), 379–387. https://doi.org/10.1002/hyp.11420.

Dansgaard, W., 1964. Stable isotopes in precipitation. Tellus 4, 436-467. https://doi.org/10.1111/j.2153-3490.1964.tb00181.x.

Dibike, Y., Prowse, T., Bonsal, B., O'Neil, H., 2017. Implications of future climate on water availability in the western Canadian river basins: Implications of future climate on water availability in western Canada. Int. J. Climatol. 37 (7), 3247–3263. https://doi.org/10.1002/joc.4912.

Downing, J.A., Prairie, Y.T., Cole, J.J., Duarte, C.M., Tranvik, L.J., Striegl, R.G., McDowell, W.H., Kortelainen, P., Caraco, N.F., Melack, J.M., Middelburg, J.J., 2006. The global abundance and size distribution of lakes, ponds, and impoundments. Limnol. Oceanogr. 51 (5), 2388–2397. https://doi.org/10.4319/lo.2006.51.5.2388.

Fang, X., Minke, A., Pomeroy, J., Brown, T., Westbrook, C., Guo, X., Guangul, S., 2007. A review of Canadian prairie hydrology: Principles, modelling and response to land use and drainage change (No. 2), Center for Hydrology. University of Saskatchewan, Saskatoon, Canada.

Ferguson, C.A., Carvalho, L., Scott, E.M., Bowman, A.W., Kirika, A., 2008. Assessing ecological responses to environmental change using statistical models. J. Appl. Ecol. 45, 193–203. https://doi.org/10.1111/j.1365-2664.2007.01428.x.

Gan, T., Tanzeeba, S., 2012. Potential impact of climate change on the water availability of South Saskatchewan River Basin. Clim. Change 112, 355–386. https://doi.org/10.1007/s10584-011-0221-7.

Gao, Z., Niu, F., Lin, Z., Luo, J., Yin, G., Wang, Y., 2018. Evaluation of thermokarst lake water balance in the Qinghai–Tibet Plateau via isotope tracers. Sci. Total Environ. 636, 1–11. https://doi.org/10.1016/j.scitotenv.2018.04.103.

Genest, C., Favre, A.-C., 2007. Everything you always wanted to know about copula modeling but were afraid to ask. J. Hydrol. Eng. 12 (4), 347–368. https://doi.org/10.1061/(ASCE)1084-0699(2007)12:4(347).

Gibson, J.J., 2002. A new conceptual model for predicting isotopic enrichment of lakes in seasonal climates. PAGES News 10 (2), 10–11. https://doi.org/10.22498/pages.10.2.10.

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

Gibson, J.J., Birks, S.J., Yi, Y., Moncur, M.C., McEachern, P.M., 2016b. Stable isotope mass balance of fifty lakes in central Alberta: Assessing the role of water balance parameters in determining trophic status and lake level. J. Hydrol. Reg. Stud. 6, 13–25. https://doi.org/10.1016/j.ejrh.2016.01.034.

Gibson, J.J., Birks, S.J., Yi, Y., Shaw, P., Moncur, M.C., 2018. Isotopic and geochemical surveys of lakes in coastal B.C.: Insights into regional water balance and water quality controls. J. Hydrol. Reg. Stud. 17, 47–63. https://doi.org/10.1016/j.ejrh.2018.04.006.

Gibson, J.J., Edwards, T.W.D., 2002. Regional water balance trends and evaporation-transpiration partitioning from a stable isotope survey of lakes in northern Canada: Regional water balance using stable isotopes. Global Biogeochem. Cycles 16 (2), 10-1-10-14. https://doi.org/10.1029/2001GB001839.

Gibson, J.J., Reid, R., 2010. Stable isotope fingerprint of open-water evaporation losses and effective drainage area fluctuations in a subarctic shield watershed. J. Hydrol. 381 (1-2), 142–150. https://doi.org/10.1016/j.jhydrol.2009.11.036.

Gibson, J.J., Reid, R., 2014. Water balance along a chain of tundra lakes: A 20-year isotopic perspective. J. Hydrol. 519, 2148–2164. https://doi.org/10.1016/j.jhydrol.2014.10.011.

Gibson, J.J., Birks, S.J., Moncur, M., 2019. Mapping water yield distribution across the South Athabasca Oil Sands (SAO) area: Baseline surveys applying isotope mass balance of lakes. J. Hydrol. Reg. Stud. 21, 1–13. https://doi.org/10.1016/j.ejrh.2018.11.001.

Gonfiantini, R., 1986. Environmental isotopes in lake studies, in: Handbook of Environmental Isotope Geochemistry. Elsevier, pp. 113–168.

Haig, H.A., 2019. Analysis of lakewater isotopes in the northern Great Plains: Insights from long-term monitoring and spatial surveys. PhD Thesis University of Regina, Regina, Canada, p. 179 pp.

Haig, H.A., Hayes, N.M., Simpson, G.L., Wissel, B., Hodder, K.R., Leavitt, P.R., 2020. Comparison of isotopic mass balance and instrumental techniques as estimates of basin hydrology in seven connected lakes over 12 years. J. Hydrol. X 6, 100046. https://doi.org/10.1016/j.hydroa.2019.100046.

Halder, J., Decrouy, L., Vennemann, T., 2013. Mixing of Rhone River water in lake geneva (Switzerland-France) inferred from stable hydrogen and oxygen isotope profiles. J. Hydrol. 477, 152–164. https://doi.org/10.1016/j.jhydrol.2012.11.026.

Hayashi, M., van der Kamp, G., Rosenberry, D.O., 2016. Hydrology of prairie wetlands: understanding the integrated surface-water and groundwater processes. Wetlands 36 (S2), 237–254. https://doi.org/10.1007/s13157-016-0797-9.

Hughes, C.E, Crawford, J., 2012. A new precipitation weighted method for determining the meteoric water line for hydrological applications demonstrated using Australian and global GNIP data. J. Hydrol. 464–465, 344–351. https://doi.org/10.1016/j.jhydrol.2012.07.029.

IAEA, WMO, Global Network of Isotopes in Precipitation 2019 The GNIP Database Accessible at: https://nucleus.iaea.org/wiser.

Jasechko, S., Perrone, D., Befus, K.M., Bayani Cardenas, M., Ferguson, G., Gleeson, T., Luijendijk, E., McDonnell, J.J., Taylor, R.G., Wada, Y., Kirchner, J.W., 2017. Global aquifers dominated by fossil groundwaters but wells vulnerable to modern contamination. Nature Geosci. 10, 425–429. https://doi.org/10.1038/ngeo2943.

Joe, H., 2015. Deposition Modeling with Copulas. Monographs on Statistics and Applied Probability 134, Published by CRC Press. 462 pp. ISBN: 978-1-4665-8322-1.

Jones, M.D., Cuthbert, M.O., Leng, M.J., McGowan, S., Mariethoz, G., Arrowsmith, C., Sloane, H.J., Humphrey, K.K., Cross, I., 2016. Comparisons of observed and modelled lake δ¹⁸O variability. Quat. Sci. Rev. 131, 329–340. https://doi.org/10.1016/j.quascirev.2015.09.012.

Khaliq, M.N., Sushama, L., Monette, A., Wheater, H., 2015. Seasonal and extreme precipitation characteristics for the watersheds of the Canadian Prairies Provinces as simulated by the NARCCAP multi-RCM ensemble. Clim. Dyn. 44 (1-2), 255–277. https://doi.org/10.1007/s00382-014-2235-0.

Kirchner, J.W., 2006. Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology: Getting the right answers for the right reasons. Water Resour. Res. 42 (3) https://doi.org/10.1029/2005WR004362.

Lapp, S.L., Jacques, J.-M., Sauchyn, D., Vanstone, J., 2013. Forcing of hydroclimatic variability in the northwestern Great Plains since AD 1406. Quat. Int. 310, 47–61. https://doi.org/10.1016/j.quaint.2012.09.011.

Liu, J., Stewart, R.E., 2003. Water vapor fluxes over the Saskatchewan River basin. J. Hydrometeorol. 4, 944–959. https://doi.org/10.1175/1525-7541(2003)004-0944:WVFOTS-2.0:CO2;

MacDonald, L.A., Wolfe, B., Turner, K.W., Anderson, L., Arp, C.D., Birks, S.J., Bouchard, F., Edwards, T.W.D., Farquharson, N., Hall, R.I., McDonald, I., Narancic, B., Ouimet, C., Pienitz, R., Tondu, J., White, H., 2017. A synthesis of thermokarst lake water balance in high-latitude regions of North America from isotope tracers. Arct. Sci. 3 (2), 118–149. https://doi.org/10.1139/as-2016-0019.

MacKinnon, B.D., Sagin, J., Baulch, H.M., Lindenschmidt, K.-E., Jardine, T.D., Vinebroke, R., 2016. Influence of hydrological connectivity on winter limnology in floodplain lakes of the Saskatchewan River Delta, Saskatchewan. Can. J. Fish. Aquat. Sci. 73 (1), 140–152. https://doi.org/10.1139/cjas-2015-0210.

Majουbe, M., 1971. Fracturation en oxygène 18 et en deuterium entre l'eau et sa vapeur. J. Chimie Phys. 68, 1423–1436.

Marra, G., Radice, R., 2017. Bivariate copula additive models for location, scale and shape. Comput. Stat. Data Anal. 112, 99–113. https://doi.org/10.1016/j.csdta.2017.03.004.

Masub, M.B., Khaliq, M.N., Wheater, H.S., 2017. Projected changes to short- and long-duration precipitation extremes over the Canadian Prairies Provinces. Clim. Dyn. 49 (5-6), 1597–1616. https://doi.org/10.1007/s00382-016-3404-0.

McCullough, G.K., Page, S.J., Hesslein, R.H., Stainton, M.P., Kling, H.J., Salki, A.G., Barber, D.G., 2012. Hydrological forcing of a recent trophic surge in Lake Winnipeg. J. Great Lakes Res. 38, 95–105. https://doi.org/10.1016/j.gibr.2011.12.012.

Nelsen, R.B., 2006. An Introduction to Copulas. Springer Series in Statistics. Nature Springer. https://doi.org/10.1007/0-387-28678-0.

Orr, H.G., Simpson, G.L., Clers, S., Watts, G., Hughes, M., Hannaford, J., Dunbar, M.J., Laizé, C.L.R., Wilby, R.L., Battarbee, R.W., Evans, R., 2015. Detecting changing river temperatures in England and Wales. Hydrol. Process. 29 (5), 752–766. https://doi.org/10.1002/hyp.10181.

PaiMazumder, D., Sushama, L., Laprise, R., Khaliq, M.N., Sauchyn, D., 2013. Canadian RCM projected changes to short- and long-term drought characteristics over the Canadian Prairies. Int. J. Climatol. 33 (6), 1409–1423. https://doi.org/10.1002/joc.3521.

Pedersen, E.J., Miller, D.L., Simpson, G.L., Ross, N. 2019. Hierarchical generalized additive models in ecology: An introduction with mgcv. PeerJ 7, e6876. https://doi.org/10.7717/peerj.6876.

Pham, S.V., Leavitt, P.R., McGowan, S., Wissel, B., Wassenaar, L.I., 2009. Spatial and temporal variability of prairie lake hydrology as revealed by stable isotopes of hydrogen and oxygen. Limnol. Oceanogr. 54 (1), 101–118. https://doi.org/10.4319/lo.2009.54.1.0101.

Pomeroy, J.W., Gray, D.M., Brown, T., Hedstrom, N.R., Quinton, W.L., Granger, R.J., Carey, S.K., 2007. The cold regions hydrological model: A platform for basing process representation and model structure on physical evidence. Hydrol. Process. 21 (19), 2650–2667. https://doi.org/10.1002/hyp.6787.

Reckerth, A., Stichler, W., Schmidt, A., Stumpp, C., 2017. Long-term data set analysis of stable isotopic composition in German rivers. J. Hydrol. 552, 718–731. https://doi.org/10.1016/j.jhydrol.2017.07.022.

Remmer, C.R., Klemt, W.H., Wolfe, B.B., Hall, R.I., 2018. Inconsequential effects of flooding in 2014 on lakes in the Peace-Athabasca Delta (Canada) due to long-term drying. Limnol. Oceanogr. 63 (4), 1502–1518. https://doi.org/10.1002/ino.10787.

Rosa, E., Hillaire-Marcel, C., Hélie, J.-F., Myre, A., 2016. Processes governing the stable isotope composition of water in the St. Lawrence river system, Canada. Isot. Environ. Health Stud. 52 (4-5), 370–379. https://doi.org/10.1080/10256016.2015.1135138.

Saskatchewan Water Security Agency, SWSA, 2008. Getting to the Source: Upper Qu'Appelle River and Wascana Creek watersheds source water protection plan. Saskatchewan Water Security Agency, Regina, Canada, p. 69.

Saskatchewan Water Security Agency, SWSA, 2013. Lower Qu'Appelle River Watershed Plan. Saskatchewan Water Security Agency, Regina, Canadsa, p. 40.

Sauchyn, D.J., St, J., Jacques, E., Barrow, M.W., Nemeth, R.J., MacDonald, M.A., Sheer, D.P. Sheer, 2016. Adaptive water resource planning in the South Saskatchewan River Basin: Use of scenarios of hydroclimatic variability and extremes. J. Am. Water Resour. Assoc. 52, 222–240. https://doi.org/10.1111/1752-1688.12378.

Seibert, J., McDonnell, J.J., 2002. On the dialog between experimentalist and modeller in catchment hydrology: Use of soft data for multi-criteria model calibration. Water Resour. Res. 38, 23-1-23-14. https://doi.org/10.1029/2001wr000978.
Shook, K., Pomeroy, J., 2010. Hydrological effects of the temporal variability of the multi-scaling of snowfall on the Canadian Prairies. Hydrol. Earth Syst. Sci. 14, 1195-1203. https://doi.org/10.5194/hess-14-1195-2010.

Shook, K., Pomeroy, J., van der Kamp, G., 2015. The transformation of frequency distributions of winter precipitation to spring streamflow probabilities in cold regions; case studies from the Canadian Prairies. J. Hydrol. 521, 395-409. https://doi.org/10.1016/j.jhydrol.2014.12.014.

Shook, K., Pomeroy, J.W., Spence, C., Boychuk, L., 2013. Storage dynamics simulations in prairie wetland hydrology models: Evaluation and parameterization. Hydrol. Process. 27, 1875-1889. https://doi.org/10.1002/hyp.9867.

Smith, A., Delavau, C., Stadnyk, T., 2015. Identification of geographical influences and flow regime characteristics using regional water isotope surveys in the lower Nelson River, Canada. Can. Water Resour. J. / Revue canadienne des ressources hydriques 40 (1), 23-35. https://doi.org/10.1080/07011784.2015.1008048.

Stadnyk, Tricia, Dow, Karen, Wazney, Lucas, Blais, Eric-Lorne, 2016. The 2011 flood event in the Red River Basin: Causes, assessment and damages. Can. Water Resour. J. Rev. Canad. Ressour. Hydriq. 41 (1-2), 65-73. https://doi.org/10.1080/07011784.2015.1008048.

Stasinopoulos, M.D., Rigby, R.A., Heller, G.Z., Voudouris, V., Bastiani, F., 2017. Flexible regression and smoothing: Using GAMLSS in R, Chapman and Hall/CRC Press. Chapman and Hall/CRC Press. 549 pp. ISBN 978-1-138-19790-9.

Steinman, B.A., Rosenmeier, M.F., Abbott, M.B., 2010. The isotopic and hydrologic response of small, closed-basin lakes to climate forcing from predictive models: Simulations of stochastic and mean state precipitation variations. Limnol. Oceanogr. 55, 2246-2261. https://doi.org/10.4319/lo.2010.55.6.2246.

Thornhwaite, C.W., 1948. An approach toward a rational classification of climate. Geogr. Rev. 38, 55-94. https://doi.org/10.2307/210739.

Tondu, J.M.E., Turner, K.W., Wolfe, B.B., Hall, R.I., Edwards, T.W.D., McDonald, I., 2013. Using water isotope tracers to develop the hydrological component of a long-term aquatic ecosystem monitoring program for a northern lake-rich landscape. Arct. Antarct. Alp. Res. 45 (4), 594-614. https://doi.org/10.1657/1938-4246-45.4.594.

Tunaley, C., Tetzlaff, D., Birkel, C., Soulsby, C., 2017. Using high-resolution isotope data and alternative calibration strategies for a tracer-aided runoff model in a nested catchment. Hydrol. Process. 31 (22), 3962-3978. https://doi.org/10.1002/hyp.11313.

Turner, K.W., Wolfe, B.B., Edwards, T.W.D., Lantz, T.C., Hall, R.I., Larocque, G., 2014. Controls on water balance of shallow thermokarst lakes and their relations with catchment characteristics: A multi-year, landscape-scale assessment based on water isotope tracers and remote sensing in Old Crow Flats, Yukon (Canada). Glob. Change Biol. 20 (5), 1585-1603. https://doi.org/10.1111/gcb.12465.

Tyler, J.J., Leng, M.J., Arrowsmith, C., 2007. Seasonality and the isotope hydrology of Lochnagar, a Scottish mountain lake: implications for palaeoclimate research. Holocene 17 (6), 717-727. https://doi.org/10.1177/0959683607080513.

Vogt, R.J., Rusak, J.A., Patoine, A., Leavitt, P.R., 2011. Differential effects of energy and mass influx on the landscape synchrony of lake ecosystems. Ecology 92 (5), 1104-1114. https://doi.org/10.1890/10-1846.1.

Vogt, R.J., Sharma, S., Leavitt, P.R., 2018. Direct and interactive effects of climate, meteorology, river hydrology, and lake characteristics on water quality in productive lakes of the Canadian Prairies. Can. J. Fish. Aquat. Sci. 75 (1), 47-59. https://doi.org/10.1139/cjfas-2016-0520.

Wang, J., Song, C., Reager, J.T., Yao, F., Famiglietti, J.S., Sheng, Y., MacDonald, G.M., Brun, F., Schmied, H.M., Marston, R.A., Wada, Y., 2018. Recent global decline in endorheic basin water storages. Nature Geosci. 11 (12), 926-932. https://doi.org/10.1038/s41561-018-0265-7.

Welch, C., Smith, A.A., Stadnyk, T.A., 2018. Linking physiography and evaporation using the isotopic composition of river water in 16 Canadian boreal catchments. Hydrol. Process. 32 (2), 170-184. https://doi.org/10.1002/hyp.11396.

Wheater, H., Gober, P., 2013. Water security in the Canadian Prairies: science and management challenges. Philos. Trans. R. Soc. A 371 (2002), 20120409. https://doi.org/10.1098/rsta.2012.0409.

Wolfe, B.B., Karst-Riddoch, T.L., Hall, R.I., Edwards, T.W.D., English, M.C., Palmini, R., McGowan, S., Leavitt, P.R., Vardy, S.R., 2007. Classification of hydrological regimes of northern floodplain basins (Peace-Athabasca Delta, Canada) from analysis of stable isotopes (δ¹⁸O, δ²H) and water chemistry. Hydrol. Process. 21 (2), 151-168. https://doi.org/10.1002/hyp.6229.

Wood, E.F., Roundy, J.K., Troy, T.J., van Beek, L.P.H., Bierkens, M.F.P., Blyth, E., de Roo, A., Döll, P., Ek, M., Famiglietti, J., Gochis, D., van de Giesen, N., Houser, P., Jaffé, P.R, Kollet, S., Lehner, B., Lettenmaier, D.P., Peters-Lidard, C., Sivapalan, M., Sheffield, J., Wade, A., Whitehead, P., 2011. Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth's terrestrial water: OPINION. Water Resour. Res. 47 (5) https://doi.org/10.1029/2010WR010090.

Wood, S.N., 2011. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. Roy. Stat. Soc. B Stat. Methodol. 73, 3-36. https://doi.org/10.1111/j.1467-9868.2010.00749.x.

Wood, S.N., Pya, N., Säfken, B., 2016. Smoothing parameter and model selection for general smooth models. J. Am. Stat. Assoc. 111 (516), 1548-1563. https://doi.org/10.1080/01621459.2016.1180986.

Wu, W., Li, L., Li, L., Jiang, J., Li, L., Liu, L., 2015. Evaporative enrichment of stable isotopes (¹⁸O and D) in lake water and the relation to lake-level change of Lake Qinghai, Northeast Tibetan Plateau of China. J. Arid. Land. 7, 623-635. https://doi.org/10.1007/s40333-015-0048-6.

Yi, Y., Brock, B.E., Falcone, M.D., Wolfe, B.B., Edwards, T.W.D., 2008. A coupled isotope tracer method to characterize input water to lakes. J. Hydrol. 350 (1-2), 1-13. https://doi.org/10.1016/j.jhydrol.2007.11.008.

Yu, Z., Ito, E., Engstrom, D.R., 2002. Water isotopic and hydrochemical evolution of a lake chain in the northern Great Plains and its paleoclimatic implications. J. Paleolimnol. 28, 207-217. https://doi.org/10.1023/A:102165182.

Yuan, F., Sheng, Y., Yao, T., Fan, C., Li, J., Zhao, H., Lei, Y., 2011. Evaporative enrichment of oxygen-18 and deuterium in lake waters on the Tibetan Plateau. J. Paleolimnol. 46 (2), 291-307. https://doi.org/10.1007/s10933-011-9540-y.