Assessing the influence of lake and watershed attributes on snowmelt bypass at thermokarst lakes

Snow represents the largest potential source of water for thermokarst lakes, but the runoff generated by snowmelt (freshet) can flow beneath lake ice and via the outlet without mixing with and replacing pre-snowmelt lake water. Although this phenomenon, called “snowmelt bypass”, is common in ice-covered lakes, it is unknown which lake and watershed properties cause variation in snowmelt bypass among lakes. Understanding the variability of snowmelt bypass is important because the amount of freshet that is mixed into a lake affects the hydrological and biogeochemical properties of the lake. To explore lake and watershed attributes that influence snowmelt bypass, we sampled 17 open-drainage thermokarst lakes for isotope analysis before and after snowmelt. Isotope data were used to estimate the amount of lake water replaced by freshet and to observe how the water sources of lakes changed in response to the freshet. Among the lakes, a median of 25.2 % of lake water was replaced by freshet, with values ranging widely from 5.2 % to 52.8 %. For every metre that lake depth increased, the portion of lake water replaced by freshet decreased by an average of 13 %, regardless of the size of the lake's watershed. The thickness of the freshet layer was not proportional to maximum lake depth, so that a relatively larger portion of pre-snowmelt lake water remained isolated in deeper lakes. We expect that a similar relationship between increasing lake depth and greater snowmelt bypass could be present at all ice-covered open-drainage lakes that are partially mixed during the freshet. The water source of freshet that was mixed into lakes was not exclusively snowmelt but a combination of snowmelt mixed with rain-sourced water that was released as the soil thawed after snowmelt. As climate warming increases rainfall and shrubification causes earlier snowmelt timing relative to lake ice melt, snowmelt bypass may become more prevalent, with the water remaining in thermokarst lakes post-freshet becoming increasingly rainfall sourced. However, if climate change causes lake levels to fall below the outlet level (i.e., lakes become closed-drainage), more freshet may be retained by thermokarst lakes as snowmelt bypass will not be able to occur until lakes reach their outlet level.

1 Introduction

In the continuous permafrost zone of the Arctic, regions with thermokarst lakes have formed where ice-rich permafrost has thawed and the ground surface has subsided. Thermokarst lakes typically range from 1 to 5 m in depth and from 0.01 to 1000 ha in area, can cover over 25 % of the land area (Grosse et al., 2008; Burn and Kokelj, 2009; Turner et al., 2014; Farquharson et al., 2016) and mostly formed during a brief warm period following the last deglaciation of the Northern Hemisphere (Brosius et al., 2021). Comparison of aerial photography from the mid-1900s with more recent satellite imagery has revealed both increases and decreases in thermokarst lake area and number (Smith et al., 2005; Plug et al., 2008; Marsh et al., 2009; Jones et al., 2011; Finger Higgens et al., 2019). These changes are partially attributed to shifting thermokarst lake water balances: increased air temperatures (Woo et al., 2008), longer ice-free seasons (Surdu et al., 2014; Arp et al., 2015), permafrost thaw (Walvoord and Kurylyk, 2016), and shrub expansion leading to increased transpiration (Myers-Smith et al., 2011) and interception (Zwieback et al., 2019), all causing less inflow and more water to evaporate from thermokarst lakes. In contrast, increasing precipitation can lead to more inflow to lakes, offsetting any rise in evaporation, interception and transpiration (Walsh et al., 2011; Stuefer et al., 2017; Box et al., 2019; MacDonald et al., 2021), while shrub expansion can also increase snow accumulation in lake watersheds, resulting in more snowmelt runoff to lakes (Turner et al., 2014; MacDonald et al., 2017). Increased rainfall has also been linked to decreases in lake surface area because lakes are more likely to experience rapid drainage due to permafrost thaw during wet years (Webb et al., 2022).

Runoff generated by snowmelt in lake watersheds represents a large potential water source for lakes, as snowfall comprises 40 % to 80 % of total precipitation in the Arctic (Bintanja and Andry, 2017). When snow melts in spring, the volume of snowmelt-driven runoff flowing into lakes (freshet) generally results in the highest lake levels of the year (Woo, 1980; Roulet and Woo, 1988; Hardy, 1996; Pohl et al., 2009). When freshet is low, thermokarst lakes are prone to desiccation (Marsh and Bigras, 1988; Marsh and Lesack, 1996; Bouchard et al., 2013). It is a reasonable expectation that lakes which receive more freshet will also contain more freshet by the end of the snowmelt if they remain below their outlet level (i.e., closed-drainage lakes). However, for non-bedfast ice-covered lakes at or near their outlet level (i.e., open-drainage lakes), freshet may flow into and out of a lake without mixing with and replacing the pre-freshet lake water, resulting in “snowmelt bypass” (Bergmann and Welch, 1985) (Fig. 1). While lake ice inhibits wind-driven mixing of lake water, the cooler, less dense freshet (∼0 ^∘C) cannot mix with the deeper, warmer and denser lake waters (<4 ^∘C). As a result, freshet water will flow into and out of an open-drainage lake without replacing the deeper, pre-snowmelt lake water until vertical mixing within the lake begins, which is initiated by the warming of lake waters from solar radiation penetrating through snow-free ice and wind-driven mixing after the lake becomes ice-free (Cortés and MacIntyre, 2020). Snowmelt bypass is a common occurrence that has been observed in a wide variety of ice-covered open-drainage lakes around the world (Henriksen and Wright, 1977; Jeffries et al., 1979; Hendrey et al., 1980; Bergmann and Welch, 1985; Schiff and English, 1988; Edwards and McAndrews, 1989; Cortés et al., 2017).

Figure 1 A conceptual cross section of an open-drainage lake when freshet has begun. Freshet initially flows into the lake at the edge where lake ice has melted. A layer of snowmelt runoff mixed with lake water then remains buoyant on top of the warmer lake water before flowing through the outlet (i.e., “snowmelt bypass”). Limited mixing occurs due to density differences between runoff and deeper lake water and the lack of wind-driven mixing due to the presence of lake ice.

Although previous studies have established the mechanisms and conditions that cause snowmelt bypass, no studies have examined how lake and watershed characteristics affect snowmelt bypass. Given that snowmelt bypass depends on the mixing conditions under lake ice, we hypothesize that lake and watershed characteristics that impact lake mixing may cause variability in snowmelt bypass among lakes in a given region. Understanding the factors that influence the amount of freshet retained by thermokarst lakes is important because of subsequent influence on lake water balance, pH, nutrient composition, and suspended sediment, among other limnological variables (Henriksen and Wright, 1977; Marsh and Pomeroy, 1999; Finlay et al., 2006; Turner et al., 2014; Balasubramaniam et al., 2015).

In this study, we determine factors influencing the magnitude of snowmelt bypass for 17 open-drainage thermokarst lakes in the lake-rich tundra uplands east of the Mackenzie Delta in the Northwest Territories, Canada, during the freshet of 2018. This area contains thousands of thermokarst lakes that constitute up to 25 % of the landscape surface area and have changed in area and number during the past several decades in response to changing precipitation and permafrost thaw (Plug et al., 2008; Marsh et al., 2009). We used isotopic composition of lake water and precipitation from before and after snowmelt to estimate the proportion of lake water replaced by freshet during spring 2018 and evaluated relationships with lake and watershed characteristics. We selected lake and watershed characteristics that had the potential to impact under-ice mixing through their influence on the water temperature profile (e.g., lake depth) or by displacing pre-snowmelt lake water (e.g., watershed area). Isotope tracers were also used to assess whether the freshet is sourced solely from snowmelt or whether other water sources contributed to freshet. Future assessments of hydrological and biogeochemical properties of thermokarst lakes can use the lake and watershed attributes we identify to affect snowmelt bypass and lake water sources to inform their results, given the distinct biogeochemical properties of freshet runoff (Finlay et al., 2006; Balasubramaniam et al., 2015) and the influence of snowmelt bypass on the amount of freshet runoff retained by lakes.

