The variability of weather conditions is expected to increase under ongoing climate change with more extreme events occurring, including, for example, heat waves, droughts, and intensive precipitation events (e.g., Diffenbaugh et al. 2013;Pendergrass et al. 2017;Cook et al. 2018). Extreme events have deleterious effects on ecosystem goods and services such as storm surges (e.g., Karim and Mimura 2008) or decreasing food security (e.g., Thornton et al. 2014). Similarly, phenological observations in lakes, such as the timing and duration of lake ice cover have been predicted to increase in variability under climate change (e.g., Weyhenmeyer et al. 2011). However, phenological changes cannot continue interminably as a new stable state might be reached, that is, lakes might turn from being ice-covered to becoming ice-free (Sharma et al. 2019). Increasing variability may provide an early-warning signal for reaching a new stable state (Scheffer et al. 2009). Thus, documenting changes in the variability of ice cover is critical for understanding how lakes are responding to climate change (Rühland et al. 2023), as ice on lakes plays an important role in numerous physical and ecological lake processes in winter and throughout the rest of each year (Hampton et al. 2017;Hébert et al. 2021;Jansen et al. 2021).

Changes in lake ice phenology (timing of ice-on and iceoff) have shortened lake ice duration over the last century because of climatic variation (Magnuson et al. 2000;Newton and Mullan 2021). Despite the consistent decrease in ice duration in lakes around the world, year-to-year variability in the length of ice cover remains high (Duguay et al. 2006;Wang et al. 2012) with linear trends explaining < 30% of the overall variation (e.g., Wynne 2000;Benson et al. 2012).

The extreme ice seasons could be driven by late freezes, early melts, multiple freeze-melt events, or even no ice cover at all (Bernhardt et al. 2012;Higgins et al. 2021;Sharma et al. 2021b). These extremes, including ice-free seasons, are predicted to increase dramatically in the future for individual lakes (Robertson et al. 1992;Magee and Wu 2017) and regions of lakes in the Northern Hemisphere (Sharma et al. 2021a;Wang et al. 2022). However, it is not yet clear which lakes are most sensitive to high interannual variability with the recent rapid increase in ice loss and which factors are driving interannual variability in lake ice (Brown and Duguay 2010).

Global anthropogenic climate change and teleconnections, large-scale climate linkages, can affect local and regional weather patterns, especially, air temperature, which is integrally related to lake ice (Ghanbari et al. 2009;Filazzola et al. 2020;Imrit and Sharma 2021). With synergistic interactions between climate change and teleconnections, extremes and interannual variability of air temperature are predicted to increase (IPCC 2021); thus, it is likely that the duration of ice cover will also become increasingly variable with periodicity related to teleconnections (Wang et al. 2012). In past research, the interannual variability of ice has been identified as predominantly increasing with shorter ice cover when examined at the annual, decadal, and 20-yr time scales (Kratz et al. 2000;Weyhenmeyer et al. 2011;Benson et al. 2012). One exception is that when broken into two 50-yr periods, ice duration variability decreased in many lakes, especially across Europe (Benson et al. 2012). Ice duration has a finite limit with the complete loss of ice, indicative of a nonlinear relationship that supports previous inconsistent results. Therefore, it is critical to understand the relationship between ice duration and variability when trying to understand and predict the response of lake ice to global drivers of regional weather like climate change and teleconnections.

Here, we explored patterns and drivers of lake ice variability in 43 Northern Hemisphere lakes over the last 87 yr, using a recently compiled database on lake ice phenology (Sharma et al. 2022). We define interannual variability in ice as the calculated standard deviation or variance of ice phenology duration over a series of years in a single lake. We asked three main questions: (1) What patterns emerge when examining the trends in ice variability over the past 87 yr?; (2) Is there a consistent relationship between aspects of ice phenology (ice-on, ice-off, and duration) and the variability observed in ice phenology across different lakes?; and (3) To what extent can climate anomalies and teleconnections, recognized as global drivers of regional weather, explain the fluctuations in ice duration amidst the observed decreasing ice trends? We hypothesized that the interannual variability of ice phenology no longer significantly increases if ice duration becomes too short, following a nonlinear relationship. The hypothesis implies that lake in lakes in colder geographic regions would experience increasing interannual variability while lakes in warmer geographical regions will experience a decrease in interannual variability. We also hypothesized that warmer global temperatures in the Northern Hemisphere winter and teleconnection indices, such as North Atlantic Oscillation (NAO) and El Niño-Southern Oscillation (ENSO), will significantly be related to the year-to-year variability in ice duration but with distinct geographical differences (Livingstone 2000;Ghanbari et al. 2009;Bai et al. 2012;Imrit and Sharma 2021).

