Modeling δD-δ18O Steady-State of Well-Sealed Perennially Ice-Covered Lakes and Their Recharge Source: Examples From Lake Untersee and Lake Vostok, Antarctica

Perennially ice-covered lakes that are tightly sealed from the atmosphere represent a unique group of polar lakes. In these lakes, the δD-δ18O evolution of the water column and steady-state conditions are controlled by rates of recharge and freezing at the bottom of the ice cover. We developed a recursive model (FREEZCH9) that takes into account the changing salinity in the water column as a result of freezing and mixes the recharge water to the residual water in well-sealed perennially ice-covered lakes. Our model is tested against datasets from Lake Vostok and is used to assess the δD-δ18O mass balance of Lake Untersee and evaluate if the lake is in isotopic steady-state. Our FREEZCH9 simulations fit well with the predicted δD-δ18O values of Lake Vostok’s upper water column and the overlying accreted ice. Simulations with FREEZCH9 also suggests that Lake Untersee is in isotopic steady-state and that its two input sources (i.e., subaqueous terminus melting of the Anuchin Glacier and subglacial meltwater) have similar δD-δ18O composition. Our modeling demonstrates that Lake Untersee most likely did not receive additional input from surface streams during the last 300–500 years. FREEZCH9 may be also used to determine if any groundwater systems of the McMurdo Dry Valleys are fully or partially recharged by subglacial lakes.


INTRODUCTION
In Antarctica, several hundred lakes have been identified and classified according to their ice cover regime ( Figure 1A; Morgan et al., 2007;Vincent et al., 2008). Seasonally icecovered lakes are found along the warmer Antarctic Peninsula and other coastal regions. These lakes lose their ice cover in austral summer and are typically recharged by a combination of local precipitation and glacial meltwater (e.g., Gibson, 1999;Lyons et al., 2013). Perennially ice-covered lakes are found in polar desert regions and exist as two main types: those that develop moats in summer and are fed by glacial melt streams, and those that are tightly sealed from the atmosphere and solely fed by subaqueous  Touzeau et al., 2016). (B) Stable water isotope (δD-δ 18 O) ratios of perennially ice-covered lakes in the Antarctic: Princess Astrid Coast (Lake Untersee), Southern Victoria land region (Lake Miers, Lake Vanda, Lake Bonney, Lake Fryxell, Lake Vida) (Matsubaya et al., 1979;Miller and Aiken, 1996;Dugan et al., 2015) and Soya Coast (Skallen Oike) (Matsubaya et al., 1979) glacial melt and groundwater (e.g., Priscu, 1998;Faucher et al., 2019). Further, many subglacial lakes are found beneath the Antarctic ice sheet and are recharged by a combination of subglacial and englacial waters (e.g., Wingham et al., 2006;Priscu and Foreman, 2009;Wright and Siegert, 2012).
The source of water recharging the various types of Antarctic lakes and its evolution can influence key physicochemical conditions, the composition of microbial communities and biogeochemical cycling (e.g., Lyons and Finlay, 2008;Wand et al., 2006). The analysis of the δD-δ 18 O composition of polar lakes has proven to be a reliable tool to determine the source of water, the process that evolved the lake chemistry, and water balance parameters, such as the amount of evaporation (Jeffries et al., 1984;Gibson, 1999;Gooseff et al., 2006;Horita, 2009). In Antarctica, the δ 18 O values of lake waters range from c. 0 to −60 ( Figure 1B), a range that reflects the main sources of water recharging the lakes. Lakes with δ 18 O values near 0 are a heritage of isostatic rebound and their formation is attributed to the freezing of trapped seawater following decrease in sea level (Bird et al., 1991). Lakes with δ 18 O values near −15 to −10 have been attributed to mixing of seawater and local coastal precipitation (Matsubaya et al., 1979;Richter and Strauch, 1983;Bird et al., 1991). Finally, lakes with δ 18 O values <−25 are recharged mainly from glacial meltwater, with those with very low δ 18 O composition likely receiving input from late Pleistocene meltwater (Matsubaya et al., 1979;Richter and Strauch, 1983;Bird et al., 1991). The ice-cover regime of the lake (i.e., seasonally or perennially ice-covered) determines whether evaporation or freezing then drives the evolution of the lake isotope chemistry (Miller and Aiken, 1996). Both seasonally ice-covered lakes and perennially icecovered lakes that develop summer moats have lake waters with δD-δ 18 O composition that are typically distributed below the Antarctic Meteoric Water Line (AMWL) (Figure 1B). The isotopic behavior of these lakes is strongly influenced by evaporation that enriches the residual water and distributes it along an evaporative line with lower slope than the AMWL. The δ 18 O composition of evaporative lakes is determined by atmospheric parameters such as humidity and isotopic composition of the water vapor, and numerical models of the evolving salinity in the water column have been developed to estimate a lake's δD-δ 18 O steady-state (i.e., Gibson, 1999;Horita, 2009).
A few Antarctic lakes have δD-δ 18 O values distributed above the AMWL ( Figure 1B): two of these are Lake Untersee, a wellsealed perennially ice-covered lake, and Lake Vostok, a subglacial lake beneath the East Antarctic Ice Sheet. Consequently, the water columns of these lakes do not interact with the atmosphere (i.e., no fractionation due to evaporation) and the δD-δ 18 O evolution of each water column is controlled by recharge and the freezing at the bottom of the ice cover. In the case of Lake Vostok, the recharge water source, the freezing process and its isotopic steady-state are well established (e.g., Royston-Bishop et al., 2004;Souchez et al., 2004;Ekaykin et al., 2010). However, the isotopic mass-balance of Lake Untersee remains unknown. In addition, unlike evaporative lakes (i.e., Gibson, 1999;Horita, 2009), no studies have yet developed a model capable of predicting the δD-δ 18 O steady-state of the water column that includes the evolution of the ice-water isotope fractionation factor as freezing changes the salinity of the residual water and recursive recharge that mixes with the residual lake water.
The objective of this study is to develop a model of δDδ 18 O evolution of well-sealed perennial ice-covered lakes that accounts for changing chemistry in the residual water and annual recharge that mixes with residual water. This objective was achieved by using FREEZCH5, an isotope-augmented version of the FREZCHEM low temperature geochemistry model (Marion et al., 2010;Fisher et al., 2020) and making the model recursive. The model is tested against datasets from Lake Vostok and is then used to assess the δD-δ 18 O mass balance of Lake Untersee and evaluate if the lake is in isotopic steady-state. The results are presented in terms of residence times of water with different mixing scenarios, which allow to apply the findings to other well-sealed perennial ice-covered lakes or subglacial lakes.

