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 what 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 limnological properties of the lake. To explore lake 5 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 source 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 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 10 freshet layer was not proportional to lake depth, which isolated a larger portion of pre-snowmelt lake water from mixing at deeper lakes. We expect a similar relationship between increasing lake depth and greater snowmelt bypass could be present at all ice-covered open-drainage lakes that are poorly 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 15 ice melt, snowmelt bypass may become more prevalent with the water remaining in thermokarst

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 cause less inflow and more water to evaporate from thermokarst lakes. Contrarily, increasing precipitation 30 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). 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.
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 flow of runoff into lakes (freshet) generally results in the highest lake levels of the year (Woo, 1980;Roulet and Woo, 1988;Hardy, 1996; Pohl 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 Delta is in grey. In the inset, key locations near the Trail Valley Creek field station are shown. and many lakes have defined channelized inflows from their watersheds in the form of small streams or ice-wedge polygon troughs.
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 80 lay between mineral earth hummocks (Quinton and Marsh, 1998). Lake watersheds contain tall shrub (>1 m), low-shrub (~0.5 m), and shrub-free landcover types comprising lichen, moss, and tussocks 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). 85 The 2018   The mean air temperature at Trail Valley Creek during the sampling period from April 26 -June 15, 2018 was 0.4 • C, which was cooler than the average of 1.7 • C during 1999-2019 ( Figure 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 May 8, which was similar to the average timing of the first above 0 • C day during 1999-2019 ( Figure 3). 95 3 Methods 3.1 Lake water and precipitation sampling for isotope analysis Lake water samples for isotope analysis were first collected from the 17 study lakes while they were ice covered (April 26 -May 1, 2018) and again soon after lakes became ice-free (June 15, 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 100 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 isotope composition of precipitation (δ P ) in the study 105 region, which are useful references for the interpretation of lake water isotope 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 where R represents the ratio of 18 O/ 16 O or 2 H/ 1 H. Every fifth sample was analyzed a second time to determine the analytical uncertainties, which were ±0.1‰ for δ 18 O and ±0.6‰ for δ 2 H, calculated as two standard deviations away from the difference 120 between the duplicate samples. All isotope data from lakes are presented in Table 2.
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: where δ L-Pre is the lake water isotope composition before snowmelt begins, δ L-Post is the lake water isotope composition 125 after snowmelt is complete, and δ I-Post is the isotope 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 (Figure 2), previous water temperature profiles made at lakes within 10 km of the Inuvik- Max 4.14 1. 32 -13.98 -120.75 -12.02 -110.51 -17.44 -133.22 -13.59 -117.58 -17.7 -135.07 52.8 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 2018-06-15 (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 (June 15 -October 15) collected at the Trail Valley Creek meteorological station located 45 km NNE of Inuvik ( Figure 2). Data from 2017 were used because it was the most recent period where lakes were exposed to meteorological conditions for an entire open-water season. The coupled isotope 135 tracer approach assumes all lakes under the same meteorological conditions will evolve towards the same isotope composition (δ * , the isotope 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 (Figure 4). We calculated δ I for pre-snowmelt and post-snowmelt lake water isotope compositions to identify whether the isotope composition of the source water changed after freshet. The percentage of lake 140 water replaced was calculated using both δ 18 O and δ 2 H using Equation 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 -3.7%.
A sensitivity analysis was performed to evaluate uncertainty in the δ I values and subsequent % lake water replaced by runoff.
Confidence in the interpretation of δ I values with respect to rainfall-or snowmelt-sourced waters depends on an accurate 145 estimate of δ P , which was determined using the average δ Rain and δ Snow values. Also, the calculation of % lake water replaced by runoff is sensitive to changes in δ I , which is sensitive to the average δ Rain value because this paramter is used to determine δ As (Equation B5), which is subsequently used to determine δ * (Equation 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 150 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 % 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 % lake water replacement values were minimal (Table B2, Figure   B1). Details of the equations and variables used in the isotope framework and the sensitivity analysis are given in Appendix B.

155
As ice forms and preferentially incorporates water containing the heavy isotopes 18 O or 2 H, the lake water beneath the ice becomes increasingly depleted in 18 O and 2 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: where δ L-BelowIce is the isotope composition of the water beneath the lake ice, α eff is the effective fractionation factor between ice and water, defined as α eff = R Ice /R 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 165 et al., 1987;Bowser and Gat, 1995;Ferrick et al., 2002) and boundary layer thickness (Ferrick et al., 1998(Ferrick et al., , 2002Gibson and Prowse, 2002). Using this information, we estimated values of α eff that produced δ L-Pre values that closely match lake water isotope compositions measured at the same lakes in August and September of 2018 ( Figure 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 , Equation 3), bathymetry was collected at Big Bear Lake 170 ( Figure 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 2 = 0.9997) was: 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: in order to satisfy these requirements, resulting in a mean offset of 1.7% between the measured bathymetric data and the ad-180 justed 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 waterbody 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.

Quantifying lake and watershed properties 185
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.

Results
Correcting for ice fractionation using Equation 3 Table 2). The corrected pre-snowmelt lake water isotope compositions were distributed across a large range of the 195 predicted LEL, spanning from near the LMWL to near δ * (Figure 5a), reflecting that the lake waters were variably influenced by evaporation. Corrected pre-snowmelt lake water isotope compositions also tightly cluster along the LEL, indicating that the predicted LEL is well characterized.
The change in lake water isotope composition from pre-snowmelt to post-snowmelt was characterized by a small (~1.5‰  (Figure 5b). The small change in lake water isotope composition 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 (Figure 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 (Figure 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 isotope  away from end-of-winter snow signify that a non-snow source of water, with a more enriched isotope composition than δ Snow , was present in freshet.
Replacement of lake water by freshet ranged widely from 5.2 -52.8%, with a median of 25.2% (19.1% to 32.4% IQR, Figure 6). A substantial proportion of this variation was explained by lake depth: deeper lakes had significantly less of their 210 water replaced by freshet, with a reduction in lake water replacement of 13% for each additional metre of lake depth (R 2 = 0.53, p < 0.001, Figure 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). Table 3. Results for a linear regression between total lake water replacement with multiple lake and watershed properties. The adjusted R 2 and p-value are shown for each isotope. Linear regressions were performed using the 'lm' function using R 4.0.2 (R Core Team, 2021) . % Lake Volume Replaced Characterization of the influence of snowmelt bypass on the replacement of pre-snowmelt lake water required accurate determination of lake water isotope compositions 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 isotope compositions required correction for the isotope fractionation caused by ice formation. Our novel approach to 220 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 ( Figure 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 225 0 20 40 1 2 3 4 Lake Depth (m) Lake Volume Replaced by Freshet (%) 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 side of the plot. A linear regression is also displayed on the plot (R 2 = 0.53, p < 0.001). The shaded grey area represents the 95% confidence interval of the linear regression.
autumn. Prior to correction, most δ L-Pre values plotted below the LEL (Figure 5a), indicating 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 230 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 ( Figure 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 235 compared to deeper lakes (Figure 7a), a larger portion of lake water was able to be replaced with freshet in shallower lakes than in deeper lakes. Shallower lakes likely had colder lakebed temperatures during freshet (Burn, 2005), which allowed more mixing between pre-snowmelt lake water and freshet inflow due to the reduction in water density gradient between the bottom of the lake and the top of the lake. To our knowledge, the relationship between lake depth and freshet retention has not been described in previous literature, although there has been little possibility to observe this relationship because estimates of 240 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 (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.
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 lakes <1 ha, which are common in the Arctic (Pointner et al., 2019), likely do not experience snowmelt bypass because freshet is able to displace the pre-snowmelt 245 lake water due to the small volume of the lake (Jansen et al., 2019; Cortés and MacIntyre, 2020). Therefore, a relationship between lake depth and retention of freshet runoff likely does not exist at smaller arctic lakes.

Sources of freshet
Following the freshet, the δ I of lakes did not shift towards the isotope composition of snow (δ Snow ) as one may expect, but instead shifted towards the average isotope composition of precipitation (δ P , Figure 5c). Other than the 21.3 mm rainfall that 250 fell during the six-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 255 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 δ 2 H values between -140‰ and -160‰ from August to September of 2014 in water samples taken at 10 cm soil depth at Siksik Creek, a watershed directly adjacent to the Trail Valley Creek camp (Figure 2). This range of soil water δ 2 H values is more enriched relative to δ P (δ 2 H = -160.1 ‰, Figure 5c), indicating that a mixture of snowmelt runoff with this active layer water could result in a water source similar in isotope composition to δ P .

260
Freshet flowing into lakes later during the snowmelt period likely had a more rainfall-sourced isotope composition, replacing the more snow-sourced runoff that had entered the lake earlier during the freshet (Figure 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 265 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 (Figure 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 snowsourced, and that incorporation of active layer runoff into lakes is more important than the volume of freshet delivered to lakes, 270 for the open-drainage lakes of this study.

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 Equation 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 275 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 isotope composition at different points along the lake ice core. These isotope measurements could then be used to estimate the α eff value used in Equation 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 280 bathymetry data or a depth-to-volume relationship, such as the one we derived using bathymetry data from Big Bear Lake (Equation 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 our lakes have the typical thermal 285 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 (Figure 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 290 typically occurs in arctic lakes (Hille, 2015;Cortés et al., 2017;Kirillin et al., 2018), the 'typical lake' became fully mixed 1 day after becoming ice-free while the small, deep lake was simulated to become fully mixed 3 days 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 295 varied between lakes of different depth, 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.

Climate change and snowmelt bypass
Future changes in snowmelt bypass are dependent on whether climate change allows open-drainage lakes to persist, or causes 300 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 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 315 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 impli-320 cations 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). 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 325 snowmelt occurs earlier in the Arctic, lakes may experience decreases in DOC, and conductivity, and increases in pH. 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 such lakes will retain more freshet runoff than comparable open-drainage lakes, because 330 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 lakes with relatively small WA/LA ratios will be more prone to becoming closed-drainage, relying on 335 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 340 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 about how shifting from bedfast to floating ice could 345 affect the amount of freshet retained by a lake.

Conclusions
The large volume of freshet that flows into lakes every year is likely to 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 350 show that as lake depth increases the amount of lake water replaced by freshet decreases, because freshet is unable to mix with deeper lake water when lakes are ice-covered and the water column is stratified. Additionally, the volume of freshet flowing into the lakes seems to have no observable impact on the amount of lake water replaced by freshet. Estimation of the isotope 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 355 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 360 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

370
The application of Equation 2 assumes that lakes are well-mixed at the time that δ L-Post water samples were corrected. We investigated if 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).

375
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 05-06-2018, followed by daily fluctuations between 2.3 -4.1°C, before continuing a warming trend again on 13-06-2018 and reaching 6.8°C on 15-06-2018 ( Figure A1). Water temperature between 0 and 4 m was measured by Hille 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, Figure 3.4). Based on these observations, we assume these two lakes were well-mixed at the time they became ice-free. We also ran two model scenarios using FLake-online (Kirillin et al., 2011, http://flake.igb-berlin.de/model/run) to gather further information about the mixing status of lakes in this region at the time they become ice-free. The first scenario represents 385 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 Dempster Highway south of Inuvik. Both scenarios were run with 'perpetual year' meteorological forcing data, whereby meteorological data from 01-11-2015 to 31-10-2016 forces 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 one day in the typical lake scenario and three days 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 395 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 2018-06-15. Days from Ice-free to Fully-mixed 1 3 Table A1. Input parameters and results for the two scenarios used to evaluate the mixing conditions of lakes during the ice-off period.
Appendix B: Isotope framework and sensitivity analysis  (Horita and Wesolowski, 1994) 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 400 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. The equation for δ * , which represents the isotope composition of a lake at the point of desiccation, is as follows (Gonfiantini, 1986): where α * is the fractionation factor between the liquid and vapour phase of water (Horita and Wesolowski, 1994), calculated 405 for δ 18 O as: The term h represents the relative humidity of the air above the water and δ As is the isotope composition of atmospheric moisture during the open water season defined as: where δ Ps is the average isotope composition of precipitation (i.e., rainfall) during the open water season. The term ε k is the 415 kinetic fractionation separation term, defined as where x = 0.0142 for δ 18 O and x = 0.0125 for δ 2 H (Gonfiantini, 1986).
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 % lake water replacement sufficently to change our interpreta-420 tion of the results. To conduct the sensitivity analysis, we calculated the standard error of the mean (SEM) for δ Snow and δ Rain : 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 425 lower bound δ Snow and δ Rain values were then used to calculate upper and lower bound δ P , δ Ps and δ * values (Table B2)  In order to determine α eff , two variables must be taken into account: the thickness of the 18 O or 2 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: where α * L−S is the equilibrium fractionation factor between ice and water (1.002909 for 18 O/ 16 O, 1.02093 for 2 H/ 1 H (Wang and Meijer, 2018)), z is the 18 O or 2 H boundary layer thickness between the ice and water (mm), v is the velocity of ice growth (cm 2 day -1 ), and D i is the self-diffusion coefficient of 1 H 2 18 O or 1 H 2 H 16 O at 0°C (cm 2 day -1 ). As the boundary layer and the velocity of ice growth increase, α eff moves from the value of α * L−S towards a value of 1 (i.e., no fractionation).

440
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(Ferrick et al., , 2002Gibson and Prowse, 2002). The boundary layer thickness between water and freezing ice in a lake was estimated to be between 1 mm and 6 mm by Ferrick  (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 day -1 had a boundary layer of~1.8 mm (Ferrick et al., 2002). The mean 18 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 day -1 . 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.  (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 2 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 18 O in a single 50 cm ice core (Bowser and Gat, 1995). The δ L-Pre values calculated using α eff = 470 1.0015 for 18 O and 1.010 for 2 H also closely match the distribution of δ L values from August and September of 2018, giving an indication that these α eff values generate realistic lake water isotope compositions ( Figure C1). Therefore, we chose α eff = 1.0015 for 18 O and α eff = 1.010 for 2 H as our α eff values, as they correspond well with other estimates α eff , are within a range of probable ice-growth rates and lake water replacement by freshet and generate pre-ice formation isotope compositions that closely match the following summer's lake water isotope compositions.

475
Appendix D: Bathymetric data and volume -depth relationship Table D1. Big Bear Lake bathymetric data and fit between modelled relationship between lake volume and lake depth as a percentage of total lake volume and depth. The modelled relationship between lake depth and volume is: where % LakeVolume is the cumulative lake volume as a percent of total and % LakeDepth is the cumulative lake depth as a percent of total (Table D1). Appendix E: Freshet layer thickness We computed the thickness of the freshet layer using the relationship between lake depth and lake volume, the % of lake volume replaced by freshet, and lake depth measurements made at each sample lake. The freshet layer thickness was calculated by rearranging Equation 5: f reshet layer thickness (m) = lake depth −10 * 100 − %lake water replaced − 10 (E1)

485
This calculation represents the thickness the freshet layer would be on 2018-06-15 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 m to 0.52 m with a standard deviation of 0.11 m (Table E1).