Using a database of 78 lakes with ice phenology records extending over 100 yr (Sharma et al. 2022), we selected 43 lakes based on records that included ice duration with more than 65% of years with ice data, even if one or more winters were noted as ice-free (Supporting Information Table S1). These lakes were found between 42.50 N and 65.60 N latitude, spanning nine different countries (Supporting Information Fig. S1). We chose to examine records between 1931 and 2018 to encapsulate contemporary ice patterns in the Northern Hemisphere with a sufficiently long time series for as many lakes as possible (Supporting Information Table S1). Missing values for ice duration were uncommon in recent decades, although a few of the lakes were missing ice duration in the years typically surrounding world or local events (e.g., wars) that prevented data collection (Supporting Information Table S1; Sharma et al. 2022).

The ice phenology records included the duration of ice cover (in days), the geospatial coordinates of the survey point (latitude and longitude), the lake name, and the winter year of ice cover, that is, a lake that froze in January 2000 would be assigned the winter year of 1999 as winter encompasses two calendar years (i.e., 1999-2000). The database we used for ice phenology records also included information on lake morphometry, such as surface area, maximum lake depth, and elevation (Sharma et al. 2022).

We obtained the maximum winter air temperatures for December, January, and February (DJF) from the Climatic Research Unit (CRU) of East Anglia (Harris et al. 2020), which were downscaled to 0.5 Â 0.5 grid cells. We acknowledge that the available climate data has limitations in terms of resolution, which may result in lakes that are close together having the same temperature value. However, we selected the CRU dataset as having the finest spatial resolution while also providing annual climate patterns. Monthly temperature values were extracted for each year at every lake where data on ice duration was available, including years with no ice present. We obtained global climate and teleconnection indices monthly for October through May, spanning the time frame of ice cover from the lakes in this dataset. Global annual temperature anomalies (GTA) were obtained from the National Oceanic and Atmospheric Administration and averaged over land and ocean (NCEI 2023). We also considered two teleconnection indices as potential drivers of local winter weather conditions. We downloaded both NAO and ENSO monthly indices from the National Weather Service Climate Prediction Center (National Weather Service 2023).

We chose ice duration for these analyses because we could appropriately quantify ice duration when a lake did not freeze (ice duration = 0 d), which is not possible with ice-on or iceoff dates when a lake did not freeze. First, for visualization, we calculated a 10-yr moving average and Bollinger Bands, one rolling 10-yr standard deviation above and below the moving average, that can indicate the volatility of a time series (Bollinger 1992). We used standard deviations to quantify variability patterns in ice duration. We applied 10-yr rolling standard deviations to account for variations included in major climate oscillations and teleconnection patterns that happen periodically (Sharma and Magnuson 2014;Imrit and Sharma 2021). All analyses and visualizations were completed using R version 4.1.2 (R Core Team 2022) for this section and the rest of the manuscript.

While simple moving averages and rolling standard deviations can help understand trends, the overlapping nature of the rolling windows results in high autocorrelation. In addition, choosing a single window for calculating variability can result in different conclusions (e.g., Benson et al. 2012). As an alternative, we identified all sequential windows between 4 and 30 yr in length (26 versions of sequential windows), starting with 2018 and moving backward to 1931. For example, a 16-yr sequential window would encapsulate nonoverlapping sets of 16 yr (e.g., 2018-2003; 2002-1986), while a 4-yr sequential window would encapsulate non-overlapping sets of 4 yr (e.g., 2018-2015; 2014-2011). For each sequential window, we required a minimum of 75% of years having duration data; for those windows, we calculated the mean (hereafter, duration mean), standard deviation (hereafter, duration SD), and coefficient of variation (duration SD Â 100/ duration mean). We also calculated the year for each sequential window as the median of the start years in that window.