Background on Water-Ice Fractionation During Freezing
The freezing of water under equilibrium condition follows a Rayleigh-type fractionation:  (Suzuoki and Kimura, 1973). During freezing, the forming ice is 18 O enriched but the residual water is progressively depleted and evolves above the Global Meteoric Water Line (GMWL) because the freezing line has a lower slope (Horita, 2009;Lacelle, 2011). Freezing of water also imparts an ionic segregation that increases the solute content in the residual water (Terwilliger and Dizio, 1970;Killawee et al., 1998). Because freezing reduces the water activity, a correction to the D/H and 18 O/ 16 O water-ice fractionation factors should be applied, especially at higher salinities where solutes may bind water molecules in hydration spheres around the cations and following the precipitation of hydrated salts (i.e., NaCl·2H 2 O, MgCl 2 ·12H 2 O, CaCl 2 ·6H 2 O). The theory of changing ice-water isotope fractionation under changing water salinities was described in Gat (1972, 1975) and Stewart and Friedman (1975) and they showed that the isotopic salt effect is linearly correlated for different salt solutions at varying concentrations.
In addition to evolving isotope salt effect, Jouzel and Souchez (1982) observed that the initial isotope composition of water also affects the isotope composition of ice during freezing. Consequently, Souchez and Jouzel (1984) suggested that the influence of the initial isotope composition of water on the water-ice fractionation factor (α δO ) could be estimated from: where δ O is the initial isotope composition of water, α i−w is the equilibrium fractionation factor. The equation shows that α δO is reduced with decreasing initial isotope composition of water; for most waters, this correction can be ignored as it only exceeds the analytical error when δ 18 O < −40 . Finally, the water-ice isotope fractionation during freezing is also dependent on the kinetics of the reaction, which includes the freezing rate, boundary-layer thickness and isotope diffusion at the ice-water interface. The kinetic water-ice isotope fractionation factor (α FR ) can be approximated by a modified version of Burton et al.'s (1953) equation for solute segregation: where h is the thickness of the boundary layer (cm), v is the freezing rate (cm s −1 ), D * is the diffusion coefficient in water. At 0C, D * values for H 2 18 O and HDO are 1.1 10 −5 cm 2 s −1 and 1.33 10 −5 cm 2 s −1 , respectively (Wang et al., 1953).
When the freezing rate is slower than the diffusion coefficient, α FR approximates the equilibrium fractionation factor, but under increasing freezing rate, α FR is reduced towards unity (i.e., Lacelle, 2011).