2 Study area

The 17 study lakes are situated in the taiga–tundra uplands east of the Mackenzie Delta, in the northwestern region of the Northwest Territories, Canada (Fig. 2). The landscape is comprised of rolling hills and is strongly influenced by permafrost thaw, as evidenced by the thousands of thermokarst lakes which formed between 13 000 and 8000 years ago (Rampton, 1988; Burn and Kokelj, 2009) that are typically 2–4 m in depth with a surface area from 10 to 1000 ha (Pienitz et al., 1997). The study lakes are situated along a ∼70 km stretch of the Inuvik–Tuktoyaktuk Highway north of the town of Inuvik (Fig. 2). The average area of the lakes is 14.2 ha (0.9–90.5 ha), and the average maximum depth is 2.2 m (1.0–4.1 m) (Table 1). All lakes have a defined outlet channel observed to be active during the spring melt, thus classifying them hydrologically as open-drainage, and many lakes have defined channelized inflows from their watersheds in the form of small streams or ice-wedge polygon troughs.

Figure 2 Blue circles indicate lakes that were sampled before and after snowmelt in 2018. Tundra uplands are in white, while the Mackenzie Delta is in grey. In the inset, key locations near the Trail Valley Creek field station are shown.

Table 1 Lake and watershed properties for sampled lakes. Lake locations are shown in Fig. 2.

Soils in the region have evolved from fine-grained morainal tills, ice-contact sediment, and lacustrine deposits (Rampton and Wecke, 1987). Subsurface flow is conveyed efficiently by a network of interconnected peat channels 0.3–1.0 m wide that lay between mineral earth hummocks (Quinton and Marsh, 1998). Lake watersheds contain tall-shrub (>1 m), low-shrub (∼0.5 m), and shrub-free land cover types comprising lichen, moss, and tussocks (Lantz et al., 2010; Grünberg et al., 2020). Mean annual air temperature in Inuvik is −8.2 ^∘C and mean annual precipitation is 241 mm, of which 66 % is snow, based on 1981–2010 climate normals (Environment and Climate Change Canada, 2019). Snowmelt usually begins in mid-May, and lakes typically become ice-free in June and freeze up in mid-October (Burn and Kokelj, 2009).

The 2018 snowmelt season was typical in comparison to recent decades. End-of-winter snow surveys conducted in the 58 km² watershed of Trail Valley Creek in 2018 (Fig. 2) recorded an average snow water equivalent (SWE) of 141 mm, close to the average SWE of 147±35 mm based on surveys from 1991 to 2019 (Marsh et al., 2019). At Trail Valley Creek, snowmelt began about 1 May, with snow-free areas beginning to appear by 8 May, while only the remnants of large snow drifts remained by 3 June. Lake ice near Trail Valley Creek became snow-free by 10 May, and lakes became completely ice-free on 14 June. The mean air temperature at Trail Valley Creek during the sampling period from 26 April to 15 June 2018 was 0.4 ^∘C, which was cooler than the average of 1.7 ^∘C during 1999–2019 (Fig. 3). Air temperatures roughly followed the average minimum and maximum daily air temperatures, with some temporal variability which can be expected for any given year. Maximum daily air temperatures were mostly above 0 ^∘C after 8 May, which was similar to the average timing of the first above-0 ^∘C day during 1999–2019 (Fig. 3).

Figure 3 Maximum and minimum daily air temperature in comparison to the average and extreme values at the Trail Valley Creek field station for the period of 1999–2019.

3 Methods

3.1 Lake water and precipitation sampling for isotope analysis

Lake water samples for stable isotope analysis were first collected from the 17 study lakes while they were ice-covered (26 April–1 May 2018) and again soon after the lakes became ice-free (15 June 2018). Pre-snowmelt water samples were obtained from a hole augured through the ice near the centre of each lake. These water samples were taken 10 cm below the water surface in the augured hole. Lake depth, snow depth on the ice, and ice thickness were recorded at the same time water samples were collected. Snow depth was typically very uniform across the lake ice. Water samples were then collected post-snowmelt at the shore of each lake shortly after the lakes became fully ice-free. Isotope data were then used to estimate the portion of lake water that was replaced by freshet between the two sampling dates.

To estimate the local meteoric water line (LMWL) and the average isotopic composition of precipitation (δ_P) in the study region, which are useful references for the interpretation of lake water isotopic compositions, samples of end-of-winter snow on the ground in April 2018 and rainfall for the period May to September 2018 were obtained. Snow samples (n=11) were collected from the study area by taking a vertical core of snow using a tube, completely melting the snow in a sealed plastic bag, and then filling a sample bottle with the meltwater. Rainfall (n=13) was collected between May and September in Inuvik using a clean high-density polyethylene container, which was then transferred to a sample bottle shortly after the rain had stopped. The midpoint between the average isotopic composition of snow samples and rain samples was used to define δ_P. All samples were collected in 30 mL high-density polyethylene bottles and were measured using an isotope analyzer to determine the ratio of ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ and ${}^{\mathrm{2}}\mathrm{H}{/}^{\mathrm{1}}\mathrm{H}$ in each sample. Isotope concentrations were measured using a Los Gatos Research (LGR) Liquid Water Isotope Analyser, model T-LWIA-45-EP at the Environmental Isotope Laboratory at the University of Waterloo. The instrument was calibrated using Vienna Mean Standard Ocean Water (VSMOW) and Vienna Standard Light Antarctic Precipitation (VSLAP) standards provided by LGR. Calibration of the instrument was checked during the analysis using the VSMOW and VSLAP standards. Isotopic compositions are expressed in standard δ notation, such that

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}_{\mathrm{sample}}={\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{VSMOW}}}}-\mathrm{1},\end{array}$$

where R represents the ratio of ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ or ${}^{\mathrm{2}}\mathrm{H}{/}^{\mathrm{1}}\mathrm{H}$. Every fifth sample was analyzed a second time to determine the analytical uncertainties, which were ±0.1 ‰ for δ¹⁸O and ±0.6 ‰ for δ²H, calculated as 2 standard deviations away from the difference between the duplicate samples. All isotope data from lakes are presented in Table 2.

Table 2 Snow and ice thickness, isotope, δ_I, and lake water replacement values for all lakes. Ice-corrected δ_L values were used to calculate ice-corrected δ_I values following Yi et al. (2008).

3.2 Estimating the replacement of lake water by freshet and lake source waters

The percentage of a lake's volume that has been replaced by a given water source can be estimated as follows:

$$\begin{array}{}\text{(2)}& \mathit{\%}\text{lake}\text{water}\text{replaced}={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{L}\text{-}\mathrm{Post}}-{\mathit{\delta }}_{\mathrm{L}\text{-}\mathrm{Pre}}}{{\mathit{\delta }}_{\mathrm{I}\text{-}\mathrm{Post}}-{\mathit{\delta }}_{\mathrm{L}\text{-}\mathrm{Pre}}}}\cdot \mathrm{100},\end{array}$$

