Earlier ice melt increases hypolimnetic oxygen despite regional warming in small Arctic lakes

Although trends toward earlier ice‐out have been documented globally, the links between ice‐out timing and lake thermal and biogeochemical structure vary spatially. In high‐latitude lakes where ice‐out occurs close to peak intensity of solar radiation, these links remain unclear. Using a long‐term dataset from 13 lakes in West Greenland, we investigated how changing ice‐out and weather conditions affect lake thermal structure and oxygen concentrations. In early ice‐out years, lakes reach higher temperatures across the water column and have deeper epilimnia. Summer hypolimnia are the warmest (~ 11°C) in years when cooler air temperatures follow early ice‐out, allowing full lake turnover. Due to the higher potential for substantive spring mixing in early ice‐out years, a warmer hypolimnion is associated with higher dissolved oxygen concentrations. By affecting variability in spring mixing, the consequences of shifts in ice phenology for lakes at high latitudes differ from expectations based on temperate regions.

With warming climate, lake ice-out now occurs earlier in many temperate and higher latitude regions (Magnuson et al. 2000).Ice-out phenology affects the timing of energy and mass inputs into lakes, ultimately shaping the relative duration of spring mixing and summer stratification and affecting the seasonal progression of biogeochemical and ecological processes (Cavaliere et al. 2021).Although trends toward earlier ice-out have been documented in lakes across the world, the links between ice-out timing and lake thermal and biogeochemical structure remain unclear, especially for lakes at higher latitudes that may respond differently compared to lakes in temperate regions (Klanten et al. 2023).Moreover, even though ice-out occurs earlier on average, there remains large interannual variability (Brown 2022); hence, it is important to understand lake responses along a gradient of ice-out timing.
The Arctic is warming three to five times faster than the rest of the globe, with the largest amplification happening during winter months and spring melt season (Rantanen et al. 2022), leading to nonlinear, abrupt shifts across lake ecosystems (Saros et al. 2022).One of the crucial differences between lakes in temperate and high-latitude regions that drive responses to ice phenology is that ice-out at higher latitudes occurs later in the season, close to the peak intensity of solar radiation (MacIntyre et al. 2009).As a result of the higher rate of incoming radiation at the time of ice-out, small lakes warm rapidly and experience incomplete or extremely short spring mixing periods that are largely invariable from year-to-year (Cortés and MacIntyre 2020;Pilla and Williamson 2021).When spring meromixis occurs and lakes effectively become monomictic, solutes accumulated during winter, cold water, and low oxygen conditions in lake bottoms persist into the summer (Findenegg 1937).However, with earlier ice-out, the duration of lake mixing may increase due to lower solar radiation, shorter days, and higher potential of exposure to weather events in the spring.
Long-term trends in lake thermal and biogeochemical structure are heterogeneous (Pilla et al. 2020), underscoring the need to assess the role of changing seasonal variability.Here, we used long-term monitoring data from 13 oligotrophic lakes in West Greenland to investigate (a) the relative effects of inter-annual variability in ice-out timing and weather conditions (air temperature, wind speed) on lake thermal structure and (b) how variability in thermal structure affects summer dissolved oxygen (DO) concentrations.The aim was to refine our understanding of the mechanisms by which ice phenology affects lake ecosystems and disentangle different seasonal signals reflected in interannual variability of lake responses.We expected that ice-out timing and weather immediately following ice-out will strongly affect temperature and DO in the hypolimnion as they are linked to spring mixing (Dokulil et al. 2006), while epilimnetic temperatures will be more strongly driven by summer weather conditions (Kettle et al. 2004).

