Landscape controls of surface-water/groundwater interactions on shallow outwash lakes: how the long-term groundwater signal overrides interannual variability due to evaporative effects

The spatial and temporal controls on variability of the relative contributions of groundwater within and between flow systems to shallow lakes in the low-relief glaciated Boreal Plains of Canada were evaluated. Eleven lakes located in a coarse glacial outwash, of varying topographic positions and potential groundwater contributing areas, were sampled annually for stable O and H isotope ratios over the course of 8 years. It was demonstrated that landscape position is the dominant control over relative groundwater contributions to these lakes and the spatial pattern of the long-term isotopic compositions attributed to groundwater overrides interannual variability due to evaporative effects. Lakes at low landscape positions with large potential groundwater capture areas have relatively higher and more consistent groundwater contributions and low interannual variability of isotopic composition. Isolated lakes high in the landscape experience high interannual variability as they have little to no groundwater input to buffer the volumetric or isotopic changes caused by evaporation and precipitation. An alternative explanation that lake morphometry (area and volume) control long-term isotopic compositions is tested and subsequently refuted. Landscape position within coarse outwash is a strong predictor for relative groundwater input; however, surface-water connections can short circuit groundwater pathways and confound the signal. A hydrogeological case study for three of the study lakes is used to contextualize and further demonstrate these results.


