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

Stable isotopes of hydrogen ( δ 2 H) and oxygen ( δ 18 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 2 Qu ’ Appelle River drainage basin of the Canadian Prairies. Analyses revealed that within-year variation in δ 18 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 (cid:0) 1 , 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.


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 (δ 18 O, δ 2 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, δ 18 O and δ 2 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(Gibson et al., , 2016b(Gibson et al., , 2016a;;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 δ 18 O and δ 2 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 δ 2 H and δ 18 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 δ 18 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 δ 18 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 14year record of bi-weekly (May-August) δ 18 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 − 1 ) in the study region (Bonsal et al., 2017), we hypothesized that these riverine study lakes would exhibit high seasonal and interannual variation in δ 18 O values that greatly affect isotope-derived estimates of lake hydrology (c.f., Gibson and Reid, 2014).

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 midreach 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 -1 ) 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 2 (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 − 1 (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)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(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 H.A. Haig et al. 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 (1990SK ( -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 (2013Regina ( -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, δ 18 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 δ 18 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).

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 δ 2 H and δ 18 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 δ notation in per mil units (‰) with analytical uncertainty of 0.1‰ for δ 18 O and 0.5‰ for δ 2 H.

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 = δ 2 H -8 δ 18 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, δ P ), the isotopic composition of lake waters in a closedbasin in which evaporation equals inflow (δ SS ), and the theoretical maximum (limiting) isotopic enrichment (δ*).δ P was calculated using the data from Saskatoon, while δ SS and δ* were calculated for each study lake and for each data point using the following formulae: and Leaney, 1982).Here m represents the enrichment slope, h is the relative humidity in decimal notation, δ A is the isotopic composition of atmospheric moisture (assumed to be in equilibrium with local precipitation), ε k and ε* are kinetic and equilibrium isotopic separations, respectively, and α* is the liquid-vapour equilibrium isotopic fractionation factor, which is related to ε* as α* = 1 + ε*/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.

Table 1
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).Δδ 18 O(‰),Δδ 2 H(‰),ΔE/I 18 O(%) 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.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: where I, Q, and E are lake inflow, outflow, and evaporation (m 3 ) 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 icefree 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): 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 Fig. 1.Map showing location of the Qu'Appelle watershed and study sites, Saskatchewan, Canada.Dominate 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.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 3 , as; Expanding Eq. ( 7) to consider the isotopic signatures of these sources, yielded the following model: where δ I , δ P , δ R, and δ J are the isotopic composition (‰) of inflow, precipitation, runoff, and upstream source waters, respectively.To estimate runoff (m 3 ), 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); Substituting Eqs. ( 9) and (6) into Eq.( 8) allowed calculation of δ I for lakes receiving evaporated inflow as: 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): 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); 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: Finally, the kinetic enrichment factor (ε K ), was calculated as; 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.

Statistical analysis
This study used a bivariate copula additive model (Marra and Radice, 2017) to simultaneously model the trends within, and correlation between, δ L 18 O and δ I 18 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, f(), 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 δ L 18 O and δ I 18 O time series are provided in the Supplementary Information.

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 H.A. Haig et al. (Figs. 2 and 3, Table 1).For example, Lake Diefenbaker (δ 18 O range 2.5‰) and peripheral sub-saline Last Mountain Lake (δ 18 O range 3.3‰) both exhibited a narrow range in δ L values.Wider δ 18 O seasonal ranges were observed for Buffalo Pound (5.9‰), Pasqua (8.6‰), Katepwa (5.3‰) and Crooked Lakes (8.6‰), with δ 18 O values generally becoming more enriched as a function of distance downstream (Fig. 2, Table 1).The most dynamic range of lakewater δ 18 O was observed in rapidlyflushing 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 δ 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).

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 δ L 18 O and δ I 18 O (Fig. 3, Supplementary Figs. 5 and 6).Some sites exhibited approximately linear increase in δ L 18 O values during the ice-free period (e.g., Buffalo Pound), while isotopic values in other lakes increased to a plateau in latesummer (e.g., Wascana).The extent of seasonal variation in both δ L and δ 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 δ 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 δ L and δ 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 δ 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 δ 18 O values (δ L; red line) and that of inflow (δ I ; blue line) was often strong (Fig. 3), as indicated by the highly positive θ values in many lake-year pairs (Fig. 4).Specifically, relationships between δ L and δ I were strongest in Wascana Lake (mean μ θ = 0.67), a basin with a short water residence time (annual mean < 3 months, Table 1) and weakest in Last Mountain Lake (μ θ = -0.04,residence time ~ 7 year).The temporal pattern in θ 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 δ L and δ I .In contrast, for lakes with sub-decadal residence times, agreement between δ L and δ I was greatest during years of high river discharge (μ θ 2010 = 0.71, μ θ 2011 = 0.68), whereas δ L appeared to be more strongly influenced by evaporative fractionation in years with slower water    renewal, leading to divergence of δ L and δ I (Figs. 3 and 4, Supplementary Fig. 7) and weak relationships (μ θ 2008 = 0.25) between parameters.

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, σ = 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).

Seasonal patterns in isotope variability
Analysis of bi-weekly records of lakewater δ 18 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 δ 18 O values with those observed in the following May (Fig. 3) demonstrated only non-significant relationships (r 2 < 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 2 > 0.6, p < 0.001) predictive relationships between δ L August and the following May in Last Mountain, Diefenbaker, and Katepwa lakes.Finally, substantial meteorological events, such as flooding in April 2011, harmonized δ 18 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).

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.

Isotopic variability among lakes
Longitudinal variation in δ L 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 δ 18 O values and modelled δ I (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 δ L 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 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 nonhydrological 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.

H.A. Haig et al.
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 andReid, 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 2 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).

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 − 1 (Supplementary Fig. 2b).In contrast, the hydrologicallyopen (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, δ 18 O = -14.97‰,δ 2 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 midwinter 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 -1 ) relative to other Qu'Appelle lakes (Vogt et al., 2011(Vogt et al., , 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 1in-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).

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 (Δδ 18 O < 1.5‰), smaller lakes were much more variable (Wascana Δδ 18 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 rainfalldominated 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., 2015Shook et al., , 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.

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 δ 2 H 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.

Fig. 2 .
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 (δ 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.
H.A.Haig et al.

Fig. 3 .
Fig.3.Seasonal trends by lake and year for isotopes from lakewater (δ L ; red) and inflow (δ I ; blue) as δ 18 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.