Site
We analyzed thermal structure shifts in 13 lakes located in West Greenland along a $ 50 km gradient between Kangerlussuaq (67.1 N, 50.7 W; Supporting Information Fig. S1) and the Greenland Ice Sheet.These lakes are small, chemically dilute, glacial kettle lakes (Fig. S1B).Seven lakes are located inland and six are near the ice sheet (Supporting Information Table S1).Summer (June-July) daily mean temperatures vary across the landscape gradient, with 7-14 C inland and 5-9 C near the ice sheet.Precipitation is approximately 150-200 mm annually, with sublimation and evaporation playing a significant role in local hydrology (Johansson et al. 2015).Data from two different weather stations were used to characterize conditions inland and by the ice sheet (Supporting Information).North Atlantic Oscillation (NAO) data  were accessed from the National Oceanic and Atmospheric Administration (https://ftp.cpc.ncep.noaa.gov/cwlinks/).

Ice-out model
Lake ice-out, the first day of year when lakes became completely ice-free, was determined by personal observations, high-frequency temperature sensors, and visual inspection of available remote sensing imagery (Landsat 4-8 and Sentinel 1-2, cover the period between 1982 and 2008; however, the record is patchy due to the lack of cloudless imagery during spring).To hindcast ice-out timing using weather data, we compared multiple linear regression models using Akaike Information Criterion (AIC) with 0 C air temperature isotherm date, air temperature and wind speed in the months preceding ice-out, and rate of spring warming as predictors (Supporting Information Table S2).The best model included the rate of spring warming from April to May, better capturing the dynamic weather conditions during spring (Supporting Information Fig. S3).
To provide a historical perspective of climate shifts in the West Greenland region over the past eight decades, we analyzed spring NAO, and air temperature time-series.Previous work in the region has established the importance of abrupt changes in the local climate and environmental time-series (Saros et al. 2019).Hence, we employed breakpoint analysis to identify shifts in the mean of each time-series using function breakpoints (strucchange, Zeileis et al. 2002Zeileis et al. , 2003)).Moreover, we were especially interested in the tendencies of ice-out timing beyond shifts in its mean.To assess trends in early and late ice-out timing, we used quantile regression (quantreg, Koenker 2023).

Thermal structure analysis
Thermal structure responses to changes in ice-out timing were primarily modeled based on 123 vertical profiles collected with a multiprobe (HydroLab, until 2015, YSI EXO from 2016, YSI EXO and Turner C3 in 2019) from 2011 to 2022 (no observations in 2012).However, considering that vertical profiles only represent snapshot conditions, we also separately analyzed 44 high-frequency datasets (Onset Hobo Pendant) from a subset of 7 lakes collected in years 2003 and 2014-2022 (Supporting Information Fig. S1B; Table S3).For each vertical profile, we calculated five different thermal metrics: (a) volumetrically averaged epilimnion temperature; (b) volumetrically averaged hypolimnion temperature; (c) temperature difference between (a) and (b); (d) mixing depth; and (e) Schmidt stability (rLakeAnalyzer, Supporting Information, Winslow et al. 2019).For each time series, we aggregated water temperatures from 19 June to 19 July (the most frequently sampled open-water period across lakes) to daily means and calculated daily thermal structure metrics (a-e) and averaged them over the sampled summer period.
To determine the effects of ice-out timing and weather conditions on each calculated thermal metric (a-e), we used hierarchical linear mixed models (function lme in nlme package; Pinheiro et al. 2021).For each thermal metric, we compared the performance of a null model, x $ 1 + (1jweather region/ lake), where x represents thermal metrics (a-e), and models that included combinations of ice-out timing and (a) weather conditions characterizing air temperatures and wind speed during summer or (b) immediately after ice-out (hereafter referred to as post-ice-out) as fixed effects.All predictor variables were z-scored so this approach allowed us to directly compare the effect sizes (β) of predictor variables (ice-out timing vs. summer or post-ice weather conditions) and identify their relative contributions to interannual variability in thermal structure metrics.Accounting for repeated, nonindependent weather (two weather stations, one for inland lakes and one for ice-sheet lakes) and lake observations, weather region and lake identity were incorporated as a random intercept term (1jweather region/lake).
Model performance was compared based on the AIC and marginal and conditional R 2 (MuMIN, Barto n 2020).Marginal R 2 (R 2 m ) characterizes variability explained by fixed effects Hazukov a et al.
Earlier ice-melt in Arctic lakes alone and conditional R 2 (R 2 c ) represents the total variability explained by both fixed and random effects.To account for temporal autocorrelation, we compared top models with and without first-order autoregressive structure (CAR1) using AIC and selected the most parsimonious model.Model diagnostics including variance inflation factors, standardized residual plots, and autocorrelation function plots were examined.To check the stability of coefficient estimates, we generated bootstrap (n = 1000) standard errors and confidence intervals (CIs) accounting for lake clusters (Supporting Information Table S4; Fig. S4).Performance of models was assessed by 10-fold crossvalidation (Table S4).Results of all models falling within AIC difference < 5 are presented in Supporting Information Tables S5 and S6.Here, we present results based on vertical profiles as more data are available from more lakes in a balanced design, but model outcomes using averaged highfrequency thermal structure data provide qualitatively the same results (Supporting Information Table S6; Fig. S5), buttressing our interpretation of thermal structure responses.