Introduction
Shallow lakes in the sub-humid climate of Canada's Boreal Plains (BP) exist in a fine balance between precipitation and evaporation. However, in coarse-textured landscapes, the connection to larger-scale groundwater flow systems may buffer these lakes from most precipitation-deficit conditions and can decrease a lake's sensitivity to anthropogenic impacts, both of which are common in the BP region of Canada (Arnoux et al. 2017;Smerdon et al. 2005). Yet, although groundwater influxes to lakes, even at low rates, are often important to both hydrologic and nutrient budgets, the process is poorly understood and often ignored (LaBaugh et al. 1995;Plach et al. 2016;Rosenberry et al. 2015;Schuster et al. 2003).
The BP is characterized by numerous shallow lakes and lake-wetland complexes underlain by thick glacial deposits, which result in interdependent surface-water processes and groundwater flow systems with varying spatiotemporal controls (Lennox et al. 1988;Lissey 1971;National Wetlands Working Group 1988;Winter 1999Winter , 2001. These shallow lakes act as important ecological, biodiversity, and biogeochemical hotspots (Beklioglu et al. 2016;Cheng and Basu 2017), and are critical for continental waterfowl populations (Blancher and Wells 2005;Mack and Morrison. 2006;Slattery et al. 2011), indigenous traditional uses and recreational activities (Parlee et al. 2012). The subhumid climate, where evapotranspiration often equals or exceeds precipitation, results in a delicate water balance with large evaporative demands (Bothe and Abraham 1993;Devito et al. 2005a, b). The BP also experiences multiyear wet and dry cycles, with the potential for increases in frequency and extremes of multidecadal drought cycles with climate change (Bonsal et al. 2013;Ireson et al. 2015;Mwale et al. 2009). During these extended dry periods, some isolated shallow lakes dry out completely, while other lakes, better connected to groundwater sources, are permanent fixtures in the landscape (Smerdon et al. 2005; Thompson et al. 2017). Impacts of climate and geology on surface-water/groundwater interactions are largely unknown and are difficult to predict due to the large interannual variability in precipitation and evaporation and subsequent variability in shallow lake water budgets observed in this region Ireson et al. 2015). The hydrology, biogeochemistry, and ecology of boreal shallow lakes are at an ever-increasing risk to impacts due to anthropogenic activity, such as forestry, oil and gas, agriculture, increased wildfire frequency and magnitude, and climate change (ESTR 2014;Ireson et al. 2015;Schindler 1998). Shallow lakes, especially those with large catchments areas, have been shown to be particularly sensitive to eutrophication, salinization, and alkalinization as anthropogenic and climatological impacts effects these lakes (Tammelin and Kauppila 2018).
The role of groundwater is an important factor in determining lake permanence and resilience to change. Although groundwater flow is much slower than surface flow, the area of interaction between a lake and an adjacent groundwater flow system is large and the flow is usually more uniform, which results in appreciable volumes of water entering and/ or leaving a lake on an annual basis (Freeze and Cherry 1979;Tóth 1999). Understanding and anticipating surfacewater/groundwater interactions and dynamics for these lakes are essential for forecasting water quantity and quality in the BP; however, very few lakes in this region have been gauged, monitored, or sampled. Consequently, techniques and frameworks to evaluate surface-water/groundwater connectivity are needed to better predict lake sensitivity and response to climate change.
There have been many intensive studies in the BP on singular lakes or connected lake systems (e.g., Ferone and Devito 2004;Gibson et al. 2002Gibson et al. , 2015Gibson et al. , 2016Riddell 2008;Shaw and Prepas 1990;Smerdon et al. 2005Smerdon et al. , 2008Thompson et al. 2015); however, these studies have focused on a single lake, deeper lakes, or purely isolated or headwater lakes. The long-term dynamics and patterns of surfacewater/groundwater interactions with shallow BP lakes are not well studied. While it has been shown that there is minimal lateral subsurface hydrologic connectivity through fine-textured mineral substrates in the BP, there is evidence of local to large-scale groundwater connectivity in coarsetextured landscapes Ferone and Devito 2004;Hokanson et al. 2019Hokanson et al. , 2020Smerdon et al. 2005).
Coarser-textured substrates tend to have deeper water tables and higher rates of groundwater recharge to provide base flows during dry periods Haitjema and Mitchell-Bruker 2005;Smerdon et al. 2008).
Extensive research conducted in the glacial outwashes of the American Midwest shows that lake-groundwater connectivity can be dynamic and operate at many different scales resulting in variable groundwater contributions to lakes (Cherkauer and Zager 1989;Jaquet 1976;Kenoyer and Anderson 1989;Krabbenhoft et al. 1990;Krabbenhoft and Webster 1995;LaBaugh 1986;LaBaugh et al. 1997;Siegel and Winter 1980;Winter 1986). Studies in northern Wisconsin showed that landscape position had a strong control on groundwater contribution within the same flow system, where lower lakes had the highest contributions from groundwater and vice versa (Kratz et al. 1997;Webster et al. 1996). Additionally, lakes that received greater contributions from groundwater were more resistant to drought, had higher biodiversity, and higher concentrations of base cations. However, the annual precipitation (i.e., potential groundwater recharge) at these Wisconsin lakes is up to 50% higher and the lakes are much deeper than in the BP. In such humid regions, groundwater flow is characteristically driven by topography and flow is therefore more predictable and temporally stable than in the BP, where the water table is often not a subdued replica of surface topography (Hokanson et al. 2019).
In the subhumid climate of the BP, evaporation is commonly the largest component of a lake's water budget (Devito et al. 2005a, b;Smerdon et al. 2005;Thompson et al. 2015). Accordingly, the onset and duration of evaporation is an important control on lake stage. Precipitation and groundwater inputs are therefore critical fluxes for shallow lake water level maintenance. Although local-scale groundwater flow systems adjacent to lakes can be both temporally and spatially dynamic, lakes positioned low in regional-scale flow systems may have more temporally stable groundwater flow (Webster et al. 1996).
The objective of this study was to evaluate and compare the relative groundwater contributions to 11 lakes located in a coarse-textured glaciofluvial outwash in the subhumid BP to determine the landscape controls of lake-groundwater interactions. Annual mid-summer stable O and H isotope ratios over eight consecutive years were used for the analysis. Specifically, the study aims to determine if (1) the spatial patterning between lake summer stable water isotope ratios was temporally persistent, (2) lake landscape position controls interannual variability of isotopic ratios, and (3) lakes lower in their respective flow system do indeed have greater groundwater contributions. Alternatively, the study (4) tests whether lakes with lower area-to-depth ratios are merely systematically associated with less evaporative influence on the isotopic mass balance leading to similar conclusions regarding groundwater inputs, albeit through a different mechanism. It is expected that isotopic compositions of lakes positioned lower in the landscape will indicate relatively large nonevaporative water inputs from groundwater with low interannual variability, while those of lakes with higher landscape positions will indicate more evaporative effects with notably higher interannual variability.