To determine whether a change in duration mean and SD occurred over the time series, we calculated linear models based on the duration mean or SD for each sequential window size. For example, with 10-yr windows, there would be up to 9 duration means and duration SD incorporated into the linear model. We used Theil-Sen median regressions (Komsta 2019) with the duration mean or duration SD as the response variable and median year as the predictor. We used a median-based regression because these methods are relatively robust to outliers, repeated measures, and changes in the distributions as the SD would become right-skewed with an increased number of years with no ice cover (Siegel 1982).

We calculated a slope for duration mean and SD for all sequential window sizes.

To determine which drivers related to trends in duration SD, we chose trends calculated with 17-yr sequential windows because 17-yr windows were the most represented when evaluating median trends in duration SD. We modeled trends in duration SD using generalized additive models (GAMs; Hastie and Tibshirani 1990;Wood 2017). We built candidate models based on ice characteristics, winter air temperature, geomorphometry, and geography established for each lake. For ice characteristics, we calculated the percent of ice-free years and the mean duration length in days for each lake. For winter air temperature, we used the annual average daily maximum temperature from DJF for each lake. Over all the years, we calculated the median DJF annual daily maximum temperature. We averaged across the three winter months to use the mean winter temperature for all analyses. We chose to summarize winter temperatures here to encapsulate the time period when most of these geographically and morphologically diverse lakes are frozen in a year. For geomorphometry, we used the surface area and maximum depth; both geomorphometry variables were log-transformed because of the several orders of magnitude spread (e.g., Lake Suwa is 7.6 m deep while Lake Baikal is 1642 m deep). For geography, we used latitude, longitude, and elevation. We fit increasingly complex GAMs using the "mgcv" package (version 1.8-40; Wood 2017) and ultimately selected the models that had statistically lower AIC and maximized deviance explained using the compareML function in the "itsadug" package (van Rij et al. 2022). We extracted all significant smooths for the selected GAM using the confint function in the "schoenberg" package (Simpson 2018), visualized the smooths using the "ggplot2" package (Wickham 2016), and arranged the plots with "patchwork" package (Pedersen 2022).

We examined the difference in variability between the two different ice phenology metrics (ice-on and ice-off) that are used to calculate ice duration. For each lake, we applied the Theil-Sen median regressions (Komsta 2019) for both ice-on and ice-off and calculated the residuals for each year. We used those residuals to calculate two overall variances (ice-on and ice-off) and compared those two variances using an F-test.

To determine the relationship between ice phenology and variability, we calculated the day of the year for ice-on and ice-off for each lake. We ignored years when the lakes did not freeze for the winter since there are no ice-on or ice-off dates recorded for that year. We used ice-on and ice-off means and SDs calculated for every lake for all sequential windows (n = 4-30 yr). To examine the shape of the relationship between mean and SD for each ice phenology variable, we fit GAM models (model: SD $ mean with k = 7 knots possible) using the "mgcv" package (version 1.8-40; Wood 2017) for each of the sequential window sizes (n = 4-30 yr). We assessed the effective degrees of freedom (edf), which reflects the degree of nonlinearity of a curve: edf = 1 indicates a linear relationship, edf up to 2 indicates a weak nonlinear relationship, and edf > 2 indicates a highly nonlinear relationship. We also assessed the mean ice-on or ice-off date when the GAM curve was at a maximum.