DO analysis
Changes in thermal structure shape vertical habitat gradients, including those of DO.We explored how variability in thermal structure metrics affects summer DO concentrations, modeling average epilimnion and hypolimnion DO concentrations as responses using each thermal structure metric (a-e) as predictors.Average epilimnion and hypolimnion DO concentrations were calculated from 116 vertical profiles measured by a multiprobe (HydroLab, until 2015, YSI EXO from 2016) for each lake-day.Here we used linear mixed models with random intercept and slopes (x j lake), where x is the predictor variable, as they allow for different magnitude of responses among lakes (Supporting Information Fig. S6).Models were selected as described above.To examine amonglake DO variability, we calculated Spearman rank correlations with mean dissolved organic carbon (DOC) concentrations and mean specific conductance (SC) for both epilimnion and hypolimnion.Vertical profiles of SC were measured by a multiprobe and DOC samples from epilimnion and hypolimnion were analyzed as in Saros et al. (2015).All data and code are available in Hazukova and Saros (2024).

Results
Partially in response to spring NAO fluctuations (Fig. 1A), both spring (April-May, Fig. 1B) and summer (June-July, Fig. 1D) air temperatures have increased in a nonlinear, abrupt, fashion over the past 80 yr.The timing of lake ice-out is driven by the rate of warming during spring (R 2 = 0.95, root mean squared error = 2.3 d).After the 1993 breakpoint (95% CI: 1985-2001, p = 0.001) median April to May temperatures increased from À1.8 C to 0.1 C and ice-out became on average 6 days earlier (Fig. 1C).While no temporal trends were detected in late ice-out (80 th percentile, F 1,118 = 1.5; p = 0.2), early ice-out is becoming earlier (20 th percentile, F 1,118 = 8.5; p = 0.004) (Fig. 1C).
When ice-out occurs early, in mid-May to late May (20 th percentile, day of year (DOY): 132-152), lakes are more likely to be exposed to relatively cooler air temperatures following ice-out (Fig. 2A) compared to when ice-out happens late during mid-June to late June (80 th percentile, DOY: 164-177) (Fig. 2A).In late ice-out years, lakes are exposed to median air temperature of 9.1 C with lower variability (interquartile Hazukov a et al. Earlier ice-melt in Arctic lakes range (IQR) = 5.7-12.5 C).In contrast, during early ice-out years, air temperature is much lower, 6.7 C, with year-to-year variability being almost twice as high (IQR = 1.6-11.7 C).Also, the daily average downwelling shortwave (Fig. 2B) and longwave (Fig. 2C) radiations increase by $ 25% and $ 15%, respectively, between mid-May and late June (Fig. 2B,C).Taken together, the likelihood of rapid onset of stratification increases when ice-out occurs late, closer to the peak solar insolation.
Epilimnetic temperatures increase with summer air temperatures (β = 1.54 AE 0.22) and in years with earlier ice-out (β = À1.00AE 0.18).The top model suggests an interaction between ice-out timing and summer air temperatures (β = 0.43 AE 0.19; R 2 c = 0.51, R 2 m = 0.54), muting the effect of air temperatures in early ice-out years, with on average warmer epilimnia (Fig. 3A).Model based on averaged thermistor data concurs with the results based on vertical profiles (R 2 c = 0.69, R 2 m = 0.95; Supporting Information Fig. S5A; Table S6).
Models based on vertical profiles suggest that mixing depth deepens with lower post-ice-out air temperature (β = À0.78AE 0.24), especially in late ice-out years (β = À0.30AE 0.19; R 2 c = 0.16, R 2 m = 0.37) (Fig. 3D).Models based on thermistor summer averages highlight the importance of ice-out timing (β = À0.88AE 0.17) and show that mixing depth becomes shallower with higher summer air temperature (β = À0.42AE 0.24, β int = À0.49AE 0.18, R 2 c = 0.42, R 2 m = 0.55; Supporting Information Fig. S5D; Table S6).Despite disagreement between the models, in both cases, early ice-out is linked to deeper mixing and the effect of air temperatures has the same direction: deeper mixing is associated with cooler air temperatures.