Study area and hydrogeology
The Utikuma Region Study Area (URSA; 56°N, 115°W) is located 370 km north of Edmonton, Alberta, in the BP ecozone of Canada (Fig. 1, inset). The climate is subhumid with historic annual potential evapotranspiration (ET = 517 mm; Bothe and Abraham 1993) often exceeding the historic average annual precipitation (P = 481 mm; Marshall et al. 1999). The region is characterized by low topographic relief and thick (45-240 m) heterogeneous glacial substrates overlying the Smoky Group, a marine shale of Cretaceous age (Vogwill 1978). This study focuses on 11 lakes on a coarse-textured glaciofluvial outwash that are shallow, cold polymictic lakes, but vary in landscape position, potential groundwater contributing area and lake area.
While all 11 lakes are in the coarse outwash, there is some variation in texture within the outwash and they can be separated into four distinct systems. Three lakes (7, 11, and 19) are located at the transition between coarse outwash deposits and a fine-textured hummocky stagnant ice moraine; they have closed basins with no channelized inflows. Lake 19 is in silt rich glacial fluvial material deposited on a thick lowpermeability clay layer 12 m above the regional water table in the underlying sand aquifer, which is part of the coarse outwash (Hokanson et al. 2019). Lake 19 has been shown to have limited to no local groundwater input and the water budget is controlled by precipitation and evaporation (Riddell 2008). Lakes 7 and 11 have both disintegration moraine and glacial fluvial ice-contact deposits near the basin and in the surface topographic catchments. Based on surficial geology maps (Fenton et al. 2013) and elevation, it is suspected that lakes 7 and 11 are functionally similar to lake 19; however, no deep wells are situated adjacent to the lakes to confirm this.
Lakes 208, 206, and 201 are in a region with clean, wellsorted sands, while lake 206 has first order channelized inflow from a local catchment with mixed disintegration moraine and glacial fluvial deposits. There are no channelized surface connections between lakes 206 and 201 or 208, but diffuse flow through wide peatlands occurs at lake 201 and 208 outlets. Lake 206 is adjacent to fine-textured hummocky stagnant ice moraine deposits; however, boreholes confirm that it is not perched on fine-textured sediment and is well connected to the coarse outwash.
Lakes 1, 17, 16, 5, and 2 are located between the two regions (i.e., the low permeability clays in the eastern region of the study area and the clean sands in the western region) and are characterized by sands and gravels with notable laterally discontinuous silt lenses. No surface outflow or channelized inflow is visible between lakes 1, 17, 16 and 5. Diffuse flow through wide alluvial fens occurs at the outflow area of lake 16. Diffuse flow and occasional channelized flow through a fen wetland connect an upstream lake to the north with lake 1. Surface flow is present from lake 5 to lake 2, and channelized outflow is visible at lake 2. All near surface lateral diffuse inflow or outflow in wetlands and channelized flow have been historically or are currently influenced by beaver activity. Lake 19, and presumably lakes 11 and 7, which are perched on fine textured sediment, contribute recharge to the underlying coarse sand aquifer which connects lakes 17, 16, and 5; however, they are hydraulically separate from the adjacent groundwater flow system.
All the lakes examined are within a 50-km 2 study area and are therefore subject to similar weather and climate conditions . Disregarding the possible effects of wind sheltering due to surrounding vegetation, precipitation and evaporation fluxes have been assumed similar amongst the lakes. The primary dissimilarities between these lakes are therefore variable surface-water/groundwater interactions and lake morphometry (i.e., area and depth). Daily precipitation amounts were collected throughout the study period (2012-2019) using two or three tipping bucket rain gauges adapted for snowfall using an antifreeze reservoir. Lake area was determined from air photos. Lake depth measurements were taken by boat or during ice cover during average climatic conditions (i.e., not during prolonged drought or deluge; Leader 2021). Lake volume was approximated as the product of lake area and average lake depth, which is a reasonable representation for the purpose of comparing lakes within the same region (Hayashi and van der Kamp 2000).
Landscape position was determined for each lake following the approach of Kratz et al. (1997) and Webster et al. (1996), which considers the relative location of each lake within their local hydrologic flow system in addition to the topographic elevation. Given a specific locale, landscape position should provide a method for predicting lake responses to drought, water quality, and ecological viability for aquatic communities (Kratz et al. 1997). Lake elevation and potential groundwater capture areas were determined from a DEM (Montgomery et al. 2019).
Groundwater flow systems were estimated for the study area using regional topographic relief, and further informed by previous studies (Hokanson et al. 2019;Riddell 2008;Smerdon et al. 2005). Potential groundwater capture areas for each lake were determined by delineating flow systems using regional highs and lows. The high permeability substrates and low groundwater recharge typical of this BP outwash result in water tables that may not mimic local topography and extend far beyond the local hillslopes adjacent to each lake . Regional-scale ridges and valleys were used to constrain the groundwater capture areas and local topographic highs (i.e., 5-10-m relief) were disregarded. Delineating groundwater flow systems in high permeability, low relief, locally hummocky regions such as this can be difficult. The groundwater contributing area to lakes in this type of terrain and geology most certainly extends beyond local surface watersheds and changes through time . Consequently, it is assumed that the flow system extents ( Fig. 1) are conservative estimates. Additionally, the focus is on the position of lakes within each flow system and comparisons within flow systems where possible.