We hypothesized that the relationship between duration and interannual variability of ice phenology would follow a nonlinear Shepherd equation (eq. 1, fig. 1a, Shepherd 1982). To determine the relationship between ice duration and variability that matches our proposed hypothesis (Fig. 1a), we used duration means and duration SD calculated for every lake for all sequential windows (n = 4-30 yr). For each sequential window size, we fit a Shepherd function (Eq. 1) between variables for duration mean (mean window ) and duration SD (SD window ), which is the generalized form of the Michaelis-Menten function with three different parameters (A, B, and C) that permits the function to be domed or unbounded with a non-zero asymptote (eq. 1, Iles 1994). The Shepherd function appeared to be a good fit from the ice phenology GAM results, given that we could now include ice-free years (duration = 0 d). We estimated the three parameters using nonlinear least-squares estimates. We calculated the peak of the curve using the root of the first derivative of the Shepherd function and the inflection point using the root of the second derivative (Iles 1994). To match with the hypothetical groups proposed in Fig. 1a, we used a k-means clustering algorithm to identify clusters across all the individual sequential window sizes. We ran the algorithm for 1 cluster up to 9 clusters and examined the declining pattern of "within sums of squares" with an increasing number of clusters to look for an elbow, indicating that additional clusters have little added explanatory value (Tibshirani et al. 2001). Using the five identified clusters, we labeled each sequential window based on group (Fig. 1a).

We hypothesized that lakes would cluster into groups along the nonlinear relationship (Fig. 1). In lakes with no ice, interannual variability is 0; those lakes are consistently frozen (Fig. 1a,b: region i). Lakes in the warmest region with the shortest ice cover would experience decreasing variability (Fig. 1a,b: region ii). In slightly cooler regions, lakes would shift to high and stable variability (Fig. 1a,b: region iii). Lakes in colder regions would experience intermediate and increasing interannual variability (Fig. 1a,b: region iv). To identify which lake characteristics predicted each lake group located on the Shepherd function (Fig. 1), we selected the window size (16 yr) that was the best fit, based on AIC and R 2 , out of each of the Shepherd model fits. We selected the most recent 16-yr sequence (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) and identified the cluster assigned by cluster analysis for each lake. We used groups assigned for the five clusters as identified above (i, ii; iii; iv.1; iv.2; iv.3) and also used three groups (i, ii; iii; iv) to match Fig. 1a as categorical response variables. We used a regression tree with morphometric variables (max depth, surface area) and geography (latitude, longitude, elevation) to explain the assigned group. A parsimonious regression tree was selected by pruning the tree to the level where the complexity parameter minimized the cross-validation error. We calculated the percent variation explained by the regression tree (R 2 ) as R 2 = 1 À relative error (Sharma et al. 2012). Regression trees were completed using the "rpart" and "rpart.plot" packages (Milborrow 2019;Therneau and Atkinson 2019).

We examined the effects of global climate and teleconnection factors on year-to-year variability, measured as residuals from a Thiel-Sen slope line fit to all data (1931-2018) as above. Given the spatial distribution of our lakes, mostly in North America and Europe, and the timing of ice phenology, spanning October to May, we collapsed all three variables (GTA, ENSO, and NAO) to bimonthly averages for October/November (ON), December/January (DJ), February/March (FM), and April/May (AM) resulting in 12 unique predictor variables. We used these 12 variables scaled to bimonthly means to capture seasonal differences between variability in the timing of ice on our study lakes while also avoiding over-parameterizing models with too many explanatory variables. We removed 4 lakes with < 5 yr of non-zero ice cover as residuals were all close or equal to 0. For the remaining 39 lakes, we modeled the annual residuals of ice duration using GAMs with the same 12 explanatory climatic variables and fixed the number of basis functions for each smoothed term to 4 for each parameter. For each lake, we estimated GAMs using automatic parameter selection by penalizing each smooth using the "select = TRUE" option in the "mgcv" package (version 1.8-40; Wood 2017). We extracted all significant smooths for the selected GAM as above.

The duration of lake ice varied considerably among years and between lakes (Fig. 2; Supporting Information Fig. S2). The average duration of ice cover for the entire dataset was 112 d, ranging from a minimum of 0 to a maximum of 236 d (Supporting Information Table S2). Some lakes that were almost entirely ice-free for the duration of their time series had little to no interannual variation, such as Aergerisee (Fig. 2a). Lakes with a high frequency of ice-free years tended to have fluctuating standard deviations with many years close to 0, such as Greifensee (Fig. 2b). Other lakes had ice durations lasting around 2 months (e.g., Balaton, Fig. 2c) or longer ice durations lasting over 100 d (e.g., Otsego, Fig. 2d), both with less frequent ice-free years over the entire record.