Discussion
Using a long-term monitoring dataset from 13 remote Arctic lakes, we provide evidence that ice-cover phenology is an important driver of summer thermal structure and DO variability in West Greenland lakes.Our findings align with predictions that ice phenology has a lasting influence on lake ecosystems in regions with long-periods of ice-cover (Dugan 2021) and arid conditions (Smits et al. 2020).However, we show that the impact of changing ice-out timing is not always linear and needs to be interpreted in the context of air temperatures.Hypolimnion temperatures are on average higher in early ice-out years, but the increase is amplified when early ice-out is coupled with lower post-ice-out air temperatures (Fig. 3B; Supporting Information Fig. S5B).Epilimnion temperatures are strongly driven by summer air temperatures, but when ice-out occurs early, the epilimnion is on average warmer and accumulates more heat given the simultaneous increase in mixing depth during early ice-out years (Fig. 3D; Supporting Information Fig. S5D; Warner et al. 2018).These results align with the expectation that surface waters will be more temperature and dissolved oxygen is positive in the hypolimnion.Lakes with higher DOC concentration exhibit supersaturated conditions in the epilimnion and a steeper increase along the temperature gradient in the hypolimnion.Thick red line shows the expected relationship based on solubility of oxygen in water along temperature gradient when at equilibrium with air at standard surface pressure of 1012 hPa, thinner red lines show the same relationship given pressure at 300 and 450 m of altitude (data based on Benson and Krause 1980;Mortimer 1981).Note that ranges of x-axes are different.
responsive to summer weather conditions while deep waters will be more strongly affected by the timing of ice-out and post-ice-out weather conditions, together shaping the extent of spring mixing (MacIntyre et al. 2009).
The counterintuitive association of higher DO concentration with a warmer hypolimnion provides further evidence of the importance of ice phenology for deep waters.During early iceout years when open-water conditions begin under a lower solar radiation regime, the likelihood of intensive turbulent mixing increases, allowing the deep waters to simultaneously warm up and re-oxygenate after winter (Fig. 5).Even though the best models did not include the effect of wind speed on hypolimnion temperatures, wind speed plays an important role in inducing mixing at ice-out, especially when the temperature difference at the air-water interface is lower.In Alaskan lakes, instabilities caused by wind-induced nonlinear waves have been linked to rapid transport of large amounts of heat into deep layers and deeper epilimnia (Cortés and MacIntyre 2020).In years with late ice-out, rapid warming following longer and more severe winters limits mixing and allows hypoxia and cold water to persist in the hypolimnion.The variability in DO responses among lakes is high but lakes with higher concentrations of DOC tend to respond more strongly along the hypolimnion temperature gradient (Fig. 4B).In the study lakes, DOC is strongly correlated with SC (ρ = 0.95) as a result of evapoconcentration (Anderson and Stedmon 2007), suggesting that lakes with higher concentrations of solutes may be forming chemical stratification under ice, making them more susceptible to spring meromixis (MacIntyre et al. 2018).
Considering that the greatest source of heat into the lakes is solar radiation, the timing of ice-out is a crucial process affecting exposure of lakes to energy fluxes, especially in extremely seasonal high latitude systems (Fig. 2; MacIntyre and Melack 2009).Moreover, trends in air temperatures vary seasonally and are highly heterogeneous across regions (Winslow et al. 2017).Hence, the implications of changing ice-cover phenology need to be considered relative to seasonality of solar radiation and climate trends that together with wind speed and humidity modify the heat budget (MacIntyre et al. 2009).For example, in West Greenland, warmer conditions during April and May that drive early ice-out (Fig. 1B,C) do not necessarily beget higher post-ice-out air temperatures (Fig. 2).In 2016 and 2019, years with the earliest ice-out on record (Fig. 1C), post-ice-out air temperatures were very different, changing the lake ecosystem conditions.In 2016, postice-out air temperatures were the lowest (1-4.5 C) since 2010, Fig. 5. Conceptual diagram demonstrating how timing of ice-out affects lake thermal and biogeochemical structure in high-latitude lakes with summer stratification.During early ice-out years, the likelihood of cooler, post-ice air temperatures that will enable spring mixing is higher.In contrast, during years when ice-out happens late, closer to the peak of solar insulation, lakes tend to stratify immediately after ice-out (sometimes even under-ice).Lack of full overturn in those years leads to the persistence of low oxygen conditions in hypolimnion.In years when ice-out occurs early, lakes tend to have higher epilimnion and hypolimnion temperatures, deeper mixing depths, and smaller differences between epilimnion and hypolimnion.
Hazukov a et al.
Earlier ice-melt in Arctic lakes leading to the highest summer hypolimnion temperatures on record (6.2-11.1 C) and high DO concentrations.In contrast, in 2019, warmer post-ice-out air temperatures induced more rapid warming that led to lower temperatures in hypolimnia (4.6-5.3C) and lower DO in the summer, suggesting limited spring mixing (Fig. 5).
Global predictions of lake responses to climate change suggest that a shift toward earlier ice-out will lead to prolonged stratification associated with negative consequences for deep water DO (Jane et al. 2021;Li et al. 2022), but our findings demonstrate that shifts to earlier ice-out may alleviate hypoxia stemming from limited spring mixing in Arctic lakes (Fig. 5).Given that spring meromixis is common among lakes across latitudes, we need to evaluate lake responses to ice phenology shifts more critically and consider weather conditions following ice-out that affect spring mixing.The interpretations of lake responses to icecover phenology are further complicated by winter precipitation: in regions where arid conditions during winter allow for the formation of clear lake ice, turnover and stratification may happen even before ice-out (Kirillin et al. 2021).On the other hand, when present, snowmelt intrusions can induce mixing during melt season (Cortés and MacIntyre 2020) and affect lake temperatures during summer (Christoffersen et al. 2008).To further our understanding of mixing processes in lakes experiencing spring meromixis, quantification of lake heat budgets will allow us to assess the duration and timing of spring turnover more accurately (Schmid and Read 2021).
Linking thermal structure metrics to specific seasonal climate signals will allow us to better disentangle long-term trends in lake responses to climate change.In West Greenland, median spring air temperatures exceeded 0 C (Fig. 1B), ice-out became earlier (Fig. 1C), and summers became warmer in the early 1990s (Fig. 1D).These recent, nonlinear shifts are linked to variability in large-scale circulation patterns expressed by the NAO index (Fig. 1A; Olsen et al. 2012) and the Greenland Blocking Index (Saros et al. 2019).Long-term trends in these indices suggest that climate is becoming more variable, bringing seasonally extreme air temperatures to the region with higher frequency (Hanna et al. 2022).However, the timing of these extremes matters.As we show, deep water temperatures carry a strong ice-out signal while surface waters respond more readily to summer air temperatures.Responses of thermal and biogeochemical structure to ice-out and air temperature variability will have further implications for the timing of ice-on (Oleksy and Richardson 2021), vertical habitat structure (Kraemer et al. 2021), and accumulation/outgassing of greenhouse gases from deep waters (MacIntyre et al. 2018).