Description of FREEZCH5 (Isotope-Augmented FREZCHEM)
The core of the FREEZCH5 model is FREZCHEMv15, an equilibrium chemical thermodynamic model over the temperature range of −73 to 25 • C (Marion and Kargel, 2008;Fisher et al., 2020). FREZCHEM simulates the freezing of water by decreasing the temperature in fixed steps and determines the presence of residual water if the water activity calculated from the Pitzer equations is less than the equilibrium constant for water-ice. At each decreasing temperature step ("i"), the model provides: (1) the amount of residual water (in g; WAT LIQ ); (2) the amount of ice formed (in g; ICE ); (3) the molalities of the major cations (Na + , Ca 2+ , Mg 2+ , K + ) and anions (Cl − , NO 3 − , SO 4 2− ); (4) the amount of hydrated salts that precipitates (in g; WATNa, WATMg, and WATCa); and (5) the total ionic strength of the solution (proportional to conductivity and salinity of the residual water). FREZCHEM completely segregates the solutes in the residual water and none are incorporated as inclusions in the ice lattice or along grain boundaries and likely over-estimates the salinity in the residual water (i.e., Killawee et al., 1998;Petrenko and Whitworth, 1999). However, properly accounting for ionic segregation in equilibrium modeling is challenging as it was found that the concentration factor (or the segregation coefficient) is dependent on initial water chemistry and freezing rates (Santibáñez et al., 2019).
The evolution of δ 18 O of the residual water (δ 18 Orw) and for the incremental ice that forms at each decreasing temperature step was added to FREZCHEM in a routine called FREEZCH5, where the initial δ 18 O composition of the water and the α i−w is specified by the user: where δ 18 O rw is given by: where α 18O−i−w is the ice-water fractionation factor (corrected for initial isotopic composition using equation 1), Note that δ 18 O rw is the isotopic composition for the total water in the solution (free water and the water attached in the hydration spheres around the positive ions). The various hydrated salts do not have reported waterprecipitate fractionation factors and we assume that they are unity or as defined by Gat (1972, 1975) for the cationshydration spheres within the liquid solution. Since FREZCHEM is an equilibrium model, the effect of freezing rates on waterice fractionation factors are not included, but these should have little effect on the freezing of perennially ice-covered lakes where freezing takes place near 0 • C at the ice-water boundary.
With δ 18 O rw , we can use the Gat (1972, 1975) isotope salt effects correction to determine the δ 18 O offset between the δ 18 O total water (free and hydration sphere waters) and δ 18 O free water only (i.e., water activity), a correction that depends on the molality of each cation (Na + , Ca 2+ , Mg 2+ , K + ): (6) where δ 18 Ocations are functions of the remaining liquid water molality, which are calculated at each temperature step. The

Description of FREEZCH9 (a Recursive FREEZCH5)
For perennially ice-covered lakes under hydrological steady-state condition, the annual rate of water freezing at the bottom of the ice cover equals the annual ablation rate at the ice surface (Wilson and Wellman, 1962;McKay et al., 1985). This implies that to maintain the lake water volume constant ( V = 0), the annual accretion and ablation rates of the ice cover must be equal. We developed a recursive version of FREEZCH5 (FREEZCH9) that maintains V = 0 by mixing the evolved chemistry and δD-δ 18 O composition of the residual water with an input of water of fixed chemistry and δD-δ 18 O composition, which then produces the new input chemistry and δD-δ 18 O for the next loop. The hydrated precipitates that may form are accumulated in the series. Supplementary Appendix A provides details of FREEZCH9.