Stable O and H isotope ratios in water
As part of the Hydrology, Ecology and Disturbance (HEAD) long-term research and monitoring project in URSA , precipitation (n = 696), lakes (n = 759), and groundwater wells (n = 614) were sampled throughout the year for stable O and H isotope ratios in water from 1999 to 2019. Precipitation samples were collected at three different meteorological stations (Fig. 1). Rain was sampled using collectors specifically built to prevent fractionation, either that of Gröning et al. (2012) or using mineral oil to prevent evaporation from the collector, while snow samples were collected as soon after snowfall as possible as described in Clark and Fritz (1997). Due to the remote nature of the study sites, the rain collectors captured multiple storm events before being sampled and therefore the resultant local meteoric water line (LMWL) is not weighted. Sampled groundwater wells ranged from 0.5 to 5.5 m below the ground surface; details regarding groundwater well installation and sampling are outlined in Hokanson et al. (2019).
For this study, the 11 study lakes were sampled between day-of-year (DOY) 184 and 193 annually from 2012 to 2019 for stable isotope ratios of O and H in water. Water samples were collected at a depth of 30 cm from the lake surface to ensure a representative sample. All lakes are polymictic and well mixed. Each yearly sampling campaign was completed within a 2-or 3-day period. Therefore, on a year-to-year basis, the lakes were "ice-free" and free to evaporate for comparable periods of time each year. Lakes that could be safely accessed by snowmobile (19,206,201,208,17,16,5,2,1) were sampled on 20 March 2020, near the end of the winter lake ice season. During this winter campaign, samples were collected from the top and bottom of the water column as well as near the bottom of the lake ice. Ice thickness was uniform (51.5 ± 5.2 cm) between lakes.
Measured isotope ratios are given in per mil (‰), relative to Vienna Standard Mean Ocean Water (VSMOW), in both a dual-isotope space and as the line-conditioned excess (lc-excess). All isotopic analyses were performed at BASL laboratory, University of Alberta, Edmonton, Alberta, Canada.
The influences of evaporation and sources of water to a lake can be inferred using a secondary isotope-derived parameter, line-conditioned excess (Esquivel-Hernández et al. 2018;Landwehr and Coplen 2006;Sprenger et al. 2017). Line-conditioned excess is a locally relevant metric, compared to the global meteoric water line, which is used in calculating deuterium excess (Dansgaard 1964). Lineconditioned excess is defined as: where a and b are the slope and intercept of the LMWL, respectively (Landwehr and Coplen 2006). The lc-excess describes the offset of a water sample from the LMWL in dual isotope space, where a negative lc-excess value is indicative of evaporative fractionation. If a body of water such as a lake, is initially formed by meteoric water (i.e., precipitation), it will plot directly on the LMWL and have an lc-excess value of zero. If that lake is allowed to freely evaporate, the lighter isotopologues of water will preferentially escape through the process of kinetic fractionation and the lake water will plot below the LMWL and have a negative lc-excess value.
Due to the shallow nature of the lakes studied here, they receive shallow poorly mixed groundwater, the isotopic signature of which can vary temporally and spatially. Additionally, the volume, precipitation input, groundwater outflow, and isotopic composition of boreal lakes may experience high variability due to summer convective storms, seasonal ice cover, and evaporation-induced fractionation; therefore, a transient or steady-state isotopic mass balance was not performed (Gibson et al. 2002;Krabbenhoft et al. 1994). Rather, long-term spatial trends were used to compare differences in isotopic compositions to infer relative groundwater contributions to the study lakes.

Groundwater flow systems
Four separate hydrologic flow systems were delineated for this study. Groundwater flow is from areas and lakes at higher elevations to lakes at lower elevations within the same flow system, without significant groundwater exchange between flow systems. The four flow systems include: (1) isolated lakes 11, 7, and 19, which are at the highest landscape positions, are isolated from groundwater inputs and have no visible surface inputs; (2) the western groundwater flow system, where lake 206 is in the highest position and lakes 201 and 208 are in similar, lower positions; (3) the eastern groundwater flow system, where lake 17 is in the highest position and lake 2 is the lowest; and (4) a central groundwater flow system, where lake 1 is in a low position before the flow systems joins the western system at lake 2 (Table 1).
Conceptualized groundwater contributing areas for each lake are shown in Fig. 1. Within these areas, the gross groundwater flow is assumed to be from high elevations to low elevations. Due to the coarse texture of the outwash sediments, the contributing area for lakes low in the flow system are large and extend well beyond local ridges (Winter 2001). Alternatively, the potential groundwater capture areas are small for lakes located near the upper ends of flow systems. Potential groundwater contributing areas were calculated (Table 1); however, these boundaries extend, for some lakes, into regions with fine-textured sediments (i.e., predominately silts and clays). Due to their low hydraulic conductivity (Hokanson et al. 2019) these fine-textured regions most likely do not generate appreciable groundwater flow when compared to the high-conductivity coarse-textured outwash; areal estimates both with and without estimated areas of fine-textured sediments are shown in Table 1.