Fig. 1 .
Fig. 1.Climate and ice phenology in Kangerlussuaq, West Greenland exhibit nonlinear trends and large variability since the beginning of the instrumental record in 1942 until 2022.Vertical gray lines depict breakpoints with 95% CIs shown as shaded gray boxes; solid lines show breakpoints with p < 0.05, and dotted lines show breakpoints with p = 0.07.White (or black in C) horizontal lines show averages before and after breakpoints.(A) Spring (MAM) NAO index with breakpoints demarcating a 1985-1994 period of consistently positive index bringing cold and dry weather to Greenland.Data only available since 1950.(B) Median April to May air temperatures shift from À1.75 C to 0.14 C after 1993.Green ribbon highlights the range of minimum and maximum temperature for each year and light gray horizontal line shows 0 C isotherm.(C) Hindcasted ice-out dates exhibit shifts in average ice-out timing from June 9 to June 3 (160-154 DOY).Crosses show observed iceout data, black points are modeled with gray ribbon showing 95% CI around the prediction.Lines in green and bordeuax represent quantile regressions (τ = 0.2, and 0.8, respectively), solid line designates significant trend ( p < 0.05).(D) Median June to July air temperatures shift from 9.4 C to 10.4 C after 1994.Bordeuax ribbon highlights the range of minimum and maximum temperature for each year and light gray horizontal line shows 0 C isotherm.Average minimum temperature has increased from À0.4 C to 0.7 C after 1994.Data between 1971 and 1973 are missing.