where δ_L-Pre is the isotopic composition of lake water before snowmelt begins, δ_L-Post is the isotopic composition of lake water after snowmelt is complete, and δ_I-Post is the isotopic composition of the source water post-snowmelt. Application of this equation assumes minimal to no change in volume, which is reasonable given the lakes we sampled are all open-drainage. We also assume that lakes are well mixed as they become ice-free when δ_L-Post water samples were collected. Water temperature measurements at Big Bear Lake (Fig. 2), previous water temperature profiles made at lakes within 10 km of the Inuvik–Tuktoyaktuk Highway and lake temperature modelling using FLake-online (Kirillin et al., 2011) all suggest that lakes were well mixed at the time of δ_L-Post sampling on 15 June 2018 (Appendix A).

We calculated δ_I following the coupled isotope tracer approach outlined by Yi et al. (2008), using an isotope framework based on 2017 air temperature and humidity data for the typical ice-free period (15 June–15 October) collected at the Trail Valley Creek meteorological station located 45 km NNE of Inuvik (Fig. 2). Data from 2017 were used because it was the most recent period when lakes were exposed to meteorological conditions for an entire open-water season. The coupled isotope tracer approach assumes all lakes under the same meteorological conditions will evolve towards the same isotopic composition (δ^*, the isotopic composition of a lake at the moment of desiccation) as lakes evaporate along lake-specific evaporation lines. These lake-specific evaporation lines are defined by extrapolating from δ^* through δ_L until intersection with the local meteoric water line, which is used to estimate δ_I (Fig. 4). We calculated δ_I for pre-snowmelt and post-snowmelt isotopic composition of lake water to identify whether the isotopic composition of the source water changed after freshet. The percentage of lake water replaced was calculated using both δ¹⁸O and δ²H using Eq. (2), and average values are reported. The difference obtained using the two isotopes in the estimate of the percentage of lake volume replaced by runoff was minimal, averaging 1.8 % and ranging from 0.6 % to 3.7 %.

A sensitivity analysis was performed to evaluate uncertainty in the δ_I values and subsequent percentage of lake water replaced by runoff. Confidence in the interpretation of δ_I values with respect to rainfall- or snowmelt-sourced waters depends on an accurate estimate of δ_P, which was determined using the average δ_Rain and δ_Snow values. Also, the calculation of percentage lake water replaced by runoff is sensitive to changes in δ_I, which is sensitive to the average δ_Rain value because this parameter is used to determine δ_As (Eq. B5), which is subsequently used to determine δ^* (Eq. B1). Since there was some variability in the δ_Rain and δ_Snow values from samples we collected, we tested the sensitivity of our estimates of δ_I and the percentage of lake water replaced by runoff to variation in average δ_Rain and δ_Snow values. We calculated the probable upper-bound and lower-bound limits of δ_Rain and δ_Snow values by adding and subtracting the standard error from the average of δ_Rain and δ_Snow values (Appendix B). Upper- and lower-bound cases were propagated through the isotope framework to calculate upper- and lower-bound δ_I and percentage lake water replacement values to evaluate whether the standard error caused enough deviation to meaningfully change the results. Overall, differences between upper- and lower-bound δ_I and percentage lake water replacement values were minimal (Table B2, Fig. B1). Details of the equations and variables used in the isotope framework and the sensitivity analysis are given in Appendix B.

Figure 4 A hypothetical change in the isotopic composition of lake water from pre-snowmelt to post-snowmelt is shown. A visualization of how δ_I is calculated for an individual lake using a lake-specific evaporation line for both pre-snowmelt and post-snowmelt is also shown, where each lake's evaporation line (dashed line) extrapolates from δ^* through δ_L until intersecting the local meteoric water line to give δ_I. The local evaporation line (LEL) is defined by the line between δ_P and δ^* (not shown).

As ice forms and preferentially incorporates water containing the heavy isotopes ¹⁸O or ²H, the lake water beneath the ice becomes increasingly depleted in ¹⁸O and ²H. Consequently, the water samples we collected pre-snowmelt were systematically isotopically depleted relative to pre-freeze-up lake water, and the magnitude of depletion depends on the fraction of lake water that had frozen into ice. We corrected δ_L-Pre for the fractionation of freezing water into ice using an equation developed by Gibson and Prowse (2002) that describes the fractionation of isotopes between water and freezing ice in a closed system:

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}_{\mathrm{L}\text{-}\mathrm{Pre}}=-f^{{\mathit{\alpha }}_{\mathrm{eff}}}\left(\mathrm{1000}\cdot f^{{\mathit{\alpha }}_{\mathrm{eff}}}-f\cdot {\mathit{\delta }}_{\mathrm{L}\text{-}\mathrm{BelowIce}}-\mathrm{1000}\cdot f\right),\end{array}$$

where δ_L-BelowIce is the isotopic composition of the water beneath the lake ice, α_eff is the effective fractionation factor between ice and water, defined as ${\mathit{\alpha }}_{\mathrm{eff}}=R_{\mathrm{Ice}}/R_{\mathrm{L}}$, and f is the fraction of unfrozen water remaining in the lake. α_eff is dependent on the thickness of the boundary layer between the forming ice and freezing water and the downwards freezing velocity of the ice. Since we did not have measurements of either of these parameters, we relied on previously estimated values of α_eff (Souchez et al., 1987; Bowser and Gat, 1995; Ferrick et al., 2002) and boundary layer thickness (Ferrick et al., 1998, 2002; Gibson and Prowse, 2002). Using this information, we estimated values of α_eff that produced δ_L-Pre values that closely match lake water isotopic compositions measured at the same lakes in August and September of 2018 (Fig. C1). Additional information about the determination of α_eff values is provided in Appendix C.

To estimate the fraction of unfrozen water remaining in lakes (f, Eq. 3), bathymetry was collected at Big Bear Lake (Fig. 2), a typical bowl-shaped thermokarst lake near the Trail Valley Creek meteorological station in June 2017 using a Garmin echoMAP CHIRP 42dv fish finder. Bathymetry data were used to determine a relationship between lake volume and lake depth. We fit a quadratic equation to the bathymetric data to estimate the fraction of lake volume relative to the fraction of lake depth. The best-fit quadratic equation (r²=0.9997) was

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\%}}_{\mathrm{LakeVolume}}=& {\left(-\mathrm{0.0115}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}}\right)}^{\mathrm{2}}\\ & +\mathrm{2.1508}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}}-\mathrm{0.4857},\end{array}\end{array}$$

where %_LakeVolume is the fraction of total lake volume and %_LakeDepth is the fraction of total lake depth. However, this fitted equation does not reach 100 %_LakeVolume at 100 %_LakeDepth or 0 %_LakeVolume at 0 %_LakeDepth, which is required to realistically represent the relationship between lake depth and lake volume. The equation was slightly adjusted to

$$\begin{array}{}\text{(5)}& {\mathit{\%}}_{\mathrm{LakeVolume}}={\left(-\mathrm{0.01}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}}\right)}^{\mathrm{2}}+\mathrm{2}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}}\end{array}$$

in order to satisfy these requirements, resulting in a mean offset of 1.7 % between the measured bathymetric data and the adjusted equation. Most lakes in this region have a bowl-shaped bathymetry because they were formed by thermokarst processes (Rampton, 1988; Burn and Kokelj, 2009), where subsidence caused by the thaw of ice-rich permafrost results in a water body which then expands outward radially in all directions. Bathymetric data for Big Bear Lake and a comparison of the equation between lake volume and depth are provided in Appendix D.

3.3 Quantifying lake and watershed properties