Most lakes (79%) displayed decreasing duration means but trends in duration SD were less consistent looking across sequential windows of 4-40 yr (Fig. 3). Duration SD significantly increased for 49%, decreased for 7%, and had no significant trend for 44% of lakes (Fig. 3). Trends in the standard deviation of 17-yr sequential windows were best explained by ice characteristics, winter air temperature, and lake depth in a GAM that explained 85.5% of overall deviance (Fig. 4; Supporting Information Table S3). Lakes with no ice-free years had increasing trends in duration SD, but lakes with an increasing number of ice-free years were more likely to have decreasing trends in duration SD until the lake was ice-free all the time (Fig. 4a). Lakes with the coldest winter daily maximum air temperatures were more likely to have decreasing duration SD while approaching 0 C air temperatures indicated increasing trends in duration SD (Fig. 4b). When approaching 5 C, lakes were likely to have to change in duration SD (Fig. 4b). Deeper lakes had increasing trends in standard deviation (Fig. 4c). Finally, the trends in duration were most likely to switch from increasing to decreasing at an average of $ 100 d of ice cover (Fig. 4d). The trends in duration CV predominantly matched those of duration SD (data not shown) and therefore, we proceeded with using duration SD for the rest of the analyses.

Ice-on dates tended to have higher variability than ice-off dates. Ice-on variance was almost twice ice-off variance (F 2978,2936 = 1.82, p < 0.001; Fig. 5a). Later mean ice-on dates had higher ice-on SDs across all sequential windows with variability increasing by $ 40% across the range of mean ice-on dates (Fig. 5b). Ice-on SDs increased linearly with increasing ice-on mean (edf = 1) but some sequential windows had increasing quadratic or higher polynomial (edf ≥ 2) fits with maxima on 17 January. The GAM model explained 13% of the variance at most. Earlier mean ice-off dates had higher SDs across all sequential windows with variability increasing by 300% across the range of mean ice-off dates (Fig. 5c). The GAM model explained up to 77% of the variance; most of the fits were highly nonlinear (edf > 2). Maximum variance was on 16 February across all sequential windows and 23 February at the models with the downward tilt in early February (edf > 4.5) which were able to capture the Shepherd function-shaped curve proposed for ice duration (Fig. 1a).

We found a nonlinear relationship between duration standard deviation and average ice duration that was similar across all sequential windows (Fig. 6), and which supported our hypothesis (Fig. 1a). The median peak of all the models was at 26.0 d of ice duration while the median inflection point was 47.8 d (Fig. 6); this also represents the transition between increasing variability and decreasing variability (Fig. 1a). The inflection point of this relationship was at $ 1.5 months, at that boundary, there is a shift from accelerating (> 1.5 months ice duration) to decelerating (< 1.5 months ice duration) duration SD. The model with the best fit, as identified by deviance explained and AIC, was for 16-yr sequential windows (A = 474, B = 175, C = 1.7, R 2 = 0.75; Fig. 6b). Using all the data across all sequential windows and all lakes, k-means clusters were calculated for 1-9 clusters. Within sums of squares minimized at 5 clusters; therefore, we used 5 clusters to categorize each group of the duration mean vs. duration SD (Fig. 6b; Supporting Information Fig. S3). One cluster was identified at the lower end of ice duration; we labeled that as group i, ii to match with groups i and ii from the conceptual model (Fig. 1). Group iii matched the conceptual model, while group iv from the conceptual model was identified by the k-mean clustering as three distinct clusters, we labeled those as groups iv.1, iv.2, and iv.3 according to increasing ice duration (Fig. 6b) and also lumped all of those iv categories together to match our hypotheses (Fig. 1a).

Geography and depth of each lake explained different categories for the most recent 16-yr window (2012-2018) for each lake. For five groups, a tree with both elevation and latitude explained 66% of the apparent variance. For three groups, a tree with both maximum depth and latitude explained 85% of the apparent variance. Lakes at higher latitudes (> 55 N) were exclusively group iv (Fig. 7). Lakes at higher elevation (> 394 m) and latitude were group iv.3 with the longest ice duration and intermediate duration SD (Fig. 7a). Lakes at lower latitudes but higher elevations tended to be group i, ii (Fig. 7a). Lakes at a lower latitude, between 40 N and 55 N and deeper maximum depth were group i, ii, and iii while shallower maximum depth (< 29 m) were in group iv (Fig. 7b).