Field Sampling and Laboratory Analysis
The water column in the oxic basin of Lake Untersee was sampled in December 2016 at 10-m depth interval using a 2.2-L acrylic Kemmerer bottle content and samples were stored in sealed 20 ml high density polyethylene bottles (HDPE) bottles. The ice cover was sampled at the same location in December 2017 with a Cold Regions Research and Engineering Laboratory coring kit to a depth of c. 2 m to prevent mixing with lake water. Surface samples from the Anuchin Glacier were collected with a Kovacs ice corer powered by a Bosch Hammer GBH 36V-LI PLUS electrical drill, to a depth of 50 cm (e.g., Weisleitner et al., 2019). Samples from the Anuchin Glacier and the ice cover of Lake Untersee (sliced at 2 cm intervals) were melted in sealed Ziploc bags, transferred in HDPE bottles and shipped in coolers to the University of Ottawa where they were stored at 4 • C until analyzed for stable water isotopes.
The 18 O/ 16 O and D/H ratios of lake water and melted ice samples were analyzed using a Los Gatos Research liquid water analyzer coupled to a CTC LC-PAL auto sampler and verified for spectral interference contamination. The results are presented using the δ-notation (δ 18 O and δD), where δ represents the parts per thousand differences for 18 O/ 16 O or D/H in a sample with respect to Vienna Standard Mean Ocean Water (VSMOW). Analytical reproducibility for δ 18 O and δD was ±0.3 and ±1 , respectively. Deuterium-excess (d) was then calculated using the following equation (d = δD -8 δ 18 O; Dansgaard, 1964).

Modeling the Freezing of Low to High Salinity Waters
The freezing point of water is determined by the chemical composition of the solution and ranges from 0 • C for pure water to −2.3 • C for sea water (Figure 2A). Freezing under equilibrium conditions progressively decreases the amount of residual water while increasing the concentration of solutes in the residual water. This cryo-concentration of solutes on decreasing residual water is dependent on the initial chemical composition of the solution where for nearly pure water the bulk of the water is frozen over the 0 to −4 • C range, but for more saline waters, c. 15-20% of residual water may occur over the −5 to −20 • C range ( Figure 3A). For brackish to saline waters, hydrated salts can precipitate (NaCl·2H 2 O at −23 • C, MgCl 2 ·12H 2 0 at −36 • C) and the formation of these precipitates further reduces the amount of residual water (Figure 2A).
The dependence of the initial chemical composition of the solutions on the amount residual water is reflected in the δ 18 O of the residual water ( Figure 2B). For nearly pure waters, the δ 18 O of the residual water rapidly decreases as most of the water freezes over a narrow temperature range. For example, the freezing of nearly pure water with an initial δ 18 O composition of −35 over the temperature range of 0 to −4 • C will cause the δ 18 O of the residual water to rapidly decrease to value of −45 (when residual water fraction is 0.05). In contrast, the freezing FIGURE 2 | Evolution of (A) residual water fraction and (B) δ 18 O of residual water during freezing as a function of various molalities of sea water; (C) relation between δ 18 O and residual water for various molalities of sea water; and (D) δD-δ 18 O evolution of liquid water and ice during freezing, using 100 and 10% sea water salinity solutions (circles = water; squares = ice). See Table 1 for values of input parameters used in modeling.
of sea water with the same initial δ 18 O composition over the temperature range of 0 to −4 • C will cause the δ 18 O of the residual water to decreases to −37 (when residual water fraction is 0.85). A value of −45 in the residual water of sea water would only be reached following the precipitation of NaCl·2H 2 0.
These numerical simulations allow us to explore the effect of temperature on residual water content and its δ 18 O composition during equilibrium freezing of solutions with varied salinities. However, exploring the evolution of δ 18 O of residual water as a function of decreasing fractions of residual water revealed no significant differences between solutions of varying salinities ( Figure 2C). As such, a Rayleigh-type isotope fractionation during equilibrium freezing of liquid water with varied salinities using FREEZCH5 shows that, as freezing begins, the first ice to form will be enriched in δ 18 O by about 3% (and about 20% in δD), but as freezing continues, the δ 18 O composition of the ice formed progressively becomes depleted as the δ 18 O composition of the residual water trends toward lower value. The extent of δ 18 O depletion of the residual water during freezing is constrained by the evolving eutectic temperature of the residual water and temperature at the ice-water interface. Despite little differences in the evolving δD-δ 18 O of residual waters of varied salinities, a difference is observed in the δD-δ 18 O composition of the forming ice following the precipitation of hydrated minerals (Figure 2D). For perennially ice-covered lakes under hydrological steadystate, the evolution of δ 18 O in lake water and ice during equilibrium freezing was explored using FREEZCH9 for a range of recharge and mixing scenarios (Figure 3). The results are expressed in terms of residence times instead of annual loops to facilitate comparison. For all mixing scenarios, the results show that isotopic equilibrium (steady-state) in the lake water and ice is reached after more than three residence times. However, the difference between δ 18 O of input water and that of isotopic equilibrium decreases under increasing mixing scenarios. For mixing scenarios <∼5%, the δ 18 O input and equilibrium approaches ln(α 18 O i−w ); however, as the amount of annual mixing approaches 95%, the δ 18 O input and equilibrium approaches 0.