We quantified multiple lake and watershed properties to explore relations with the amount of lake water replaced by freshet. These properties included lake depth, lake volume, snow depth on the lake, ice thickness, lake area and watershed area. Lake depth, snow depth on the lake and ice thickness were measured at the same time as pre-snowmelt lake samples were collected. Lake volume was approximated by multiplying the product of lake depth and lake area by 0.7, a relation derived from the measured lake volume of Big Bear Lake. Watershed area was estimated by applying the D8 water routing algorithm (O'Callaghan and Mark, 1984) to the 2 m resolution ArcticDEM (PGC, 2018) using ArcGIS 10.7.1 (ESRI, 2019).

4 Results

Correcting for ice fractionation using Eq. (3) resulted in an increase in estimated δ_L-Pre values as expected, with the median shifting from −17.50 ‰ (−19.32 ‰ to −15.68 ‰ IQR, inter-quartile range) to −14.70 ‰ (−16.85 ‰ to −13.29 ‰ IQR) for δ¹⁸O (Fig. 5a, Table 2). The corrected pre-snowmelt isotopic compositions of lake water were distributed across a large range of the predicted local evaporation line (LEL), spanning from near the LMWL to near δ^* (Fig. 5a), reflecting the fact that the lake waters were variably influenced by evaporation. Corrected pre-snowmelt isotopic compositions of lake waters also tightly cluster along the LEL, indicating that the predicted LEL is well characterized.

Figure 5 Lake water and precipitation isotope data are displayed on δ¹⁸O–δ²H graphs. The local meteoric water line (LMWL: ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{7.1}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{10.0}$) is indicated by the solid line, while the LEL (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=\mathrm{5.2}\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-\mathrm{48.9}$) is indicated by the dashed line. δ_P represents the average value of precipitation in the region, based on 2018 sampling of end-of-winter snow and rainfall from April to September. δ_SSL is the point at which evaporation and inflow are equal ($E/I=\mathrm{1}$). (a) Uncorrected and corrected for ice fractionation pre-snowmelt lake water isotope data. (b) Corrected pre-snowmelt and post-snowmelt data. (c) The shift in δ_I from pre-snowmelt to post-snowmelt, as indicated by a circle for pre-snowmelt δ_I values and the end of the arrow for post-snowmelt δ_I values. All δ_I values are offset from the LMWL for visibility, as all δ_I values are constrained to the LMWL.

The change in the isotopic composition of lake water from pre-snowmelt to post-snowmelt was characterized by a small (∼1.5 ‰ in δ¹⁸O) shift towards δ_P, with median pre-snowmelt δ_L-Pre values of −14.70 ‰ (−16.85 ‰ to −13.29 ‰ IQR) and median δ_L-Post values of −16.15 ‰ (−17.86 ‰ to −15.45 ‰ IQR) for δ¹⁸O (Fig. 5b). The small change in the isotopic composition of lake water meant that most lakes retained an evaporated isotope signature post-snowmelt, overlapping with a substantial portion of δ_L-Pre and continuing to plot along the LEL (Fig. 5b). Post-snowmelt, about half of the lakes (9 of 17) also plotted above the LEL, indicating that the δ_I of these lakes was more similar to rainfall than snowfall (Fig. 5b). The shift in δ_I for lakes from pre-snowmelt to post-snowmelt shows a convergence of most δ_I values towards a value near δ_P and away from the isotopic composition of the end-of-winter snow (δ_Snow) or rainfall (δ_Rain) (Fig. 5c). The convergence of δ_I values towards δ_P and away from end-of-winter snow signifies that a non-snow source of water, with a more enriched isotopic composition than δ_Snow, was present in freshet.

Replacement of lake water by freshet ranged widely from 5.2 % to 52.8 %, with a median of 25.2 % (19.1 % to 32.4 % IQR, Fig. 6). A substantial proportion of this variation was explained by lake depth: deeper lakes had significantly less of their water replaced by freshet, with a reduction in lake water replacement of 13 % for each additional metre of lake depth (R²=0.53, p<0.001, Fig. 6, Table 3). Lake water replacement was not independently correlated with any other lake or watershed attribute, including watershed area, lake volume, snow depth on the lake ice, lake ice thickness, lake area, and the ratio of lake area to watershed area (Table 3).

Figure 6 The relationship between the amount of lake water replaced by freshet and lake depth. The distribution of lake water replacement by freshet is shown by the boxplot on the right-hand side of the plot. A linear regression is also displayed on the plot (R²=0.53, p<0.001). The shaded grey area represents the 95 % confidence interval of the linear regression.

Table 3 Results for a linear regression between total lake water replacement with multiple lake and watershed properties. The adjusted R² and p values are shown for each isotope. Linear regressions were performed using the “lm” function using R 4.0.2 (R Core Team, 2021).

5 Discussion

5.1 Influence of snowmelt bypass on the replacement of lake water by freshet

Characterization of the influence of snowmelt bypass on the replacement of pre-snowmelt lake water required accurate determination of the isotopic composition of lake water prior to freeze-up. Given that lake water isotope samples are unavailable from autumn 2017 and δ_L-Pre values were instead obtained from drilling through the lake ice before the lakes became ice-free, their isotopic compositions required correction for the isotope fractionation caused by ice formation. Our novel approach to correcting δ_L-Pre values for the fractionation caused by lake ice formation provides a reasonable estimate of δ_L prior to lake ice formation. While our correction of δ_L-Pre involves some uncertainty, such as having to estimate the relationship between lake depth and lake volume, corrected δ_L-Pre values closely align with the general distribution of water isotope measurements from August and September 2018 of the same lakes (Fig. C1). Corrected δ_L-Pre values are also situated near or above the LEL, reasonably indicating a more rainfall-sourced δ_I that would be present in lakes at the time of freeze-up during the previous autumn. Prior to correction, most δ_L-Pre values plotted below the LEL (Fig. 5a), indicating that lakes had the majority of their inflow sourced from snow, which would be unlikely at the time when lake ice began forming during the previous autumn. We considered using δ_L values from September 2018 instead of correcting for ice fractionation, but 2018 was a cooler and wetter year than 2017, meaning the lake-specific δ_L values in September 2018 likely differed somewhat from September 2017.

The presence of a somewhat uniformly thick layer of freshet beneath lake ice (Appendix E) likely explains the relationship between lake depth and the amount of lake water replaced by runoff, because the freshet layer represents a relatively smaller portion of lake volume at deeper lakes (Fig. 6). We calculated the freshet layer thickness to be an average of 0.28 m, ranging from 0.12 to 0.52 m with a standard deviation of 0.11 m (Table E1). Previous studies have measured the thickness of the snowmelt bypass layer at the onset of freshet inflow to be ∼ 30–200 cm (Henriksen and Wright, 1977; Bergmann and Welch, 1985). Since the mixed layer of pre-snowmelt lake water and freshet comprised a relatively larger volume in shallower lakes compared to deeper lakes (Fig. 7a), a larger portion of lake water was able to be replaced with freshet in shallower lakes than in deeper lakes. We hypothesized that because shallower lakes likely had colder lake-bed temperatures during freshet (Burn, 2005), more mixing between pre-snowmelt lake water and freshet inflow would occur due to the reduction in the water density gradient between the bottom of the lake and the top of the lake. However, the estimated thickness of the freshet layer was uniform across lakes, indicating that colder lake-bed temperatures may not have contributed to greater mixing at shallower lakes. To our knowledge, the relationship between lake depth and freshet retention has not been described in the previous literature, although there has been little possibility of observing this relationship because estimates of freshet recharge in more than one lake are scarce (Falcone, 2007; Brock et al., 2008). We expect that a similar relationship between lake depth and snowmelt bypass could be present in other open-drainage lakes that experience snowmelt bypass, since the relationship between increasing lake depth and greater snowmelt bypass is caused by the typical water temperature gradient which is present in ice-covered lakes at the onset of freshet. However, smaller lakes <1 ha, which are common in the Arctic (Pointner et al., 2019), likely do not experience as strong a snowmelt bypass effect because freshet is able to displace the pre-snowmelt lake water due to the small volume of the lake (Jansen et al., 2019; Cortés and MacIntyre, 2020). However, smaller lakes also typically have smaller under-ice chemical gradients that increase density stratification and limit freshet mixing and retention (Cortés and MacIntyre, 2020). Therefore, a relationship between lake depth and retention of freshet runoff likely does not exist at smaller Arctic lakes.

Figure 7 (a) A conceptual model showing the relative differences in snowmelt bypass between a shallow lake and a deep lake. Shallower lakes have a larger portion of their volume replaced by the runoff layer that flows beneath the ice, while a larger portion of water is isolated from mixing with runoff in deeper lakes. The estimated freshet layer thickness was somewhat uniform among all the lakes (Appendix E). (b) A conceptual model of how snowmelt bypass occurred over the course of the snowmelt period. Pre-snowmelt samples were taken from water beneath lake ice which was isotopically depleted in comparison to the lake ice. As time goes on, the source of freshet shifts from snow-sourced to active layer-sourced, while mixing increases beneath the lake ice simultaneously.

5.2 Sources of freshet

Following the freshet, the δ_I of lakes did not shift towards the isotopic composition of snow (δ_Snow), as one may expect, but instead shifted towards the average isotopic composition of precipitation (δ_P, Fig. 5c). Other than the 21.3 mm rainfall that fell during the 6-week period between the two isotope sampling dates, the only other potential source of water during this time period is water stored in the active layer, which likely mixed with snowmelt runoff as the soil thawed throughout the spring. The high infiltration capacity of the peat channels that convey runoff causes nearly all snowmelt to flow through subsurface routes, as was observed in the field and reported by Quinton and Marsh (1998). As water near the surface of the active layer is most likely to be comprised of rainfall from the previous year, we expect that much of the water stored in the top of the active layer would have been largely sourced from rainfall from the previous summer and autumn. In support of this inference, Tetzlaff et al. (2018) reported δ²H values between −140 ‰ and −160 ‰ from August to September 2014 in water samples taken at 10 cm soil depth at Siksik Creek, a watershed directly adjacent to the Trail Valley Creek camp (Fig. 2). This range of soil water δ²H values is more enriched relative to δ_P (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}=-\mathrm{160.1}\mathrm{‰}$, Fig. 5c), indicating that a mixture of snowmelt runoff with this active layer water could result in a water source similar in isotopic composition to δ_P.

