Loss of Ice Cover, Shifting Phenology, and More Extreme Events in Northern Hemisphere Lakes

Long‐term lake ice phenological records from around the Northern Hemisphere provide unique sensitive indicators of climatic variations, even prior to the existence of physical meteorological measurement stations. Here, we updated ice phenology records for 60 lakes with time‐series ranging from 107–204 years to provide the first re‐assessment of Northern Hemispheric ice trends since 2004 by adding 15 additional years of ice phenology records and 40 lakes to our study. We found that, on average, ice‐on was 11.0 days later, ice‐off was 6.8 days earlier, and ice duration was 17.0 days shorter per century over the entire record for each lake. Trends in ice‐on and ice duration were six times faster in the last 25‐year period (1992–2016) than previous quarter centuries. More extreme events in recent decades, including late ice‐on, early ice‐off, shorter periods of ice cover, or no ice cover at all, contribute to the increasing rate of lake ice loss. Reductions in greenhouse gas emissions could limit increases in air temperature and abate losses in lake ice cover that would subsequently limit ecological, cultural, and socioeconomic consequences, such as increased evaporation rates, warmer water temperatures, degraded water quality, and the formation of toxic algal blooms.

variations, with changes in the timing of seasonal ice formation and loss (collectively referred to as ice phenology) responding sensitively to changes in the climate, notably air temperature (Marszelewski & Skowron, 2006;Palecki & Barry, 1986;Ruosteenoja, 1986). In turn, lake ice phenology is considered an important sentinel of climate change , as well as essential climate variables (Woolway et al., 2020). Furthermore, the dynamics of lake ice cover, as well as its temporal variation, is known to drive important physical, chemical and biological events in lakes, including cross-seasonal cascades that involve both biotic and abiotic processes (Blenckner et al., 2007;O'Reilly et al., 2015;Straile et al., 2003;W. Wang et al., 2018). A detailed understanding of the timing of lake ice responses to climatic variations is, therefore, essential for climate change impact studies.
In a seminal paper on lake ice phenology, Magnuson et al. (2000) observed long-term changes in lake ice cover in 20 spatially and morphologically heterogeneous lakes from 1855 to 1995. These long-term records of lake ice phenology demonstrated that many Northern Hemisphere lakes experienced earlier ice breakup, later ice freeze-up, and shorter ice duration within a warming world. In a follow-up study, Benson et al. (2012) updated these Northern Hemisphere lake ice records by 10 years to 2004 and observed that, over 150 years, and the most recent 100 and 30-year periods, lake ice-on was further delayed, ice-off was earlier, and ice duration was considerably shorter than previously calculated. Benson et al. (2012) also expanded on the previous analysis of Magnuson et al. (2000) by investigating the increasing occurrence of no ice years, as well as late ice formation and/or early ice break-up, which can be affected by changes in both the mean timing as well as the variability of the annual formation and break-up of lake ice. Subsequent studies of lake ice cover have demonstrated the effects of lake morphometry and location on ice phenology; for example, deeper lakes with longer fetches are most sensitive to climatic change resulting in the increased intermittency of ice cover (Magee & Wu, 2017;Sharma et al., 2019Sharma et al., , 2021Woolway et al., 2020).
In the present study, we extended the database  used by Magnuson et al. (2000) that was up to 1995, and thereafter by Benson et al. (2012) that was up to 2004, by another 15 years up to 2019 and for an additional 40 Northern Hemisphere lakes (i.e., 60 lakes in total). We investigated if lakes were now losing their ice cover faster than previously anticipated, and thus continuing on the warming trajectory previously reported (Benson et al., 2012;Magnuson et al., 2000). Considering the continued increase in air temperature as well as the increased occurrence of extreme climatic events in recent decades (Sanches-Lugo et al., 2020;Seneviratne et al., 2012), we hypothesized that lakes were now losing their ice at a rate faster than previously observed. As air temperatures reach unprecedented territory, updating our knowledge of changes in lake ice cover is of critical importance. In this study, we asked: (a) What are the updated trends in ice phenology in the face of continued climate change? (b) How do the rates of warming differ between 25-year periods from 1917-2016? and (c) Are lakes experiencing more extreme ice phenology including later ice-on, earlier ice-off, shorter ice duration, and ice-free years in recent decades than previously observed?

Data Acquisition
We assembled lake ice phenology records with time-series extending over 100 years for 60 Northern Hemisphere lakes ( Figure 1); 20 of these same lakes were analyzed by Magnuson et al. (2000) and 43 of the lakes analyzed by Benson et al. (2012). Our study lakes were located between 36°N to 66°N latitudes covering the entire range of dimictic and cold polymictic lakes (Lewis, 1983). The length of the time series used ranged from 107 to 204 years, with a mean record length of 144 years. Our team updated the lake ice phenology records stored at the National Snow and Ice Data Center  to the 2018-2019 winter season, where possible (Figures 1 and 2). Forty lakes are in North America, 18 lakes in Europe, and 2 lakes in Asia. Most datasets are representative of entire lakes; however, some (e.g., Bayfield Bay in Lake Superior) report ice conditions at a specific location and are not representative of the larger lake. For ease of discussion, we refer to all study sites as "lakes" in this study. Only lakes with less than 16% of missing years of data were included in the analysis, although ∼75% of lakes had less than 5% of missing years and 90% of lakes had less than 10% of missing years. For example, 38 of the 60 lakes had only two or less missing years. For lakes with missing data, missing years sometimes coincided with wars and were, very rarely, in the most recent decades. In addition, we obtained information on lake geography, morphology, and climate variables for each study lake. Latitude, longitude, elevation, surface area, mean depth, and maximum depth were obtained from the National Snow and Ice Data Center . Missing geography and morphology data were supplemented from the published literature or from HydroLAKES (Messager et al., 2016; Table S1). Mean annual air temperatures, precipitation, and cloud cover from 1970 to 2010 were acquired for each lake from the University of East Anglia's Climatic Research Unit, available as interpolated measurements from meteorological station measurements onto 0.5° latitude/longitude grids (CRU TS4.02; Harris et al., 2020).

Data Analysis-Changes Over Time
We calculated temporal trends in ice-on, ice-off, and ice duration dates over the entire time-series using linear regressions for each lake. Linear regression slopes were calculated to provide a direct comparison to earlier studies of trends in ice phenology (Benson et al., 2012;Magnuson et al., 2000). We calculated linear regression slopes using the "lm" function in R (R Development Core Team, 2020). We used the mean iceon and ice-off date for years with missing data because the mean of the entire time series proved the best fit in filling in the gap years . Further, the latest ice-on and earliest ice-off dates were used for years when a lake did not freeze similar to the methodology used in Magnuson et al. (2000) and Benson et al. (2012). Additionally, we compared linear regression slopes with non-parametric Sen's slopes and found no significant difference between the trends (Kolmogorov-Smirnov Test, p-value >0.58 in all cases) ( Figure S1).
Next, we calculated linear regression trends for all lakes within common 25-year windows to quantify the rates of warming among lakes and time periods. We calculated linear slopes for ice-on, ice-off, and ice duration for the 25-year time periods (1917-1941; 1942-1966; 1967-1991; 1992-2016) for which we had the most consistent number of lakes with ice phenology records. To assess which trends from the four periods were significantly different from one another, we performed a Kruskal-Wallis test on the linear regression slopes on the four 25-year time periods, followed by a pairwise Wilcoxon rank sum test to identify which pairwise periods differed significantly.

Extreme Events
We characterized extreme events in lake ice phenology in three ways: (b) identifying the number of ice-free years, (b) recording the number of multiple freeze-thaw events per season, and (c) quantifying the occurrence of extreme ice seasonality events (i.e., ice-on, ice-off, or ice duration) from 2005 to 2019 that were outside the 5th and 95th percentile from all years in the ice phenology records preceding 2005.

Ice-Free Years
To determine if the number of ice-free years has changed through time, we computed the fraction of the lakes that were ice-free for each year in which at least 53 lakes (88% of all lakes in data set) had available ice cover data. We used a local regression smoothing (LOESS with a span = 0.7) on the percentage of lakes with no ice for each year. We then applied a logit transformation to the smoothed data because the many years with no ice-free lakes would have resulted in non-linearity for a logit transformation (Fowlkes, 1987). We used a linear regression to determine the relationship among years and the logit transformed data to identify if the probability of ice-free seasons increased over time. Ultimately, the probability (p) of having an ice-free lake each year is: where b is the slope and a is the intercept fit by the logistic regression, respectively.

Freeze-Thaw Events
Several of the lakes in this data set had multiple freeze-thaw events that represent the thawing of lake ice after it had been established. Multiple freeze-thaw events were recorded by the data providers as multiple and paired ice-on and ice-off dates throughout the winter. We recorded the number of lakes that experienced multiple freeze-thaw events after the first recorded ice-on and prior to the last ice-off within a winter season. We note that multiple freeze-thaw events may affect the calculation of ice duration and highlight the importance of using the same definition consistently throughout the time-series of a lake, as was done here.

Extreme Years
To evaluate if more extreme ice phenology has occurred in recent years, we used the start of the record until 1995 as the baseline data for each of the three response variables, ice-on, ice-off, and ice duration, using the original phenology as day of year and duration in length of days for each lake. We identified extreme events after 1995 as those occurring outside the 5th and 95th percentiles (McPhillips et al., 2018;Seneviratne et al., 2012) that lead to shorter ice seasons. For ice-on day, values greater than the 95th percentile were considered extreme values indicating extremely late ice-on in recent years. For ice-off and ice duration, values less than the 5 th percentile were considered extreme values; these extremes indicate that the ice-off occurred extremely early and ice duration was extremely short in recent years.
For each lake, we required a minimum of 80% of the years for the analysis. On average, lakes had 92.6% of years for ice-on, 94.9% of years for ice-off, and 94.4% of years for ice duration. If the lakes had less than 90% of values, we examined the distribution of missing years to ensure either random distribution or earlier clusters in the time series that would minimize bias in the analysis. If a lake did not freeze during a winter season, then ice-on was set as the maximum value recorded of ice-on day of year over the record plus 1 day, ice-off was set as the minimum value recorded of ice-off day of year over the entire range minus 1 day, and ice duration was 0 days. For each response variable in each lake, we calculated the proportion of years after 1995 that were extreme years. To test whether there are more extreme events recently, we performed a randomized permutation test (Adams & Anthony, 1996) for each response variable. We resampled the entire data set choosing 24 random years and calculated the proportion of years with extreme values 10,000 times to create a probability distribution. We calculated a probability (p-value) of our most recent period  occurring under the probability distribution by integrating from the probability of the sample under that probability distribution out to the most extreme tail. If the p-value was less than our a priori selected α = 0.05, we considered this result significant and that there were more than expected extreme events occurring in the most recent period.
To explore the relationship between ice phenology extremes and geomorphometric and climatic factors, we performed linear regression analyses between the median ice-on date, ice-off date, or duration and the percentage of recent (1995-2019) years that were extreme. Additionally, we created a regression tree for the same response variables but also included variables such as, latitude, longitude, elevation, mean annual air temperature, surface area, cloud cover, precipitation, distance to the nearest oceanic coast, and mean relevant monthly air temperatures. For ice-on, we included mean September, October, and November air temperatures, for ice-off, we included mean March and April air temperatures. For duration, we included all of the preceding months. Surface area was log transformed because of a right-skewed distribution. The most 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 developed in R 4.0.3 using the rpart and rpart.plot functions (Milborrow, 2019;Therneau & Atkinson, 2019).

Changes Over Time in Lake Ice Cover Dynamics
Ice phenology has changed rapidly in the 60 Northern Hemisphere study lakes over the entire time series for each lake (Figure 2). For the 31 lakes with ice-on data, ice-on was on average 11.0 days later per century (SD = 4.8 days/century) over the entire period of record for each lake. Ice-on was delayed by 6-26 days/ century, with significant (p < 0.05) changes for 30 of the 31 lakes (Table S1), suggesting a similarity in the directionality of trends between lakes to delayed ice formation. Ice formation on Lake Suwa was delayed the most by 25.6 days per century since 1897.
For the 58 lakes with ice-off data, ice-off was on average 6.8 days earlier per century (SD = 4.0 days/century, n = 58 lakes) over the entire period of record for each lake (Figure 2). The ice-off trend varied between 3 days/century later and 18 days/century earlier, although 56 of the 58 lakes had trends in the direction of earlier ice-off and 46 of the trends in ice-off were significant. Grand Traverse Bay in Lake Michigan experienced one of the fastest trends in ice-off, melting 16.4 days earlier per century.
For the 30 lakes with ice duration data, ice duration was on average 17.0 days shorter per century (SD = 6.2 days/century) across our Northern Hemisphere study lakes over each lake's entire time series (Table S2). Ice duration shortened by 3-38 days/century, with significant trends for 29 of 30 lakes in the direction of shorter ice cover (Table S2). Bayfield Bay in Lake Superior lost ice the fastest of our study lakes and has experienced a 37.8 days/century reduction in the length of its ice season (Table S2).
Trends in ice-on and duration have increased in the last 25-year window beginning in 1992 relative to the trends from the entire time-series and previous 25-year windows (Figure 3). For example, trends in ice-on dates over the latest 25-year period (1992-2016) were 72 days/century, which is over 6.5 times faster than the long-term trend (11 days/century; Figure 3a). However, for ice-off dates, the 1967-1991 25-year window had the largest changes in the past century such that ice was melting on average 45 days earlier per century ( Figure 3b) and over six times faster than the long-term trend (7 days/century). Trends in ice duration were significantly shorter in the last 25-year window (1992-2016) than any other 25-year windows since 1917 (Figure 3c). More specifically, lake ice duration was 106 days shorter per century in the last 25-year period, over 6 times faster than long-term trends in ice duration (17 days/century). In particular, rates of ice loss were faster for lakes in Sweden and Finland than any other region since 1992.
Trends in ice-on date were related to the median ice-on day of the year with lakes having later ice-on dates changing at faster rates (F 1,28 = 22, p < 0.001, Figure S2a). Median ice-off days were not related to ice-off trends ( Figure S2b). Median ice duration was not related to ice duration trends ( Figure S2c).

Ice-Free Years
Of the 60 study lakes, 9 lakes experienced at least one ice-free year. The last year that all 60 lakes froze for a specific winter was 2002. The 6 years with the highest percentage (>7%) of lakes that were ice-free all occurred between 2001 and 2016 ( Figure 4). The percentage of ice-free lakes increased significantly over time (F 1,107 = 1,370, p < 0.001, b = 0.022, a = −47.34) from 0.5% in the early 20th century to 5% by 2016, on average ( Figure 4).

Freeze-Thaw Events
Only eight of our study lakes had years with recorded multiple freezethaw events within a winter. Since 2005, three of these lakes (Minnetonka, Otsego, and Shields) have had one winter with multiple recorded freeze-thaw events. For Detroit Lake, the earliest multiple freeze-thaw event was recorded in 1949, but the remaining 5 events have occurred since 1989. Similarly, for Lake Monona, Wisconsin, six of the 7 recorded multiple freeze-thaw events have occurred since 1995. On the other hand, Lake Mendota, Wisconsin, multiple freeze-thaw events have been recorded in 1936, 1968, 1977, 1978, 1984, 2001, and 2018. For Cazenovia Lake, 25% of years had multiple freeze-thaw events from 1838 to 1889 and 13% years from 1904 to 1957. More recently, 45% of the years from 2000 to 2019 had multiple freeze-thaw events. Similarly, Oneida Lake had 70% of years between 2000 and 2019 with multiple freeze-thaw events. We note the importance of documenting freeze-thaw events and preparing a standard operating procedure for determining how freeze-thaw events affect the calculation of ice-on, ice-off, and ice duration because the frequency of freeze-thaw events is expected to increase in a warming world.

Extreme Events
Of the 31 lakes with ice-on data, 63% of lakes had significantly more extreme ice-on dates since 1995 (n = 19, Figure 5a). For example, from the baseline data, we would expect 5% of years to be extreme, but, for Mirror Lake, 38% of years after 1995 exceeded the 95th percentile. Of the 58 lakes with ice-off data, 67% had  (1917-1941; 1942-1966; 1967-1991; 1992-2016) for all the lakes (filled circles) and compared the slopes among these four windows. The boxes indicate the median ± one quartile; the whiskers extend to the 5th and 95th percentile. Bars with the same letters represent non-significant differences, whereas bars with different letters are significantly different from each other. significantly more extreme ice-off dates since 1995 (n = 39, Figure 5a). For example, 46% of years after 1995 in Bayfield Bay in Lake Superior had extremely early ice-off dates. Of the 30 lakes with data for ice duration, 88% of lakes had more extremely short ice durations after 1995 (n = 21, Figure 5b). For example, Oulujärvi had 42% of recent years with extremely short ice duration.
Lakes in colder regions at lower elevations were the most likely to experience extremely late ice-ons and early ice-offs after 1995. Lakes in warmer regions had the lowest percentage of extreme ice-on dates while lakes in cooler regions had the highest percentage of extreme years with a small difference between lower and higher elevation lakes (Figure 6a; R 2 = 33%). Lakes in cooler regions at higher elevations had the lowest percentage of extreme events, whereas lakes in cooler regions at lower elevations had the highest percentage of extreme events (Figure 5b, R 2 = 9%). Surface area and April air temperature were the two most important variables in explaining the percentage of extreme short ice durations. The largest lakes had the lowest percentage of extremely short ice seasons (Figure 6c; R 2 = 20%).

Discussion
The rates of change in ice phenology calculated here were among the fastest reported for lakes with over 100 years of lake ice observations. Over a 150-year period, between 1846 and 1995, ice-on was 5.8 days later and ice-off was 6.5 days earlier per century for 37 lakes and rivers distributed across the Northern Hemisphere . Adding 24 additional years of ice phenology records almost doubled the rate of lake ice loss, particularly because of rapid changes in the timing of ice formation. In particular, we found that in the last 25 years, the overall trends for ice-on and ice duration were over 6 times faster than earlier periods, revealing substantially higher rates of ice loss in the most recent decades. Most notably, our analysis suggests that in the most recent 25 years where ice observations are available (1995-2019), lake ice phenology has changed at a rate of 72 days/century and −32 days/century in terms of ice-on and iceoff, respectively, and ice duration has changed at a rate of −106 days/century. The slope of our calculated trends are considerably greater than those reported by Benson et al. (2012Benson et al. ( ) from 1975Benson et al. ( to 2004., the most recent 30 years in their analysis): ice-on was changing at a rate of 15.9 days/century, ice-off at a rate of −18.7 days/century, and ice duration changing at a rate of −43.4 days/century. The more rapid changes observed in our updated data set largely conformed with our expectations given the increase in air temperatures reported in recent decades (Sanches-Lugo et al., 2020) and that air temperature is one of the most important climatic drivers of lake ice dynamics, owing to its additive effects on various components of the lake energy budget (Palecki & Barry, 1986;Vavrus et al., 1996;Weyhenmeyer et al., 2011).
Our investigation also suggested a recent increase, exceeding those reported by Benson et al. (2012), in the number of extreme ice-on and ice-off dates, for the majority of lakes within the data set. This is in agreement with the observed increase in the frequency and severity of extreme climatic events (IPCC, 2014;Seneviratne et al., 2014). In addition, we found that the duration of winter ice cover has decreased, particularly since 1995, to the point where some lakes are beginning to have more winters with minimal or even no ice cover. An increase in the number of ice-free years was also reported by Benson et al. (2012) and Filazzola et al. (2020), but we find these to become even more frequent in recent years in response to warmer air temperatures (Filazzola et al., 2020). Ice-free years are forecasted to become more common as winter air temperatures are projected to continue to warm (Filazzola et al., 2020) and up to 5,700 northern lakes may permanently lose ice cover by the end of the century . In particular, our study suggests that the ice phenology observations that were rare in the pre-1995 climate are now much more common in the 1996-2016 period. Our results also suggest that lakes in colder regions, with traditionally longer ice seasons, have experienced more extreme ice phenology in recent decades. This is likely a result of colder regions typically warming at a faster rate than elsewhere in recent years for example, due to polar amplification (Post et al., 2018;Stuecker et al., 2018). Notably, the rapid increase of air temperatures at higher latitudes can increase the probability of extremes in air temperature (e.g., the Siberian heat wave of Figure 5. (a) Ice-on and ice-off phenology and (b) ice duration for all study lakes for each lake's records prior to 1995 (boxplots) with extreme data from recent 1995-2019 values as light blue points (ice-on: circles, ice-off: squares, ice duration: triangles). Dark gray filled boxes represent lakes with significantly more extreme values in recent  years. The boxes represent the median ±1 quartile, the whiskers represent the 5th and 95th percentiles for each lake's record prior to 1995. Lakes are abbreviated by their first four letters for ease of viewing (see Table S1 for lake names) and ranked by median ice-off day of year from the pre-1995 baseline. 9 of 12 2020) and, as a result, extremes in ice phenology (Diffenbaugh et al., 2017). Similarly, ice duration was most reduced in lakes where it was historically longest (i.e., colder regions). Furthermore, smaller lakes with warmer spring temperatures tended to have the most years with extreme ice durations. Conversely, for lakes with already short ice seasons (e.g., warm regions), reductions in ice cover were not anomalous relative to the long-term record. Additionally, larger lakes had less years with short ice durations, likely because of the heat capacity of water and the energetic requirement for heating or cooling such large volumes of water (J. Wang et al., 2012). With climate change projected to continue with up to 4°C warming in some regions by 2050 (Collins et al., 2013), we expect the loss of ice cover to accelerate in the coming decades as some lakes cross thresholds to intermittent ice cover .
Not only are the studied lakes experiencing more extremes in ice phenology, but our study also suggests that the likelihood of multiple freezing and thawing events has increased in recent decades. Specifically, as winter air temperatures increase in many regions, warming to near or above 0°C for a sustained period, our observations show that lake ice cover is increasingly susceptible to becoming an intermittent phenomenon . Previous modeling studies of lake ice cover, for example, in the midwestern U.S., eastern U.S., and central Europe, had already foreseen the increased prevalence of freeze-thaw events within a warmer climate (DeStasio et al., 2016;Livingstone & Adrian, 2009;Robertson et al., 1992). However, we are now at the stage where these projections are becoming a reality and are increasingly evident in our observational records.
As the variability in ice phenology increases and the intermittency of ice cover becomes the norm in many lakes, we hypothesize that this is indicative of a rapid non-linear change in the response to warming (Carpenter et al., 2011), with potential regime shifts in the physical environment of lakes likely to occur from annual continuous ice cover to intermittent ice cover to permanent loss of ice cover. Such regime shifts are projected by the end of the current century; for example, many dimictic lakes will transition to a monomictic mixing regime (Woolway & Merchant, 2019). Further, Lakes Superior, Michigan, and Suwa are among 5,700 lakes predicted to permanently lose ice cover in the 21st century without greenhouse gas mitigation . However, given the observed changes in our data set, these regime shifts to intermittent winter ice cover may occur much sooner in some lakes, possibly within decades (e.g., Sharma et al., 2019). It is thus critical to understand the implications of these rapid rates of ice loss and whether behavioral, ecological, and evolutionary adaptations are possible to offset the alarmingly fast rates of warming. We encourage Figure 6. Regression tree results for (a) ice-on, (b) ice-off, and (c) duration percentage of years recently (1995-2019) that were outside 5th percentile or 95th percentile from the pre-1995 record. The size of the point represents the number of years with extremely late ice-on dates, early ice-off dates, or short ice durations relative to the number of years with data from 1995 to 2019. The dotted lines indicate split points from optimal regression trees. The numbers represent the mean percentage from all lakes in that regression tree split. lake ice observers to record details on freeze-thaw events within a winter as they become more frequent. Such observations will be essential to observe and predict a potential shift to an ice-free winter in the future.

Conclusions and Implications
Changing ice conditions, as we have described in this study, can have a substantial effect on lake ecosystems as well as the many ecosystem services that lakes provide to society (Flaim et al., 2020;Hampton et al., 2017;Knoll et al., 2019;Powers et al., 2017;Sharma et al., 2020). For example, the loss of lake ice can influence socioeconomic and cultural activities, including winter transportation, ice fishing, recreational skating, and religious rituals (Knoll et al., 2019). The most obvious effects of ice loss for lake ecosystem processes arise due to earlier break-ups and shorter durations of ice cover on spring and summer water temperature with implications for a host of temperature dependent ecosystem processes. Notably, a warming during spring can result in an earlier onset of thermal stratification in many lakes and a subsequent advance in the seasonal succession of both phytoplankton and zooplankton during spring and summer (Adrian et al., 1999;Straile, 2000;O'Reilly et al., 2015;Winder & Schindler, 2004;Woolway et al., 2020). In mesotrophic and eutrophic lakes, warmer spring water temperatures can further affect the timing and duration of the clear-water phase, a period of marked reduction in algal biomass in late spring-early summer caused by Daphnia grazing (Matsuzaki et al., 2021). However, lake ice loss in the future could have even greater effects on local economics (Knoll et al., 2019) and ecosystem processes, including those occurring under the ice Powers et al., 2017). In turn, there is an urgent need for research focused on implications of losing lake ice cover, both economically and in terms of lake ecology. Accumulating knowledge of winter ecology can improve our capacity to understand the role of ice cover on lakes, and the people who depend on ice cover, before it is lost.

Data Availability Statement
All data used in this study are publicly available. Lake ice phenology data was originally from the National Snow and Ice Data Center  with updates to the 2019-2020 winter season published there as well. Metadata for all lakes were also accessed via the National Snow and Ice Data Center  with any missing values supplemented from the published literature or from HydroLAKES (Messager et al., 2016). The metadata is available in Table S1. Mean annual air temperatures, precipitation, and cloud cover from 1970 to 2010 were acquired for each lake from the University of East Anglia's Climatic Research Unit, available as interpolated measurements from meteorological station measurements onto 0.5° latitude/longitude grids (CRU TS4.02; Harris et al., 2020).