Testing FREEZCH9 With Lake Vostok, East Antarctic Plateau
Lake Vostok is a large subglacial lake located beneath c. 4 km of ice on the east Antarctic plateau (Lyons et al., 2016). The lake is about 250 km long and 80 km wide, has a surface area of ∼13,000 km 2 and a water volume of 4658-6350 km 3 depending on the estimation method (e.g., Studinger et al., 2004;Siegert et al., 2011;Li et al., 2019). Lake Vostok is recharged from the melting of the overlying East Antarctic Ice Sheet in the northern basin (incoming rate of c. 4 cm yr −1 ) and loses ice via ice accretion at the ice-water interface in the southern and central basins (freezing rate of c. 2 cm yr −1 ) (Siegert et al., 2000;Bell et al., 2002;Studinger et al., 2004). The annual input of water is estimated at 0.19 km 3 yr −1 , which translates to a residence time equivalent to c. 28,000 years or an annual mixing rate of c. 0.0037% (Royston-Bishop et al., 2004).
Based on the average of δ 18 O time-series during the last four glacial cycles in the Dome B record, Royston-Bishop et al. (2004) first estimated the δD-δ 18 O of input water to Lake Vostok to be −442.2 and −56.7 , respectively. Lipenkov et al. (2016) later revised the δD-δ 18 O composition of input water, based on data from the overlying accreted ice, to be −442.4 and −56.2 , respectively, with D-excess = 7.2 . Based on a closedsystem isotopic equilibrium model using the δD-δ 18 O data of re-frozen lake water in borehole, Ekaykin et al. (2016) estimated the average δD and δ 18 O composition of the Lake Vostok to be −455 ± 1 and −59.0 ± 0.3 , respectively (D-excess of 17 ± 1 ).
FREEZCH9 was used with input parameters defined as 0.0037% annual mixing (e.g., Royston-Bishop et al., 2004), the average major ions concentrations in the Vostok ice core which extends to 160 kyr BP (Legrand et al., 1988), and the equilibrium factors (α i−w ) of O'Neil (1968) and Suzuoki and Kimura (1973) also corrected for initial isotope composition of water. When using the δD-δ 18 O of the input water estimated by Royston-Bishop et al. (2004) (−442.2 and −56.7 , respectively), the model does not perfectly match the measurements in the accreted ice (Figure 4). However, if we use the α i−w corrected for initial isotope composition of water and adjust the δD-δ 18 O of input water to match the measurements in the accreted ice, we find input values for δD-δ 18 O of −443.5 and −56.4 , FIGURE 3 | (A) δ 18 O evolution of lake water and forming ice for various annual mixing scenarios as a function of number of residence time (initial δD-δ 18 O values are -270 and -35 , respectively). The colored lines represent different annual mixing rate scenarios (0.1-50%); (B) δD-δ 18 O plot of evolving lake water and ice as a function of the annual mixing rates of well-sealed perennially ice-covered lakes. The red dots indicate the δD-δ 18 O values of initial waters and accreted ice; (C) effect of annual mixing rates in perennially ice-covered lakes on the difference between the initial and steady-state δ 18 O values of lake water; and (D) effect of annual mixing rates in perenially ice-covered lakes on the difference between the initial and steady-state D-excess values of lake water.
respectively. With 0.0037% annual mixing, the time for the water column to reach isotopic equilibrium is in the order of c. 2-3 residence times and the δD-δ 18 O of Lake Vostok is estimated to be −453.5 and −59.1 , within error of the values reported by Ekaykin et al. (2016).