Freshet flowing into lakes later during the snowmelt period likely had a more rainfall-sourced isotopic composition, replacing the more snow-sourced runoff that had entered the lake earlier during the freshet (Fig. 7b). A shift from snow-sourced water towards rainfall-sourced water during the course of the snowmelt period has been observed using water isotope measurements at Siksik Creek by Tetzlaff et al. (2018). Additionally, the mixing of freshet beneath lake ice increases with time as the temperature and density gradient lessens between the top and bottom of the lake water column (Cortés et al., 2017). Based on our results and these previous studies, we conclude that the ability of the active layer to contribute runoff to lakes appears to be maximized at the same time that vertical mixing in the lake is stronger, while snowmelt runoff flows into lakes at a time when little vertical mixing is occurring and is also likely replaced by later runoff (Fig. 7b). Such interplay between timing of snowmelt runoff, lake ice melt and hydrological behaviour of the active layer explains why the source of water to lakes is not solely snow-sourced and why incorporation of active layer runoff into lakes is more important than the volume of freshet delivered to lakes for the open-drainage lakes of this study.

5.3 Assumptions and improving the utility of water isotope data from ice-covered lakes

In our estimation of lake water replacement by freshet we had to make some assumptions (Appendix C) when estimating δ_L-Pre using Eq. (3). Future studies could sample lakes in the previous autumn before ice formation begins to avoid these assumptions, as there is minimal water flow in and out of Arctic lakes during the winter months due to frozen soils and ice cover on lakes (Woo, 1980). If lakes cannot be sampled in the previous autumn, another option would be to take a lake ice core and sample the isotopic composition at different points along the lake ice core. These isotope measurements could then be used to estimate the α_eff value used in Eq. (3), as has been done by Souchez et al. (1987) and Bowser and Gat (1995). This approach would avoid the assumptions we made in estimating α_eff outlined in Appendix C. However, with this approach, one still needs to know the volume of the lake ice relative to the volume of the remaining unfrozen water and must rely on lake bathymetry data or a depth-to-volume relationship, such as the one we derived using bathymetry data from Big Bear Lake (Eq. 4). Since we only have one survey of lake bathymetry, we do not know how well our depth-to-volume relationship describes other lakes in the region, and it could be that this relationship varies as lakes increase in surface area or in areas of different surficial geology, where the intensity of thermokarst processes may have varied during lake formation.

Since we do not have measurements of the lake temperature profile, we also assume that our lakes have the typical thermal structure of cryomictic (Yang et al., 2021) ice-covered lakes that leads to snowmelt bypass. We relied on measurements of water temperature at 1.25 m depth in Big Bear Lake (Fig. 2) and lake temperature profiles modelled using FLake-online (Appendix A) to conclude that our lakes were likely well-mixed at the time lakes were sampled post-freshet. Two scenarios were established in FLake-online, one representing a typical lake from the ones we sampled and another lake representing a small, deep lake where mixing would be less likely. Even though FLake-online does not simulate under-ice mixing, which typically occurs in Arctic lakes (Hille, 2015; Cortés et al., 2017; Kirillin et al., 2018), the “typical lake” became fully mixed 1 d after becoming ice-free, while the small, deep lake was simulated to become fully mixed 3 d after becoming ice-free. Even though snowmelt bypass is a common phenomenon in many types of ice-covered lakes around the world (Henriksen and Wright, 1977; Jeffries et al., 1979; Hendrey et al., 1980; Bergmann and Welch, 1985; Schiff and English, 1988; Edwards and McAndrews, 1989; Cortés et al., 2017), knowledge about how the thermal structure of our study lakes evolved over time and varied between lakes of different depths would have informed interpretation of our results. Such data could have helped better explain why shallow lakes retain more freshet runoff than deeper lakes and also could have helped confirm our hypothesis that water flowing into lakes later during the freshet mixes more readily with lake water.

5.4 Climate change and snowmelt bypass

Future changes in snowmelt bypass are dependent on whether climate change allows open-drainage lakes to persist or causes them to become closed-drainage, given that snowmelt bypass can only occur when lakes are at or above their outlet level. There are multiple consequences of Arctic warming that will influence lake water balance: changes in rainfall (Bintanja and Andry, 2017) and snowfall (Brown and Mote, 2009; Ernakovich et al., 2014), increases in active layer thickness (Walvoord and Kurylyk, 2016; Tananaev and Lotsari, 2022; Koch et al., 2022), the proliferation of deciduous shrubs (Loranty et al., 2018), and longer lake ice-free periods (Woolway et al., 2020). Whether the combination of these changes will result in an increase or decrease in runoff to lakes is currently unknown (Blöschl et al., 2019), making it difficult to predict whether open-drainage lakes will persist or whether some open-drainage lakes may shift to being closed-drainage under future climate. Due to this uncertainty, we discuss potential future changes in snowmelt bypass under two scenarios: (a) where lakes remain open-drainage in the future and (b) where some lakes become closed-drainage in the future.