Across the 39 lakes, ice duration residuals were significantly related to a range of climate and teleconnection variables. Selected GAMs explained between 8% and 59% of the deviance in ice duration residuals (Fig. 8). Between 0 and 6 explanatory bimonthly variables (median = 2) were significant for each lake ( p < 0.05; Fig. 8). Of all the climate and teleconnection variables, NAO for ON (n = 17) and the global temperature anomaly for AM (n = 16) were the most common significant explanatory variables. In general, higher global temperatures in any bimonthly period resulted in shorter-than-expected ice durations (Fig. 8). Similarly, increasing NAO indices in ON resulted in shorter-than-expected ice durations (Fig. 8).

Not all lakes experienced increasing interannual variability in lake ice duration, despite most experiencing unprecedented rates of recent ice loss supporting our initial hypothesis. Therefore, as for each lake. The violin plot shows the distribution of the points and the lines on the violin plots represent the quartiles for each distribution with a wider spread between lines indicating more variability. Fitted GAM models between mean and standard deviation (SD) for (b) ice-on and (c) ice-off dates for each sequential window size from 4 to 30 yr.

lakes continue to warm and ice duration decreases (Sharma et al. 2021b), we can anticipate an increase in variability until ice seasons last $ 1 month (Fig. 6). After which, there are increasingly high numbers of ice-free years with decreasing variability and, eventually, lakes may cross a tipping point to either have a sequence of ice-free years or become permanently ice-free as forecasted by Sharma et al. (2021a) but will remain to be seen in the coming decades if greenhouse gas emissions are not mitigated. This suggests that year-to-year variability in ice duration will be larger when there is a short duration of ice cover. Geography, air temperatures, and lake depth were found to drive the trends of ice variability, in addition to the frequency of ice-free years (Figs. 4, 7), suggesting that there may be some lakes that are naturally more variable or sensitive to changes in climate than others. Finally, in many lakes, year-to-year variability responded to both large-scale indices of climate change and teleconnections such as NAO and ENSO.

Most lakes have been experiencing a rapid decline in ice duration (Fig. 3), consistent with other lakes and rivers in the Northern Hemisphere (e.g., Magnuson et al. 2000;Newton and Mullan 2021;Sharma et al. 2021b). Several lakes in this study did not have decreasing ice durations because they have already transitioned to predominantly ice-free lakes (e.g., Fig. 2a). On average, lakes were losing 21.7 d of ice per century using the sequential window technique in this study, which was similar to rates calculated using linear regression for these lakes in a prior study (Sharma et al. 2021b). The duration SD gained an average of 4 d per century, with many lakes increasing in variability. This reflects the potentially increasing variability of both components of ice duration, ice-in and ice-out, which is driven by regional weather conditions and the rate of change of those weather conditions at either end of the winter season (Kratz et al. 2000;Arp et al. 2013). Notably, some lakes had decreasing variability, counter to previous studies indicating only increasing or no change in variability (Weyhenmeyer et al. 2011;Benson et al. 2012;Kainz et al. 2017); this phenomenon may be a potential indicator of an ice-free future.