FREEZCH9 and Lake Untersee, Dronning Maud Land, Antarctica
Lake Untersee is a perennially ice-covered lake located in the Gruber Mountains of central Dronning Maud Land, approximately 90 km southeast of the Schirmacher Oasis. The 6.5 km long and 1.5 km wide lake (surface area = 8.73 km 2 ), with a water volume of 5.21 × 10 8 m 3 , is located in a closed-basin dammed at its northern end by the Anuchin Glacier (Faucher et al., 2019). Lake Untersee has two subbasins: (1) a large well-mixed oxic basin at the northern and central section with a maximum depth of 169 m; and (2) a shallower basin at its southern section with a maximum depth of 100 m depth (Wand et al., 1997;Andersen et al., 2011). Unlike most perennially ice-covered lakes, Lake Untersee does not develop a summer moat. Loss of water is entirely through the sublimation of the ice cover (61 ± 11 cm yr −1 ) with annual loss of 4.90 ± 0.8 × 10 6 m 3 of water (or 0.94 ± 0.17% of total lake water volume) (Faucher et al., 2019). To maintain hydrological steady-state, Lake Untersee receives an equal volume of inflow sourced from the subaqueous melting of the Anuchin Glacier (40-45% of annual contribution) and from subglacial meltwater (55-60% of annual input) (Faucher et al., 2019).  Table 1 for values of input parameters used in modeling.
The oxic water column of Lake Untersee has a narrow δDδ 18 O composition, from −288.4 to −287.3 and from −37.9 to −37.6 , respectively (Figure 5). With D-excess ranging from 12.8 to 14.7 , the water samples are distributed above the AMWL (Figure 5).The ice cover sampled at the same location has higher δD-δ 18 O values varying from −279.5 to −268.5 and from −36.2 to −33.55 , respectively, with lower D-excess ranging from −8.8 to 11.4 (Figure 5). FREEZCH9 was used with input parameters defined as 1% annual mixing (Faucher et al., 2019), α i−w of O' Neil's (1968) and Suzuoki and Kimura's (1973) corrected for initial isotope composition of water and the δD-δ 18 O of input water was adjusted to match the δD-δ 18 O measurements in the oxic water column. The initial chemistry of the input water was based on data from nearby shallow firn core (Isaksson et al., 1996). The FREEZCH9 modeling indicates that the measured δD-δ 18 O values in the water column is in isotopic equilibrium with initial input water having δD-δ 18 O values of −275.0 and −35.0 , respectively. With 1% annual mixing, the time to reach isotopic equilibrium in the water column is in the order of c. 300-500 years.