If open-drainage lakes persist as such in the future, we suspect the freshet that is incorporated into lakes may shift towards being more rainfall-sourced. Projected rainfall increases (Bintanja and Andry, 2017) will likely leave the active layer in a more saturated state when the active layer freezes in autumn, potentially providing more water to lakes during the freshet. Other studies have already established a strong positive relationship between increased rainfall in the previous summer and a more efficient conversion of the snowpack into freshet (Bowling et al., 2003; Stuefer et al., 2017), indicating the presence of a strong connection between snowmelt runoff and water stored in the thawing active layer. Increasing shrub heights will advance snowmelt timing relative to lake ice breakup (Marsh et al., 2010; Wilcox et al., 2019; Grünberg et al., 2020), causing more snowmelt to flow into lakes at a time when there is limited below-ice vertical mixing. Combining earlier snowmelt timing with a more rain-saturated active layer could result in more freshet bypassing open-drainage lakes early during the freshet, with the active layer thawing deeper and shifting freshet more towards rainfall-sourced water by the time below-ice mixing begins.

Advancement of snowmelt timing and a shift from snow-sourced to rainfall-sourced freshet may have limnological implications for lakes. Cation and anion concentrations in snowmelt water tend to decrease over the course of the snowmelt period (Marsh and Pomeroy, 1999), while the pH of snowmelt runoff increases with time (Quinton and Pomeroy, 2006). Snowmelt also tends to have higher dissolved organic carbon (DOC) concentrations than summertime runoff (Finlay et al., 2006) and typically contributes the largest input of terrestrial organic mater to lakes in organic-rich landscapes (Townsend-Small et al., 2011; Olefeldt and Roulet, 2012). Balasubramaniam et al. (2015) observed that thermokarst lakes dominated by snow-sourced water tended to have lower pH, higher conductivity and higher DOC concentrations than lakes dominated by rain-sourced water. Based on these observations, as snowmelt occurs earlier in the Arctic, lakes may experience decreases in DOC and conductivity and increases in pH. Such changes to lake biogeochemistry caused by shifts in freshet runoff retention could affect the productivity and ultimately the climate feedbacks of these lakes. Future research could combine estimates of lake water replacement by freshet with water chemistry measurements to further our understanding of the impact of snowmelt bypass on lake chemistry and other limnological properties.

In a future where some open-drainage lakes become closed-drainage due to greater evaporation under a longer ice-free season, for example, we expect that such lakes will retain more freshet runoff than comparable open-drainage lakes, because closed-drainage lakes retain the additional freshet that is required to recharge the lake to its outlet level. Since snowmelt bypass cannot occur until a closed-drainage lake is recharged to its outlet level, we expect that freshet retention by closed-drainage lakes will not be as influenced by lake depth. Lakes with smaller ratios of watershed area to lake area (WA $/$ LA) are more prone to a more negative water balance (Marsh and Pomeroy, 1996; Gibson and Edwards, 2002; Turner et al., 2014; Arp et al., 2015). Therefore, we expect that lakes with relatively small WA $/$ LA ratios will be more prone to becoming closed-drainage, relying on freshet to recharge them to their outlet level and retain more freshet as a result. A corollary of this prediction is that other ice-covered lakes which currently lie below their outlet level at the onset of the freshet (i.e., closed-drainage lakes) likely retain more freshet than open-drainage lakes of a similar lake depth. A more saturated active layer at the onset of snowmelt combined with a greater amount of snowfall should increase the ability of freshet to recharge any closed-drainage lake.

An additional complication to predicting the future of snowmelt bypass in response to climate change is caused by the shifting of lake ice regimes from bedfast ice to floating ice. Lakes that freeze completely to the bed in winter (bedfast ice) melt from the surface downwards in spring and likely do not develop the thermal stratification necessary for snowmelt bypass to occur. Remote sensing studies have already observed many bedfast ice lakes shifting to floating ice regimes during the past few decades in response to climate change (Arp et al., 2012; Surdu et al., 2014; Engram et al., 2018). We are unaware of any lake mixing studies on bedfast ice lakes, making it difficult to hypothesize how shifting from bedfast to floating ice could affect the amount of freshet retained by a lake.

6 Conclusions

A portion of the large volume of freshet that flows into lakes every year can bypass ice-covered, open-drainage lakes due to limited mixing between lake water beneath the lake ice and freshet. By estimating the percentage of lake water replaced by freshet at 17 open-drainage lakes, we have been able to explore which lake and watershed attributes affect snowmelt bypass. Our data show that, as lake depth increases, the amount of lake water replaced by freshet decreases, likely because freshet is unable to mix with deeper lake water when lakes are ice-covered and the water column is stratified; however, we lack data demonstrating the extent of mixing in the lakes we studied. Additionally, the volume of freshet flowing into the lakes seems to have minimal impact on the amount of lake water replaced by freshet, given that the ratio of watershed area to lake area was not correlated with the percentage of lake water replaced by freshet. Estimation of the isotopic composition of source waters showed that the freshet remaining in lakes was not solely snow-sourced – rainwater left in the active layer from the previous autumn had mixed with snowmelt before entering lakes. Active layer-sourced water likely flows into lakes later in spring and at a time when freshet can more easily mix with pre-snowmelt lake water, replacing the earlier more snow-sourced freshet.

Models specialized for northern environments are rapidly improving their ability to represent the complicated processes present in permafrost regions, such as the effect of shrubs on snow accumulation, snowmelt and active layer thickness (Krogh and Pomeroy, 2019; Bui et al., 2020), lake ice formation and decay (MacKay et al., 2017) and the mixing processes that lead to snowmelt bypass (MacKay et al., 2017). Such models could be used to examine how freshet water sources may change in the future, which could have significant impacts on limnological properties including water chemistry (Finlay et al., 2006; Balasubramaniam et al., 2015). Additionally, current physically based lake models can represent vertical mixing beneath lake ice (MacKay et al., 2017) and could be used to further evaluate the influence of lake depth, lake ice regime or climate change on snowmelt bypass and resulting impacts on limnology.

Appendix A: Lake mixing status at ice-off: water temperature data and modelling

The application of Eq. (2) assumes that lakes are well-mixed at the time that δ_L-Post water samples were corrected. We investigated whether the lakes we sampled were likely to have been well-mixed at ice-off, because some lakes have been observed to be not well-mixed at ice-off (Vachon et al., 2019; Wiltse et al., 2020; Cortés and MacIntyre, 2020). We rely on water temperature data from Big Bear Lake and another lake near the Inuvik–Tuktoyaktuk Highway investigated by Hille (2015) and lake temperature modelling using FLake-online (Kirillin et al., 2011).

Water temperature measurements at Big Bear Lake and a lake nearby the Inuvik–Tuktoyaktuk Highway suggest that lakes in this region become well-mixed during the ice-off period. At Big Bear Lake, water temperature at 1.25 m depth reaches 4 ^∘C initially by 5 June 2018, followed by daily fluctuations between 2.3 and 4.1 ^∘C before continuing a warming trend again on 13 June 2018 and reaching 6.8 ^∘C on 15 June 2018 (Fig. A1). Water temperature between 0 and 4 m was measured by Hille (2015) at a 10 m-deep lake approximately 10 km from the Inuvik–Tuktoyaktuk Highway in 2009. The Hille (2015) measurements show uniform warming of the water column beneath the ice from 1 to 4 m, with water temperatures reaching ∼5 ^∘C by the time the lake became ice-free (Hille, 2015, Fig. 3.4). Based on these observations, we assume that these two lakes were well-mixed at the time they became ice-free.