Fig. 2 .
Fig. 2. Comparison of the range of (A) air temperatures, (B) shortwave radiation, and (C) longwave radiation that lakes are exposed to when ice-out occurs early (20 th percentile, DOY: 132-152, green) and late (80 th percentile, DOY: 164:177, bordeuax) in Kangerlussuaq, West Greenland.Panels show the intra-annual variability with green and bordeuax vertical lines demarcating the period between the earliest and the latest ice-out on record.Each point shows a mean daily value and bars display standard deviation across years (A: 1993-2022, temperature data from the DMI station, B, C: 2011-2022, solar radiation data provided by the PROMICE/GAP program).

Fig. 3 .
Fig. 3. Critical thermal structure metrics, (A) volume-averaged epilimnion temperature, (B) volume-averaged hypolimnion temperature, (C) difference between epilimnion and hypolimnion temperatures, (D) mixing depth, and (E) Schmidt stability, show strong responses across the range of ice-out timing.Lines show how are responses to ice-out timing modulated by summer air temperatures (A, C, E; yellow-orange gradient spans temperature range of 5-15 C) and air temperatures immediately following ice-out (B, D; blue gradient spans temperature range of À1-15 C) as predicted by linear mixed models with summer vertical profile data as input (best model structure is described in each panel).Points show measured data with color corresponding to the temperature gradient modulating the response.Predictions and associated CIs are estimated based on models without autocorrelation structure.

Fig. 4 .
Fig. 4. (A) Mean dissolved oxygen concentration in the epilimnion decreases with increasing lake water temperature, while (B) the relationship between