δD-δ 18 O Composition of Input Waters to Lake Untersee
Based on a hydrological and energy mass balance, it was suggested that the well-sealed perennially ice-covered Lake Untersee is recharged by two sources: subaqueous terminus melting of the Anuchin Glacier (40-45% of annual contribution) and from subglacial meltwater and/or from groundwater (55-60% of annual input) (Faucher et al., 2019). The FREEZCH9 modeling using 1% annual recharge and mixing with lake water suggests that the δ 18 O composition of Lake Untersee (−37.9 to −37.6 ) is in equilibrium with source waters having an average δ 18 O composition of −35.0 . The δ 18 O composition of  Table 1 for values of input parameters used in modeling.
the contribution from the subaqueous melting of the Anuchin Glacier can be estimated from the thickness of the glacier at the lake-glacier interface (c. 160 m) and the age-depth relation from the Dome Fuji, Dome C, and Taylor Dome ice cores (Lorius et al., 1979;Steig et al., 2000;Uemura et al., 2012). The age of the ice in the uppermost c. 160 m at these three sites varies from 0 to 3000-4500 years ( Figure 6A); assuming a similar age-depth relation applies for the Dronning Maud Land region, the average age of the subaqueous meltwater recharging Lake Untersee is estimated at c. 1800 years. The δ 18 O composition of a series of shallow ice cores (uppermost 50 cm) of a 3.2 km long transect on both sides of the medial moraine on the Anuchin Glacier ranged from −36.4 to −30.8 (average of −34.2 ± 1 ). Given that over the past 3000-5000 years, the δ 18 O time-series of deep ice cores showed fluctuations in the order of 0.5 (Figure 6B), the average δ 18 O of near-surface ice of the Anuchin Glacier likely reflects that over a 160 m thick section (or past c. 1800 years). Therefore, the δ 18 O of Lake Untersee water is in equilibrium with subaqueous melting of the Anuchin Glacier having an average δ 18 O composition of −34.2 ± 1 and to obtain the calculated δ 18 O of input water of −35.0 , a two-component mixing model indicates that the subglacial meltwater has an average δ 18 O composition of −35.8 ± 1 (δD: −280 ; D-excess: 6 ). Considering that the δ 18 O of last glacial maximum to Holocene age glacial ice is in the order of 6-8 , the subglacial meltwater recharging Lake Untersee is likely of similar age (late Holocene) to that of the Anuchin Glacier and does not reflect subglacial meltwater sourced during the last glacial period or far in the Antarctic interior.
The oxic water column of Lake Untersee is well-mixed as a result of buoyancy-driven convection which allows for effective mixing within a time-scale of one month (Steel et al., 2015). Based on the volume of water in the lake, its surface area and ablation rate, the residence time of water in the lake was estimated at c. 100 years (Faucher et al., 2019) and δ 18 O steady-state is achieved within 300-500 years (3-5 residence times). Considering that the 2016 water samples in Lake Untersee plot above the AMWL with D-excess of 14 , this suggest that for at least the past 300-500 years: (1) subaqueous melting of the Anuchin Glacier and subglacial meltwater were the only sources of water recharging the lake (i.e., no to very little recharge from surface streams with an evaporated δ 18 O signature); (2) the δD-δ 18 O composition of the subglacial meltwater is likely distributed on the AMWL and does not reflect a subglacial meltwater that has evolved under cryo-freezing; and (3) the ice cover has remained wellsealed to atmospheric exchange (i.e., no evaporative loss of the water column). These conditions were likely made possible by the cooler air temperatures in the region since the 1800s (Isaksson et al., 1996).
Despite generally cooler temperatures over the past few hundred years, we can test for the effect of an anomalous warm summer on the potential contribution of surface stream to the δD-δ 18 O composition on Lake Untersee. As a hypothetical scenario, we mixed 5-10% input of stream having a δD = −233.0 and δ 18 O = −23 (a scenario that represents the Anuchin Glacier δD-δ 18 O composition that evaporated 40% under relative humidity of 50%) with Lake Untersee and the δD-δ 18 O composition of the lake evolved outside the range measured in the oxic water column (10%mix: δD = −281.5 δ 18 O = −36.4 ; 5%mix: δD = −284.7 δ 18 O = −37.1 ). If Lake Untersee then continues to be under hydrological steadystate and keeps being recharged by subaqueous melting and subglacial meltwater sources having an average δ 18 O composition of −35.0 , the δD-δ 18 O of Lake Untersee will evolve to the isotopic steady-state determined by input and recharge over a period of 200-300 years. Therefore, this modeling exercise shows FIGURE 6 | (A) Ice cores age-depth relations for the Dome Fuji (Uemura et al., 2012), Dome C (Lorius et al., 1979) and Taylor Dome (Steig et al., 2000). The dotted line indicates the thickness of the Anuchin glacier at the lake-glacier interface (160 m that Lake Untersee most likely did not receive additional input from streams under an anomalous warm summer as its δD-δ 18 O composition would be distributed closer to the AMWL.

Application of FREEZCH9
This study developed an isotope-augmented version of FREZCHEM (FREEZCH5) and made it recursive (FREEZCH9) to predict the δD-δ 18 O evolution and isotopic steady-state of well-sealed perennially ice-covered lakes and subglacial lakes. The model accounts for changes in the ice-water isotope fractionation factor as freezing increases the salinity of the residual water and, to maintain the lake volume in equilibrium, mixes some input water of fixed chemistry and δD-δ 18 O composition to the residual water in the lake. Lerman (1979) suggested that the δD-δ 18 O steady-state of freezing lakes is reached after more than three residence times. Although the modeling reaches the same conclusion, the δD-δ 18 O composition of freezing lakes is not only determined by the residence time of water but also by the volume of ice that freezes annually which determines the volume of water that is added and mixes in the lake to maintain hydrological equilibrium. For mixing scenarios <∼5%, the difference between δ 18 O input and equilibrium in the lake water approaches ln(α 18 O i−w ); however, as the amount of annual mixing approaches 95%, the δ 18 O input and equilibrium approaches 0. In terms of deviation from the GMWL, the difference between D-excess of input and equilibrium in the lake water approaches 8 for annual mixing scenarios <5% and is reduced to near 0 as the amount of annual mixing approaches 95% ( Figure 3D).
FREEZCH9 allows one to test for the effects of input water with differing isotopic signatures (i.e., evaporated water, meteoric water, cryo-freezing water) on the distribution of lake water isotopic composition in δD-δ 18 O space (Figure 7). Horita (2009) suggested that lakes with δD-δ 18 O composition plotting above FIGURE 7 | Effects of water input with various δD-δ 18 O signatures (i.e., evaporated water, meteoric water, cryo-freezing water) on the δD-δ 18 O composition of well-sealed ice-covered lakes. For this example, the mixing rate is 1% and isotopic steady-state is reached in c. 300 years.
the GMWL are reflective of being affected by the freezing process. The residual waters in well-sealed ice-covered lakes and subglacial lakes that are recharged by meteoric-like waters (i.e., plotting along the GMWL) should plot above the GMWL, regardless of their annual mixing rate. The extent to which the D-excess of these lakes will deviate from the GMWL will be a function of their annual mixing rate. This is the case for Lake Vostok (0.0037% annual mixing) where the δ 18 O of recharge water (−56.7 ) approximates that in the accreted ice (−56.2 ) and the water column is depleted by c. 3 (−59.0 ± 0.3 ) but with D-excess (17 ± 1 ) c. 10 higher than estimated for the input water (7.2 ).
In the case of well-sealed perennially ice-covered lakes or subglacial lakes recharged by surface or soil waters (e.g., ephemeral streams) and whose δD-δ 18 O composition is shifted below the GMWL due to evaporation, the δD-δ 18 O of the water column may not be shifted above the GMWL following the freezing process (Figure 7). The latter would be dependent on the extent of evaporation and mixing condition in the lake.
Finally, in the case of well-sealed perennially ice-covered lakes or subglacial lakes recharged by subglacial waters that have experienced freeze-concentration, the δD-δ 18 O of these lakes will strongly deviate above the GMWL for any mixing scenarios. That is because the recharge waters will have δD-δ 18 O composition plotting above the GMWL and the freezing process taking place in the lake will shift the δDδ 18 O of the lake water further above the GMWL (but within the limits shown in Figures 3C,D). Such lakes may have D-excess that evolved to values greater than 17 . It has been suggested that there might exist a regional groundwater system beneath the Antarctic Ice Sheet, as inferred for Lake Bonney in the McMurdo Dry Valleys from airborne transient electromagnetic sensor (i.e., Mikucki et al., 2015). If this regional groundwater system includes circulation through subglacial lakes, this isotopic signature could be used to infer the source of water (largely from basal melting or multiple freeze-concentration).

CONCLUSION
In this study, we developed an isotope-augmented version of FREZCHEM (FREEZCH5) and made it recursive (FREEZCH9) to predict the δD-δ 18 O evolution and isotopic steady-state while accounting for changes in the ice-water isotope fractionation factor as freezing increases the salinity of the residual water. The following conclusions can be drawn from our results: 1. Unlike previous models, FREEZCH9 takes into consideration the combined effects of residual water mixing with incoming waters and evolving geochemistry on the stable water isotopic composition of well-sealed ice-covered lakes; it can be a tool to estimate the annual mixing fraction and residence time of the lake. 2. FREEZCH9 modeling of Lake Vostok suggests that the δDδ 18 O composition of the upper water column should be −453.5 and −59.1 , respectively; these are within error of the calculations of Ekaykin et al. (2016) (−455 ± 1 and −59.0 ± 0.3 for δD-δ 18 O, respectively). 3. FREEZCH9 modeling of Lake Untersee using 1% annual recharge and mixing with lake water suggests that the δ 18 O composition (−37.9 to −37.6 ) is in equilibrium with source waters having an average δ 18 O composition of −35.0 . The lake is recharged from subaqueous melting of the Anuchin Glacier and subglacial meltwater. Based on the δ 18 O age-depth relation, the δ 18 O of the subaqueous melting of the Anuchin Glacier is estimated at −34.2 and two-component mixing indicates that the δ 18 O of subglacial meltwater is −35.8 . The distribution of δD-δ 18 O of the oxic water above the AMWL with D-excess of 14 suggest that the lake has been well-sealed for the last 300-500 years and not been fed by surface waters (recharged only by subaqueous melting of the Anuchin Glacier and from subglacial meltwater). 4. The distribution of well-sealed ice-covered lakes in δDδ 18 O space depends on the input water (i.e., evaporated water, meteoric water, cryo-freezing water) having different D-excess signatures. Hence, FREEZCH9 may be also used to determine if the hypothetical groundwater systems of the McMurdo Dry Valleys (i.e., Mikucki et al., 2015) are fully or partially recharged by subglacial lakes.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
BF: manuscript writing (main contributor), data collection, data analysis and interpretation, and editing prior to manuscript's submission. DL: manuscript writing, data collection, data interpretation, and editing prior to manuscript's submission. DF: model development and editing prior to manuscript's submission. KW and DA: data collection and editing prior to manuscript's submission.