Figure A1 Big Bear Lake water temperature at a depth of 1.25 m between 25 May and 19 June 2018.

Table A1 Input parameters and results for the two scenarios used to evaluate the mixing conditions of lakes during the ice-off period.

We also ran two model scenarios using FLake-online (Kirillin et al., 2011; http://flake.igb-berlin.de/model/run, last access: 31 July 2022) to gather further information about the mixing status of lakes in this region at the time they become ice-free. The first scenario represents a typical lake compared to the lakes we sampled in size and depth (Table A1). In the second scenario, lake depth is increased and lake area is decreased to represent the “worst-case scenario” for lake mixing after lakes became ice-free (Table A1). Water clarity was set to 2 m, based on an average Secchi depth measurement of 1.88 m based on measurements made by Vucic et al. (2020) at lakes along the Inuvik–Tuktoyaktuk Highway and the Dempster Highway south of Inuvik. Both scenarios were run with “perpetual year” meteorological forcing data, whereby meteorological data from 1 November 2015 to 31 October 2016 force the model for multiple years until a quasi-steady equilibrium is reached.

In both model scenarios, the mixing depth reached the maximum lake depth rapidly after ice-off, taking 1 d in the typical lake scenario and 3 d in the worst chance of mixing scenario. Notably, under-ice warming and mixing does not occur in FLake-online, although under-ice mixing has been observed in other Arctic lakes (Hille, 2015; Cortés et al., 2017; Kirillin et al., 2018). Given that both models indicate full mixing within a few days of lakes becoming ice-free, despite the absence of under-ice warming and mixing that was likely present at our study lakes based on temperature data from Big Bear Lake and Hille (2015), we believe that the lakes we sampled were very likely fully mixed when we sampled them on 15 June 2018.

Appendix B: Isotope framework and sensitivity analysis

The isotope framework (i.e., establishment of the predicted local evaporation line; LEL) used for this study was based on the coupled isotope tracer method developed by Yi et al. (2008), following other studies that have investigated lake water balances using water isotope tracers (Turner et al., 2014; Remmer et al., 2020; MacDonald et al., 2021). Below are the parameters and equations required to calculate δ^*, the terminal point on the LEL (Table B1). The equation for δ^*, which represents the isotopic composition of a lake at the point of desiccation, is as follows (Gonfiantini, 1986):

$$\begin{array}{}\text{(B1)}& {\mathit{\delta }}^{*}={\displaystyle \frac{h\cdot {\mathit{\delta }}_{\mathrm{As}}+{\mathit{\epsilon }}_k+\left({\mathit{\epsilon }}^{*}/{\mathit{\alpha }}^{*}\right)}{h-{\mathit{\epsilon }}_k-\left({\mathit{\epsilon }}^{*}/{\mathit{\alpha }}^{*}\right)}},\end{array}$$

where α^* is the fractionation factor between the liquid and vapour phases of water (Horita and Wesolowski, 1994), calculated for δ¹⁸O as

$$\begin{array}{}\text{(B2)}& {\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{V}}^{*}={\mathrm{2.718}}^{\left(-\mathrm{7.685}+\mathrm{6.7123}\cdot \frac{{\mathrm{10}}^{\mathrm{3}}}{T}-\mathrm{1.6664}\cdot \frac{{\mathrm{10}}^{\mathrm{6}}}{T^{\mathrm{2}}}+\mathrm{0.35041}\cdot \frac{{\mathrm{10}}^{\mathrm{9}}}{T^{\mathrm{3}}}\right)/\mathrm{1000}}\end{array}$$

and calculated for δ²H as

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{V}}^{*}\\ & ={\mathrm{2.718}}^{\left(\mathrm{1158.8}\cdot \frac{T^{\mathrm{3}}}{{\mathrm{10}}^{\mathrm{9}}}-\mathrm{1620.1}\cdot \frac{T^{\mathrm{2}}}{{\mathrm{10}}^{\mathrm{6}}}+\mathrm{794.84}\cdot \frac{T}{{\mathrm{10}}^{\mathrm{3}}}-\mathrm{161.04}+\mathrm{2.9992}\cdot \frac{{\mathrm{10}}^{\mathrm{9}}}{T^{\mathrm{3}}}\right)/\mathrm{1000}}.\end{array}\end{array}$$

The term ε^* is a separation term where

$$\begin{array}{}\text{(B4)}& {\mathit{\epsilon }}^{*}={\mathit{\alpha }}^{*}-\mathrm{1}.\end{array}$$

The term h represents the relative humidity of the air above the water and δ_As is the isotopic composition of atmospheric moisture during the open water season defined as

$$\begin{array}{}\text{(B5)}& {\mathit{\delta }}_{\mathrm{As}}={\displaystyle \frac{{\mathit{\delta }}_{\mathrm{Ps}}-{\mathit{\epsilon }}^{*}}{{\mathit{\alpha }}^{*}}},\end{array}$$

where δ_Ps is the average isotopic composition of precipitation (i.e., rainfall) during the open water season. The term ε_k is the kinetic fractionation separation term, defined as

$$\begin{array}{}\text{(B6)}& {\mathit{\epsilon }}_k=x*(\mathrm{1}-h),\end{array}$$

where x=0.0142 for δ¹⁸O and x=0.0125 for δ²H (Gonfiantini, 1986).

Gonfiantini (1986)Environment and Climate Change Canada (2019)Environment and Climate Change Canada (2019)Horita and Wesolowski (1994)Horita and Wesolowski (1994)Gonfiantini (1986)

Table B1 Variables used in the isotope framework and sources of their calculation.

Given that there is some variability in δ_Snow and δ_Rain values from samples collected, we conducted a sensitivity analysis to evaluate whether uncertainty in these values could affect δ_I and percentage lake water replacement sufficiently to change our interpretation of the results. To conduct the sensitivity analysis, we calculated the standard error of the mean (SEM) for δ_Snow and δ_Rain:

$$\begin{array}{}\text{(B7)}& \text{SEM}={\displaystyle \frac{\mathit{\sigma }}{\sqrt{n}}},\end{array}$$

where σ is the standard deviation of δ_Snow or δ_Rain values and n is the number of δ_Snow or δ_Rain samples. The SEM was added to δ_Snow and δ_Rain to calculate an “upper-bound” estimate and subtracted to calculate a “lower-bound” estimate. These upper- and lower-bound δ_Snow and δ_Rain values were then used to calculate upper- and lower-bound δ_P, δ_Ps and δ^* values (Table B2). These upper- and lower-bound values were then used to calculate upper- and lower-bound δ_I (Fig. B1a), and percentage lake water replaced values (Fig. B1b) were calculated. Overall, δ_I and percentage lake water replaced values change minimally between the upper- and lower-bound cases (Fig. B1a, b) and do not alter our interpretation of the results.

Table B2 Comparison of isotope framework parameters for upper- and lower-bound cases.

Figure B1 Comparison of upper- and lower-bound δ_I and percentage lake water replacement values against the base case.

Appendix C: Determination of α_eff values