Ice conditions, air temperature, and depth had the largest effects on trends in duration variability. Lakes experiencing icefree winters for more than half of the time experienced rapidly decreasing variability in ice duration, most rapid rates of ice loss, and are vulnerable to permanent ice loss if greenhouse gas concentrations are not mitigated (Sharma et al. 2021a,b). Air temperature is closely linked with ice duration (Palecki and Barry 1986;Robertson et al. 1992;Duguay et al. 2006) and we confirm that this extends to trends in ice variability (Fig. 4). Lakes found in the southern regions of the "slush zone" in the United States and Eurasia where daily winter air temperatures reach a maximum of around or just below 0 C have increasing variability (Fig. 4b) and are most sensitive to the increased frequency of extreme ice-free years (Filazzola et al. 2020). The deepest lakes, which are also vulnerable to short ice duration, intermittent ice cover, and some of the fastest rates of ice cover loss (Sharma et al. 2019(Sharma et al. , 2021b)), are increasing in ice duration variability. Larger and deeper lakes require consistently colder air temperatures because larger volumes of water must be cooled in the late fall and early winter (Brown and Duguay 2010;Arp et al. 2013;Magee and Wu 2017). Large lakes with long fetches are also more sensitive to wind action, breaking the skim of ice at the beginning and end of the ice season (Brown and Duguay 2010;Leppäranta 2010;Magee and Wu 2017). For example, Grand Traverse Bay in Fig. 7. Regression tree results for the most recent 16-yr window (2012-2018) for each lake (a) using 5 groups identified by k-means cluster analysis and using (b) three groups, collapsing all groups from iv.1 to iv.3 down to iv. The lines indicate split points from optimal regression trees for the explanatory variables, including latitude, elevation, and maximum depth (max. depth) for each lake.

Lake Michigan and Bayfield in Lake Superior had the highest variability in ice duration.

Ice phenology exhibited more variability at the beginning of the season than at the end (Fig. 5a), consistent with other lakes (Kratz et al. 2000;Zdorovennov et al. 2013).