Long-term isotopic compositions of URSA lakes and groundwater
The LMWL for URSA derived from the long-term precipitation samples collected from 1999-2019 (Fig. 2) is defined as: The slope of the URSA LMWL is lower than that of the global meteoric water line (GMWL; δ 2 H = 8 · δ 18 O + 10) of Craig (1961) and is similar to that of Mildred Lake (200 km NE of the URSA) in the Athabasca Oil Sands area of Alberta (δ 2 H = 7.20 · δ 18 O -10.3; Baer et al. 2016). Rain at URSA has more isotopic variability than snow, Fig. 1 Ground elevation of the Utikuma Region Study Area (URSA) and its relative location within Canada and North America (inset), with study lake locations and elevations, delineated potential groundwater capture areas and surficial geology type (Fenton et al. 2013). Bold lines indicate limits of separate flow systems (i.e., the poten-tial groundwater capture area for the lowest lake in that system); finer lines indicate potential groundwater capture area for the lakes positioned higher in each system. Cross section A-A′ is presented in Fig. 2c with the δ 18 O of rain ranging from -22.2 to -10.2 ‰ and snow ranging from -26.0 to -19.3 ‰, which represent the 5th and 95th percentiles, respectively.
Lake water samples at URSA, which includes samples from over 30 lakes across URSA sampled over the same period (i.e., 1999-2019), plot on a local evaporation line (LEL), defined by: The long-term LEL has a slope of 4.41, similar to studies in nearby regions in central and northern Alberta (Gibson et al. , 2019. The intersection of the LEL and the LMWL can be used as an estimate of the weighted mean isotopic composition of annual precipitation, which is often assumed to be the isotopic composition of input to the lake (Gibson et al. 1993;Paulsson and Widerlund 2020). Here, the intersection falls approximately between the composition of rain and snow, tending towards rain (δ 18 O ≈ -18 ‰; δ 2 H ≈ -142 ‰; Fig. 2a). However, this may not be true for nonheadwater lakes where the primary input could be groundwater. The slope of the LEL showed high interannual variability and ranged from as low as 3.4 (R 2 = 0.89; n = 42) in 2017 and as high as 4.7 (R 2 = 0.83; n = 52) in 2012 (Fig. 2a inset). This variability usually reflects the interannual variability of temperature, humidity, and wind over the evaporation season (Gibson et al. 1993) and may be influenced by strong seasonality of evaporation and ice coverage of lakes (Gibson et al. 2008). Just as the slope of the LEL varied year-to-year, so did the intercept with the LMWL as the weighted mean annual isotopic composition of precipitation change yearto-year. The intercept varied by approximately 4 ‰ and ranged from -19 to -15 ‰ in δ 18 O space.
For the most part, the isotopic ratios of the groundwater near the study lakes in the glacial outwash (Fig. 2b) plot directly on the LMWL, meaning the groundwater is of meteoric origin and has not undergone evaporative fractionation. There are, however, some exceptions where wells intersect groundwater that was recharged directly from an adjacent lake and shows a fractionated signature and plots on the LEL; these cases are discussed in a later section. The isotopic composition of the groundwater at the outwash tends more towards that of snow and shoulder season (i.e., early spring and late fall) rain.
While the groundwater samples plot on the LMWL, there is considerable spread. This is evidence of the local and episodic nature of groundwater recharge typical of the BP. Most precipitation (50-60%) falls during summer months when vegetation is active and potential evapotranspiration is at its peak; therefore, most groundwater recharge in the BP occurs outside of the growing season and is the product of snow melt and late season storms (Smerdon et al. 2008; (3) Thompson et al. 2015). Deeper groundwater, representing regional groundwater flow systems would be old and well mixed, and thus temporally uniform (Tóth 2009). Studies which utilize lake isotopic mass balances or end-membermixing analyses yield the best results when groundwater is uniform and well mixed, usually clustered at the intersection of the LEL and the LMWL (e.g., Arnoux et al. 2017;Gibson et al. 2019;Krabbenhoft et al. 1990;Petermann et al. 2018); however, this is often assumed, yet may not be true. Groundwater in this study area, especially shallow groundwater, is not well mixed and varies isotopically in both time and space (Fig. 2b). In many studies, the groundwater is either not sampled or not sampled extensively, either spatially or temporally (e.g., Jones and Imbers 2010; Shi et al. 2017;Zuber 1983).

lc-excess of study lakes
Due to the high interannual variability of the URSA LEL (slope and intercept) and groundwater isotopic ratios in the outwash, it is difficult to determine the starting composition or the composition of the inflow for the study lakes on an annual or subannual basis. While this variability makes creating an accurate mass balance difficult, the relative influence of evaporation in lakes can be inferred by using lc-excess (Eq. 1; Landwehr and Coplen 2006), where a negative value indicates evaporative fractionation. Because both groundwater and precipitation at URSA plot along the LMWL, the lc-excess parameter is more informative than δ 2 H or δ 18 O alone since it describes the offset of a water sample from its original meteoric composition, regardless of the variability due to seasonality along the LMWL (Esquivel-Hernández et al. 2018; Sprenger et al. 2017). Lakes with the smallest offset from the LMWL will have the highest (i.e., closest to zero) lc-excess values and represent lakes whose water balance is least affected by enrichment due to evaporative fractionation. These lakes are generally more flushed by meteoric-sourced water (i.e., groundwater inflow and outflow). Conversely, lakes with large offsets from the LMWL will have low, more negative lc-excess values and will represent lakes whose water balance is dominated by evaporative fractionation and have little to no external sources of water . While there are characteristic temporal changes in lake lc-excess due to the annual variations in the balance between precipitation and evaporative fractionation (Gibson et al. 1993), groundwater input can buffer these changes, where lakes with little or no groundwater input rely entirely on precipitation to mitigate volumetric losses and isotopic enrichment due to evaporation. Variations in lc-excess of the 11 study lakes in the coarse-textured glaciofluvial outwash are shown in Figs. 3 and 4. As landscape position decreases and groundwater contributing area increases, the lc-excess values generally increase (Table 1). This is evident within both the western (lakes 206, 201, and 208) and eastern (lakes 17, 16, 5, and 2) flow systems as a function of landscape position, and between flow systems as a function of groundwater contributing area (Fig. 5a). This increase is due to progressively larger volumes of groundwater, of meteoric origin, being contributed to the lakes through the permeable substrate of the outwash (Cheng and Anderson 1994;Smerdon et al. 2005;Townley and Davidson 1988;Winter et al. 2003). The isolated lakes (11, 7, and 19) had a much wider distribution of lc-excess and showed the greatest evaporative enrichment, additionally they had lc-excess values as low as -41.9, -49.9, and -38.8‰, respectively. Within each flow system (Table 1; Fig. 4) there is notable temporal persistence in the relationships of annual lake lcexcess values. For every year in the eastern flow system, the relative order of lake lc-excess is in ascending order: lakes 17, 16, 2, and 5. In the western flow system, the lc-excess of lake 201 is consistently higher than lake 206, whereas lakes 208 and 201 are very similar and track together over time. Lake 1, which is in its own central flow system, but similar landscape position to lakes 2 and 5, regularly has an lc-excess between both lakes.
The persistence in the pattern of lake lc-excess is also present during the winter sampling in March 2020 (Fig. 3). The lakes were sampled near the end of the ice-over period to test whether the spatial pattern of lc-excess values described in the preceding was present after the formation of ice and whether continued influxes of groundwater without the presence of evaporative fractionation would further raise the lc-excess values of the lake water. Heavy isotopologues of water are preferentially incorporated into ice as it forms, thereby decreasing the δ 2 H and δ 18 O values of the underlying lake water (i.e., resulting in increased lc-excess values), mimicking the effect of groundwater inputs to the lake. The shallow lakes, except for possibly the deepest, lake 17, can be assumed well mixed during the ice-off season; however, isotopic stratification through fractionation may be common during winter ice-over. The lc-excess of the winter lake ice was similar to or lower than under-ice water samples, as expected (Fig. 3). Most lakes did not show significant stratification, although lakes 16 and 2 exhibited notable differences in lc-excess values between the top and bottom of the water column. Lake 1 showed a slight difference. Due to the fractionation that takes place during ice formation it is expected that water just below the ice would be isotopically lighter than deeper lake water but, interestingly, the water sampled at the bottom of the lakes is even lighter. If the stratification were purely due to ice creation, isotopically lighter water (i.e., less negative lc-excess) would be located at the top of the column and isotopically heavier water (i.e., more negative lc-excess) at the bottom (Krabbenhoft et al. 1990); however, the opposite trend is seen at lakes 16, 2, and 1. Both the "reverse stratification" at lakes 16, 2, and 1 and the lack of stratification at the remaining lakes could be the result of isotopically light groundwater discharging to the lake throughout the winter. More rigorous sampling is required to confirm this hypothesis.
The region has low annual precipitation (~444 mm year −1 ) with significant variability between years, which in turn affects the annual isotopic signatures and water balances of these BP shallow lakes (Fig. 4). The greatest disparity between interannual variability in lc-excess can be seen  (11, 7, and 19) with the more connected lakes that are part of larger flow systems. Isolated lakes had the largest interannual variability, with lc-excess values ranging from -50 to -12 ‰ over the course of the study period, and had an average interquartile range (IQR) of 12.7 ‰. The annual summer lc-excess values of the isolated lakes had a strong linear correlation with the cumulative annual precipitation (i.e., precipitation accumulated over the previous year from the sampling date) (R 2 = 0.69; p < 0.001). Without groundwater to buffer them the isolated lakes are virtually exclusively controlled by precipitation and evaporation, and due to their shallow nature and overall low volume, small amounts of precipitation would have large controls over the isotopic mass balance. These findings support previous studies that show that Lake 19, which is located high in the landscape, had only ~20 mm of shallow lateral inflow during the 2005 and 2006 hydrologic years (Riddell 2008). In contrast, the more connected lakes had less interannual variability and the annual values were not correlated with precipitation (R 2 = 0.18). Lakes in the eastern flow systems (17, 16, 5, 2, 1) appear to have similar interannual variability (average IQR = 6.3 ‰) over the study period, and do not appear to vary systematically with topographic position. The western lakes (206,201,208) have the least interannual variability (average IQR = 3.9 ‰). This is likely due to the difference in geology between the western lakes, which are in an area with clean, uniform, high conductivity sands that promotes steady hydraulic conditions. The eastern lakes, in contrast, have notable heterogeneity in the form of sands and gravels interbedded with silt lenses, the latter of which can result in transient lateral groundwater connectivity.

Possible effects of lake morphometry on lc-excess
Lake volume has a direct effect on the isotopic mass balance, and therefore on the amount of groundwater or precipitation necessary to maintain lake levels and less evaporative isotopic compositions. Deeper lakes are less sensitive to volumetric and concentration changes for a given evaporation or precipitation event when compared to shallower lakes with the same surface area. As an alternative explanation Well LEL (long-term average) Fig. 2 Dual isotope plot of Utikuma Region Study Area (URSA) precipitation and lake waters with (inset) annual local evaporation lines (LEL) from 1999-2019; b dual isotope plot of URSA groundwater sampled from 1999-2019 and c a generalized geologic cross-section along transect A-A′ ( Fig. 1) with groundwater lc-excess (sampled from 2000 to 2019) values shown above. Each set of box and whiskers represents one groundwater well or a cluster of wells. The cross section shows the lake and well locations that lie on the transect (black numbers indicate lake ID on the cross section; gray numbers indicate lakes off the cross section). Vertical lines represent boreholes; see Hokanson et al. (2019) for more details ▸ to groundwater inputs being the primary mechanism for lake maintenance, deeper, larger lakes could simply be less sensitive to evaporative fractionation and merely be maintained by precipitation in this landscape. However, this study has both a large lake (lake 17; ~6 × 10 6 m 3 ) with very little groundwater contributing area and very low lc-excess (-26.7 ‰), and a small lake (lake 1; ~2.5 × 10 5 m 3 ) with a large groundwater contributing area and relatively high lc-excess (-15.9 ‰). If lake morphometry were a controlling factor (i.e., lakes with larger areas and shallower depths are associated with greater evaporative influences) there would be a negative correlation between lc-excess and the lake area/ depth ratio (Yang et al. 2019). At the study lakes where it was hypothesized that groundwater was a controlling factor for lc-excess, the opposite is true-there is a weakly positive correlation (R 2 = 0.25). The pattern of lc-excess appears to be independent of lake morphometry (Fig. 5b,c).

Competing roles for groundwater and surface water
There are varying degrees of interaction between surface water and groundwater in this glaciofluvial outwash. Due to the local subhumid climate, lakes act as 'evaporation windows' on groundwater and can influence isotopic compositions in both groundwater outflow and downstream lakes by reintroducing fractionated (i.e., lake-affected) water back into the flow system, which has the potential to then discharge to down-gradient lakes (Gat and Bowser 1991). This is evident in both groundwater and lake water samples at the URSA outwash. The groundwater is mostly isotopically uniform and similar to that of meteoric waters (Fig. 2). However, wells on the southwest side of lake 206 (green symbols on Figs. 2b) and the well between lakes 5 and 2 (black symbols on Figs. 2b), show evaporatively fractionated water.
Groundwater sampled adjacent to lake 206 is isotopically heavier when compared to other groundwater, although less so than lake 206 itself. This is a result of water in lake 206 mixing with ambient groundwater recharged from meteoric water. There is also evidence for lake water being affected by other 'upstream' lakes, in a "chain of lakes effect" (Gat and Bowser 1991). Lake 2 has a lower lc-excess value relative to its landscape position. Although it is in the lowest landscape position for its respective flow system, lake 2 has a lower lc-excess than lake 5, which is up-gradient and has a correspondingly smaller potential groundwater contributing area (Fig. 3). Lake 2 is one of only two lakes in this study to have surface inflow, where a semipermanent stream flows from lake 5 to lake 2 most years (Table 1; Fig. 1). Both groundwater and surface water flowing from lake 5 to lake 2 are isotopically heavier due to the fractionation that took place in lake 5. Consequently, while the sources of inflow to the remainder of the lakes are either precipitation or groundwater, both of which have a meteoric origin, lake 2 receives both slightly fractionated groundwater and fractionated lake water from lake 5.
Synoptic isotope sampling has become a popular method of characterizing major hydrologic fluxes, sources, and fates. The data presented here, however, demonstrate the need for both surface and groundwater processes to be contextualized together and that their interactions may vary in both space and time.
A hydrogeology case study: lakes 17, 16, 5 Previous works by Smerdon et al. (2005Smerdon et al. ( , 2008Smerdon et al. ( , 2012 have provided a comprehensive hydrogeological model and understanding of the hydrologic controls on lakes 17, 16, and 5. These studies included intensive monitoring of hydrometric  (Table 2). Groundwater moves from southeast to northwest at an average horizontal gradient of 0.002 ( Fig. 6; Smerdon et al. 2008). For lake 16, and most other lakes in the Boreal Plain, precipitation and evaporation are the dominant influx and outflux, respectively, on an annual basis (Ferone and Devito 2004;Smerdon et al. 2005;Thompson et al. 2017). It was shown for lake 16 for two consecutive years that the groundwater components were consistent, contributing 230 and 222 mm year −1 during the 2002 and 2003 hydrologic years, respectively. Groundwater discharge to the lake was temporally consistent throughout the year and was the dominant influx between November and April, when evaporation is limited due to low temperatures and/or ice cover. Lake 16 is 2.4 m lower than lake 17 and 1.8 m higher than lake 5. The stages at lakes 17, 16, and 5 showed progressively less variability as their landscape position decreased.
Lake 17 was shown to have little to no input from groundwater sources and to be sensitive to annual precipitation as a primary input, while lakes 16 and 5 are progressively lower in the topographic gradient and therefore capture more groundwater flow. Lake 17 was most sensitive to annual precipitation and evaporation. Figure 6 shows the groundwater flow directions around lake 16, illustrating the position of the lakes within a larger flow system (Smerdon et al. 2008). The isotope data for lakes 17, 16, and 5 presented here expand on and fully support the findings of Smerdon et al. (2005Smerdon et al. ( , 2008, where lc-exess is progressively less negative as the landscape position decreases (i.e., from lake 17 to 16 to 5) in one groundwater flow system (Fig. 3). This pattern, irrespective of interannual climate variation, is exceedingly stable through time (Fig. 6). These principles of lake landscape position, potential groundwater contributing area, and spatial changes in isotopic concentrations are expanded on, both spatially and temporally, and supported by the isotopic Relationships between median July lc-excess anda groundwater contributing area and b-c lake morphometry (area and depth). Lakes considered isolated are not included due to the wide range of summer lc-excess values. Groundwater contributing area only consid-ers coarse surficial geology (Table 1). Marker size is proportionate to lake volume (a) and lake area (b). Lake volume is approximated by multiplying lake area by average depth data presented for the nonisolated lakes (17,16,5,2,1,208,201,206) in this study.

Conclusion
Eleven shallow boreal lakes were sampled annually for 8 years. Stable O and H isotope ratios indicate that lakes with lower landscape positions, and associated larger potential groundwater-contributing area, had relatively larger groundwater inputs. This is the first study to show that the spatial pattern of long-term groundwater signals amongst lakes in the same flow system overrides the seasonal evaporative effects on stable water isotopic compositions of these shallow boreal lakes. In other words, the spatial variability of lc-excess between lakes within the same groundwater flow systems persisted through time, despite interannual variability due to climate-driven evaporative losses. This work demonstrates the feasibility of using lc-excess to qualitatively compare water sources of lakes within flow systems; future work in this complex environment should focus on performing isotopic mass balances to quantify groundwater inputs and impacts on lc-excess of lake water over time. Isolated lakes high in the landscape experience much higher interannual variability because they have little to no groundwater input to buffer the volumetric or isotopic changes due to evaporation and precipitation. Landscape position within coarse outwash landscapes is a strong predictor for relative groundwater input; however, surface-water connections can short circuit groundwater pathways and confound the signal. Due to the highly heterogeneous glacial substrates present in the Boreal Plain, it can be difficult to interpret between flow systems; a conceptual understanding of the system's surface and subsurface hydrology is therefore required. A strong conceptual model and thorough understanding of potential the LMWL (dashed line) and long-term LEL (solid line) for reference. The further along the LEL (i.e., offset from the LEL-LMWL intersection) the data point lies, the more negative the corresponding lc-excess value water sources and fates is required to explain the spatial heterogeneity in groundwater and surface-water isotopic compositions at local and regional scales.
Funding We acknowledge the support by Natural Sciences and Engineering Research Council of Canada (NSERC) grants to KJD from a Collaborative Research and Development Grant (CRDPJ477235-14) with industry partners Syncrude Canada and Canadian Natural Resources, and funding from forWater, a pan-Canadian strategic research network awarded to KJD and CAM) and numerous other partners (www. forwa ter. ca). The digital elevation model was produced from combined LiDAR point cloud datasets collected by the authors (2015)  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.