In order to determine α_eff, two variables must be taken into account: the thickness of the ¹⁸O or ²H boundary layer across which heavy isotopes are diffusing from water into ice and the downward velocity of the freezing ice (Ferrick et al., 2002). If these two parameters are known, the fractionation factor can be estimated using a linear resistance model developed by Ferrick et al. (1998), which is similar in structure to the Craig and Gordon (1965) linear resistance model for evaporation. Ferrick et al. (1998) define the effective fractionation factor between ice and water as

$$\begin{array}{}\text{(C1)}& {\mathit{\alpha }}_{\mathrm{eff}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{S}}^{*}}{{\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{S}}^{*}+\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{S}}^{*}\right)\exp \left[\frac{-zv}{D_i}\right]}},\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{S}}^{*}$ is the equilibrium fractionation factor between ice and water (1.002909 for ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$, 1.02093 for ${}^{\mathrm{2}}\mathrm{H}{/}^{\mathrm{1}}\mathrm{H}$; Wang and Meijer, 2018), z is the ¹⁸O or ²H boundary layer thickness between the ice and water (mm), v is the velocity of ice growth (cm² d⁻¹), and D_i is the self-diffusion coefficient of ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ or ¹H²H¹⁶O at 0 ^∘C (cm² d⁻¹). As the boundary layer and the velocity of ice growth increase, α_eff moves from the value of ${\mathit{\alpha }}_{\mathrm{L}\text{-}\mathrm{S}}^{*}$ towards a value of 1 (i.e., no fractionation).

As we do not know the boundary layer thickness at the ice–water interface or the exact ice growth velocity for the lakes studied here, we relied on multiple other sources of information to estimate a probable upper and lower bound of α_eff. We took into account previous estimates of α_eff for ice–water fractionation (Souchez et al., 1987; Bowser and Gat, 1995; Ferrick et al., 2002) and boundary layer thickness from other studies of lakes (Ferrick et al., 1998, 2002; Gibson and Prowse, 2002). The boundary layer thickness between water and freezing ice in a lake was estimated to be between 1 and 6 mm by Ferrick et al. (1998); however, they revised this estimate with a more rigorous diffusion model to 1±0.3 mm for ¹⁸O and 0.4±0.2 mm for ²H (Ferrick et al., 2002). They also found that the boundary layer thickness remained mostly stable across different ice growth velocities, although the lowest ice growth velocity of ∼0.9 cm d⁻¹ had a boundary layer of ∼1.8 mm (Ferrick et al., 2002). The mean ¹⁸O α_eff values for two ice cores taken from the lake studied by Ferrick et al. (2002) were 1.0021 and 1.0020, with respective ice growth velocities of 3.7 and 4.1 cm d⁻¹. A 1 mm boundary layer was also estimated by Gibson and Prowse (2002) beneath river ice in northern Canada; however, they also suggest that the boundary layer thickness can reach up to 4 mm in quiescent lake water. Therefore, we assume a minimum boundary layer thickness of 1 mm and a maximum boundary layer thickness of 4 mm.

Figure C1 Comparison of pre-lake ice formation and August/September 2018 isotopic composition of lake waters. The α_eff values (δ²H α_eff=1.010, δ¹⁸O α_eff=1.0015) match closely with the August and September 2018 isotopic composition of lake waters. This suggests that these α_eff estimates are appropriate to use for estimating the pre-ice formation isotopic composition of lake waters.

We estimated the minimum possible freezing velocity of our study lakes using the initial date of ice formation and the ice thickness we measured in spring. Based on Sentinel imagery (Sentinel Playground), all studied lakes became ice-covered by 16 October 2017. Ice thickness was measured at Big Bear Lake and Little Bear Lake (near Trail Valley Creek camp, Fig. 1) on 21 March 2018, and when ice thickness was remeasured again in late April, it had not become any thicker. Therefore, assuming ice growth began on 16 October 2017 and ceased on 21 March 2018, the ice growth velocities for our study lakes range from an average of 0.50–0.84 cm d⁻¹ (0.78–1.32 m ice thickness). This only provides a lower-bound estimate for ice growth velocity, as ice growth likely stopped earlier than 21 March 2018 and was more rapid during initial ice formation.

We further constrained our estimate of α_eff by assuming that α_eff values that result in lake water replacement estimates of >100 % or <0 % were not correct. Using all these sources of information, we calculated an upper bound of α_eff values based on the minimum possible ice freezing velocity (²H α_eff=1.0199, ¹⁸O α_eff=1.00286) and a lower bound of α_eff values which still generate lake water replacement estimates that are >0 % and <100 % (²H α_eff=1.010, ¹⁸O α_eff=1.0015).

Assuming a 2 mm boundary layer, which is within the range of our boundary layer thickness estimates, α_eff values of 1.0015 for ¹⁸O and 1.010 for ²H correspond to ice growth rates of 3.62 cm d⁻¹ for ¹⁸O and 3.34 cm d⁻¹ for ²H. Similar ice growth rates have been observed in Arctic lakes (Woo, 1980); greater ice growth rates were also estimated for a lake in a warmer climate by Ferrick et al. (2002). These α_eff values also compare well with other measured α_eff values in lakes: a range of α_eff= 1.013–1.015 for ²H was found by Souchez et al. (1987) for a 4.4 cm-thick lake ice cover; α_eff has been found to range from 1.0005 to 1.0027 for ¹⁸O in a single 50 cm ice core (Bowser and Gat, 1995). The δ_L-Pre values calculated using α_eff=1.0015 for ¹⁸O and 1.010 for ²H also closely match the distribution of δ_L values from August and September of 2018, giving an indication that these α_eff values generate realistic isotopic compositions of lake water (Fig. C1). Therefore, we chose α_eff=1.0015 for ¹⁸O and α_eff=1.010 for ²H as our α_eff values, as they correspond well to other estimates α_eff, are within a range of probable ice-growth rates and lake water replacement by freshet and generate pre-ice formation isotopic compositions that closely match the following summer's isotopic composition of lake waters.

Appendix D: Bathymetric data and volume–depth relationship

The modelled relationship between lake depth and volume is

$$\begin{array}{}\text{(D1)}& {\mathit{\%}}_{\mathrm{LakeVolume}}={\left(-\mathrm{0.01}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}}\right)}^{\mathrm{2}}+\mathrm{2}\cdot {\mathit{\%}}_{\mathrm{LakeDepth}},\end{array}$$

where %_LakeVolume is the cumulative lake volume as a percent of the total and %_LakeDepth is the cumulative lake depth as a percent of the total (Table D1).

Table D1 Big Bear Lake bathymetric data and fit between the modelled relationship between lake volume and lake depth as a percentage of total lake volume and depth.

Appendix E: Freshet layer thickness

We computed the thickness of the freshet layer using the relationship between lake depth and lake volume, the percentage of lake volume replaced by freshet, and lake depth measurements made at each sample lake. The freshet layer thickness was calculated by rearranging Eq. (5):

$$\begin{array}{}\text{(E1)}& \begin{array}{rl}& \text{freshet layer thickness (m)}\\ & ={\displaystyle \frac{\text{lake depth}}{-\mathrm{10}\cdot \sqrt{\mathrm{100}-\mathit{\%}\text{lake water replaced}}-\mathrm{10}}}.\end{array}\end{array}$$

This calculation represents the thickness the freshet layer would be on 15 June 2018 if it had not mixed with pre-snowmelt lake water. In reality, the unmixed freshet layer during the height of freshet would likely be thicker due to rises in lake level caused by freshet that are not accounted for in this equation. Layer thicknesses averaged 0.28 m, ranging from 0.12 to 0.52 m with a standard deviation of 0.11 m (Table E1).

Table E1 Calculated layer thickness for each lake using Eq. (E1).