Ice-on dates are controlled by local factors like freezing air temperatures, precipitation, and low wind that will set-up ice formation (Duguay et al. 2006;Mishra et al. 2011;Hou et al. 2022). Ice-off dates still are dependent on crossing the 0 C threshold at the end of the ice-season but also reflect the entire winter season with precipitation on ice, ice thickness, and snow cover and drive the timing of ice melt (Jensen Fig. 8. Relationships of ice duration residual and bimonthly average global temperature anomaly (GTA), North Atlantic Oscillation (NAO), and El Niño-Southern Oscillation (ENSO) for October/November (ON), December/January (DJ), February/March (FM), and April/May (AM) as determined by a general additive model (GAM). Any significant parameters were identified by a filled tile, the smooths for each relationship are plotted as a black line to see the direction and shape of the trend. The right panel indicates the percentage of deviance explained (Dev. Expl.) for each lake's GAM fit. et al. 2007;Preston et al. 2016). Despite both ice phenology metrics increasing in variability as the ice season shortens, ice-off dates exhibit more nonlinear patterns. Ice-on dates could continue to increase in variability while ice-off dates exhibit a nonlinear curve that we originally hypothesized that ice duration followed and likely drives more of the ice duration pattern. Ice duration is a better metric for understanding patterns in lake ice variability because ice duration captures ice phenology from both the start and end of the season, while also allowing for incorporation of ice-free years.

Earlier studies had suggested that variability increases with shortened ice duration (i.e., Weyhenmeyer et al. 2011;Sharma et al. 2016), yet we observed a nonlinear relationship between variability and ice duration both across and within lakes over time (Fig. 6). The previously undocumented nonlinear relationship between variability and ice duration may now be apparent because of accelerated rates of ice loss and warmer winter temperatures contributing to a higher occurrence of ice-free years in lakes around the Northern Hemisphere in recent decades (Sharma et al. 2019;Newton and Mullan 2021), a phenomenon which was not as widespread in earlier studies (Weyhenmeyer et al. 2011;Benson et al. 2012). Our new analysis with ice duration (Fig. 6) is more reflective of the current state of Northern Hemisphere lakes as they move from consistent ice cover to intermittent or no ice winters.

The critical transition points from increasing to decreasing variability at $ 1 month may portend ecological regime shifts, as variability changes can be an early-warning indicator of an impending regime shift (Scheffer et al. 2001). Once lakes cross that boundary and begin to have decreasing variability, the shift to ice-free winters may be an inevitable outcome. Within the past 90 yr, some of our study lakes have already transitioned to a new ecological state and represent the endpoints of the mathematical relationship where they are now permanently ice-free and, therefore, have no interannual variability (Fig. 6).

Our initial hypothesis was that there would be four different groups within this mathematical relationship (Fig. 1a). These groups could either represent the characteristics of a lake as a whole or represent intervals of time for a particular lake, which might not be fixed in time as ice duration declines. Because of the sharp decline in the shape of the curve, lakes in groups i and ii were lumped together by the clustering analysis (Fig. 6b) but represent high variability decreasing to completely ice-free. Geography and depth were the best predictors of the groups identified for the most recent 16-yr window (2012-2018), which is consistent with other studies (Arp et al. 2013). Lakes found at higher latitudes were consistently higher in ice duration and had moderate but increasing duration SD. The cutoff for latitudes between 50 N and 62 N is consistent with the 61 N boundary below which lakes are highly susceptible to ice loss (Weyhenmeyer et al. 2011). At the lower latitudes, the deeper lakes at higher elevations were the most likely to be in group i, ii in lakes with these lakes most sensitive to experiencing ice-free years and intermittent ice cover (Sharma et al. 2019).

Although lower elevation sites tend to be less climatically variable (Palazzi et al. 2019), we observed higher variability at low elevations, likely driven by warmer air temperatures and less winter snowpack, causing shorter ice seasons (Palecki and Barry 1986;Brown and Duguay 2010;Arp et al. 2013).

Overarching trends in lake ice decline are ultimately linked to climate change (Magnuson et al. 2000;Sharma et al. 2019). For example, higher global temperature anomalies, especially in AM, result in shorter ice seasons (Fig. 8), likely affecting spring melt for many Northern Hemisphere lakes. However, global temperature and weather patterns vary from year to year, with the effects of climate change on regional and local drivers of limnological processes like lake ice being modulated by teleconnections (Wilkinson et al. 2020). The resulting synergistic or antagonistic between climate change and teleconnections could result in extremes in ice duration; for example, variance in ice phenology has been attributed to NAO or ENSO teleconnections (Sharma and Magnuson 2014;Bai et al. 2012;Schmidt et al. 2019). In this study, many northern European lakes had their ice duration affected by ON NAO where NAO effects are strongest in the early winter (Hurrell et al. 2002). With climate change driving greater variability and extremes in some of these oscillations (e.g., ENSO, Wang et al. 2019), lakes may also experience abrupt shifts in their phenology between years in response to phase switches of teleconnection patterns or especially strong teleconnection years (Bai et al. 2012;Wang et al. 2012). Teleconnections and the global temperature might be better predictors of long-term and ecosystem-wide processes, such as lake ice duration, because they integrate direct drivers, such as meteorology, over space and time (Hallett et al. 2004).

There was a wide range in the deviances of ice duration residuals explained by the global temperature anomaly and the two teleconnection indices that we examined. Depending on the timing of ice-on and ice-off, some lakes may be less responsive to metrics averaged bimonthly. Location may play a large role as well; for example, NAO strongly affects the Atlantic basins of both North America and Europe (Hurrell et al. 2002), but lakes inland from the Atlantic Ocean might not be as responsive. Similarly, different geographic regions might respond to the teleconnections differently, positive NAO indices link to warm conditions in northeastern North America and Northern Europe but cooler conditions in southern Europe (Hallett et al. 2004). Northern European lakes in this study had a negative relationship between ice duration and NAO indices for late fall and early winter months (Fig. 8).

The effects of climate change on ecological, societal, and physical processes have frequently been identified as nonlinear processes (e.g., Grünig et al. 2020). Our results confirm nonlinear responses for ice cover dynamics, with shifting interannual ice phenology variability patterns if lake ice cover lasts for less than a month. The observed shifting patterns in lake ice variability will have consequences for both humans and ecosystems making planning for recreational opportunities, such as skating races and ice fishing tournaments, even more difficult (Magnuson and Lathrop 2014;Knoll et al. 2019). Ultimately, these recreational events will be permanently lost when lakes no longer freeze in warmer winters. The loss of ice cover for lakes can promote summer warming of lakes and harmful cyanobacterial blooms thereby reducing freshwater ecosystem goods and services such as recreational activities and access to potable water (Weyhenmeyer et al. 2008;Hampton et al. 2017). Future studies on the cryosphere should include an analysis of interannual variability to serve as early-warning indicators and identify which systems may be approaching an ice-free state with deleterious effects on freshwater ecosystem goods and services year-round.