On thin ice: Linking elevation and long‐term losses of lake ice cover

Despite long‐term analyses of lake ice phenology globally, comparatively little is known about high‐elevation lakes, for which climate shifts are thought to be occurring faster than at lower elevations. Using a 36‐yr dataset (1983–2018) on alpine lakes (> 3000 m ASL) from the Green Lakes Valley, Colorado (GLV), we found that ice‐cover duration decreased by an average of ~ 24 d, due to both earlier ice‐off (9 d) and especially later ice‐on (15 d). Spring ice thickness also decreased by 0.88 cm yr−1. By comparison, ice‐cover duration in the GLV is decreasing ~ 50% faster than for Northern Hemisphere lakes (n = 215), which translates to an increase in open water duration ~ 2.5 times more in the GLV than the average of the Northern Hemisphere. Our analytical comparison demonstrated more rapid trends in this alpine region compared to lakes more broadly, and especially emphasized the influence of elevation on lake ice phenology.

Ongoing climate shifts are associated with changes in lake ice phenology worldwide . Changes in air temperature and precipitation can alter these patterns and lead to long-term changes in ice phenology (Brown and Duguay 2010). Other factors, such as lake geography (Williams and Stefan 2006), can have important mediating effects for local and regional differences in ice-cover duration. For example, Williams and Stefan (2006) showed that higher latitudes and elevations can delay lake ice clearance and advance ice formation, such that lake ice-cover duration increases by~3.7 d per degree north and~0.04 d per meter in elevation. Although climatic change effects are greater for lakes in cold regions (Williamson et al. 2009), there are fewer long-term empirical data on patterns in lake ice phenology comparing these regions to global trends. Multiple lines of evidence suggest that ice phenology is likely to shift more rapidly and to a greater relative extent in colder regions. Processes such as Arctic amplification and snow/ice-albedo feedback (Henneman and Stefan 1999) can lead to faster changes in lakes in cold regions relative to other regions.
Due to their typically remote locations, mountain lakes are also a generally understudied class of lakes (Catalan and Donato-Rondón 2016), and the paucity of long-term empirical data on alpine lake ice-cover duration has limited efforts to compare trends in ice phenology among lakes globally. The range in values for alpine/montane systems are often defined regionally and less applicable broadly. For example, an elevation of 265 m ASL (above sea level; Sharma et al. 2019) or 608 m ASL (Kainz et al. 2017) would not be representative of high-elevation regions in the western United States, where lakes above 2000 m ASL are common (Bahls 1992). Thus, ice phenology trends for many high-elevation lakes remain unknown and are likely underestimated. Ice phenology trends are also sensitive to the timing and length of the sampling period, such that analyses comparing patterns of ice phenology across elevation gradients should be compared using the same time period (Kainz et al. 2017).
To better understand the patterns of ice-phenology change in alpine lake ecosystems and how they differ from lakes across elevational and latitudinal gradients generally, we analyzed temporal trends in ice-on date, ice-off date, ice thickness, total ice-cover duration, and heterogeneity for lakes in the Green Lakes Valley (GLV) of the Southern Rocky Mountains, which represent the longest-running high-elevation study on lake ice globally (51 yr for seven lakes at 3126-3620 m ASL). This study extends the dataset for lakes in the GLV by 3 yr for ice-off, 12 yr for ice-on, and 16 yr for ice thickness relative to any previously published research (Caine 2002;Williams and Stefan 2006;Preston et al. 2016).
In addition, we compared results from the GLV with those from natural lakes in the Northern Hemisphere along gradients in both latitude and elevation, with the aim of using a consistent analytical framework to assess the magnitude and rate of ice cover changes for the period from 1983 to 2013. Our goals were to determine the degree to which phenological shifts in ice cover associated quantitatively with lake elevation and latitude, thereby comparing trends in cold regions to those more broadly.

Green Lakes Valley
Research in the GLV represents the longest-running study (51 yr) of high-elevation lakes (> 3000 m ASL) in North America. Located in north-central Colorado within the Southern Rocky Mountains and just east of the Continental Divide ( Fig. 1), the GLV includes a series of seven montane to alpine lakes between 3126 and 3620 m ASL, with maximum depths ranging 7-16 m, surface areas 0.4-2.2 km 2 , and mean depths 3.7-8.6 m (Caine 2001). The historical mean annual air temperature for the GLV is about −4.1 C (range: −7.7 to −0.5 C), while annual precipitation averages 993 mm (range: 512-1581 mm) (Caine 2001). Since 1968, ice clearance date has been recorded in Green Lake 4 (GL4), while ice clearance and formation date have been recorded throughout the GLV since 1982. Based on weekly to biweekly observations throughout spring and fall, ice formation and clearance dates were classified visually as when each lake was 100% covered or cleared of ice, respectively. Ice thickness has also been measured since 1982 in GL4 (Caine 2018) by inserting a 320-cmlong probe into a drilled ice hole at the deepest part of the lake between 20 and 31 March, when maximum ice thickness typically occurs. These measurements were collected through "black" ice and did not include snow (see Caine 2002). Ice thickness data are available in the Niwot Ridge LTER repository at https://portal.edirepository.org/nis/ mapbrowse?packageid=knb-lter-nwt.199.2. For this study, we included Green Lakes 1-5 (3400-3620 m ASL) and years in each lake for which ice clearance and the previous year's ice formation were recorded. Ice phenology data are available in the Niwot Ridge LTER repository at https://portal. edirepository.org/nis/mapbrowse?packageid=knb-lter-nwt.106.3 (Caine 2018). This generated a dataset from 1983 to 2018. Due to the importance of open water duration for temperate lakes, we (1) analyzed trends in the duration of open water (ice-free days) using the complete GLV dataset and (2) trends in icecover duration (ice-covered days) for comparisons with an available database on Northern Hemisphere lakes (1983-2013) (Benson et al. 2000).

Northern Hemisphere comparisons
To analyze patterns in lake ice cover in the Northern Hemisphere as a function of latitude and elevation, we obtained data from the Global Lake and River Ice Phenology Database housed at the National Snow and Ice Data Center. Data and metadata are available at http://nsidc.org/data/ g01377.html (Benson et al. 2000). These data encompass centuries of ice records and represent the most comprehensive global dataset in lake ice cover. We limited our analysis to Northern Hemisphere's natural lakes that freeze annually, excluding rivers or streams. We further constrained the analysis to the period between 1983 and 2013, which maximized overlap with the GLV dataset while including the most recent update to the global database at time of submission. Further analysis of the entire Northern Hemisphere dataset can be found in Sharma et al. (2019) and Lopez et al. (2019). For each record, we gathered information on the elevation, latitude, year, and ice-cover duration of each lake, predefined in the dataset as the time between ice formation and ice clearance. In total, we obtained 3221 observations representing 215 lakes and an elevational range of 5-1768 m ASL, latitudinal range of 36.2-74.4 N (Fig. 1), and for lakes with reported area and depth, areas ranged 0.005-31,500 km 2 (mean: 445 km 2 , median: 15 km 2 ), and maximum depths 1.8-1741 m (mean: 50 m, median: 27 m).

Statistical analyses
We used linear mixed-effects models (LME) using the lme4 package (Bates et al. 2015) in R 3.6.0 (R Core Team 2019) to determine trends in ice clearance, ice formation, and open water duration in the GLV. The number of days until ice-off, ice-on, or total days of ice coverage or open water was modeled as a Gaussian distributed response variable, with elevation and sampling year as fixed effects. To account for repeated (nonindependent) observations on each lake, we incorporated lake identity as a random intercept term. We also accounted for random slopes by allowing the annual trend in ice duration to vary by lake. Similarly, the Northern Hemisphere dataset was analyzed to determine ice-cover duration trends while accounting for fixed effects of elevation, latitude, and year, again incorporating a random intercept term for lake identity and random slopes for year. Although other measures such as lake size, depth, and volume may be important, these values are less commonly available and were not found to be significant during initial analysis. Because ice thickness was observed in only a single lake (i.e., GL4), we first used a linear regression model to determine the yearly trend; after further consideration, we applied a break-point analysis to account for potential nonlinearity using the "segmented" package (Muggeo 2020) in R. Heterogeneity of ice-cover duration among years was assessed via the coefficient of variation (CV) among 10-yr moving windows for each lake and among the GLV.
For diagnostics, we examined residuals of resulting models for normality and equality of variance (Fox 2008). For the GLV-specific models, we checked for coherence in ice phenology among GLV lakes by determining the correlation in ice duration by year. For all models, we tested for remaining temporal autocorrelation in the residuals using the "acf" function in R (R Core Team 2019) and verified a lack of collinearity among predictors by calculating variance inflation factors, for which we considered values below 5 to indicate suitability (Akinwande et al. 2015).

Green Lakes Valley
Since 1983, the ice-free period in the GLV has increased by an average of 24 d (LME, β = 0.66, SE = 0.12, p < 0.001; Table 1; Fig. 2), with a range of 17-32 d (Table 1). This pattern was driven by both a 9-d average advance in ice-off date (LME, β = −0.25, SE = 0.09, p = 0.006) and especially a 15-d average delay in ice-on date (LME, β = 0.41, SE = 0.07, p < 0.001; Table 1; Fig. 2). Lake elevation negatively affected open water duration (LME, β = −0.23, SE = 0.03, p < 0.001), such that the highest lakes (e.g., GL4 and GL5) were typically ice-free for 41 fewer days than the lower lakes (GL1-3). Although there was no significant interaction between elevation and year (LME, β = 0.002, SE = 0.001, p = 0.20), estimates of the temporal trend for lakes individually highlight variation in magnitude (Table 1). Model coefficients of determination (R 2 ) were 0.62 (marginal) and 0.65 (conditional upon the random effects). Lastly, from the linear analysis, ice thickness in GL4 decreased at a rate of 0.88 cm yr −1 (Fig. 2), declining from an average 127 cm during the first three observations on record to an average 94 cm by the last three observations (26% decrease). However, the break-point analysis detected a break point occurring in 2000, with little change in ice thickness from 2000 to 2019 (Fig. 2).
Over the 36-yr time series, the average ice-cover duration of lakes in the GLV was 246 AE26 d, which ranged from 273 d Table 1. Coefficients (standard error) for all linear mixed effects models of lakes in the Northern Hemisphere and Green Lakes Valley (GLV; lakes 1-5: GL1-5) datasets. The full GLV dataset is analyzed (1983-2018) for lake ice-free duration. A truncated dataset (1983-2013) is also used for comparing ice-cover duration in the GLV to the Northern Hemisphere dataset and includes lake elevation, latitude, and respective interactions, while also accounting for random slopes (year) within the model. All parameters are significant to p < 0.05 unless noted as "NS." between GL4 and GL5 (r = 0.96). Across all lakes, the interannual heterogeneity (CV) in ice-cover duration decreased slightly (~4%) over this time period.  Hemisphere dataset as well as for the Green Lakes Valley, Colorado; GLV from 1983 to 2013, with each lake's latitude and elevation shown for comparison. Note that lakes in the GLV are the highest elevation but nearly the lowest latitude lakes, while generally, the highest latitude lakes of the Northern Hemisphere dataset are among the lowest elevation.
By comparing the patterns for GLV with those of Northern Hemisphere lakes generally, we found that the decreasing trend in ice-cover duration for GLV (20 d) was~50% faster than the Northern Hemisphere average (13 d), with ice-cover duration for lakes in the GLV decreasing faster than 87% of all lakes in the Northern Hemisphere dataset. On average, the GLV is ice-free for 121 d per year compared to an average of 211 d across the Northern Hemisphere. Stated another way, this means the GLV experienced a 20% increase in ice-free days between 1983 and 2013 compared to an 8% increase for the Northern Hemisphere dataset.

Discussion
Ice phenology in the GLV of the Southern Rocky Mountains has shifted considerably over the last 36 yr. On average, the total annual ice-free season has increased by approximately 3 weeks since 1983. Because these small, highelevation systems are frozen for > 225 d per year, this change amounts to a 20% increase in the open-water period. Statistically, this pattern was driven more strongly by shifts in the timing of ice-on rather than ice-off, with ice formation being delayed at nearly twice the rate of earlier ice clearance. By applying the same analytical approach to well-studied lakes in the Northern Hemisphere sampled over a matching time frame, we found that while ice-cover duration is decreasing for lakes broadly, consistent with past analyses (Brown and Duguay 2010), the annual rate of ice cover loss in the highelevation lakes of the GLV (0.66 d yr −1 ) is occurring~50% faster than the Northern Hemisphere average (0.41 d yr −1 ). Although increasing latitude is associated with large decreases in ice-cover duration (Latifovic and Pouliot 2007), our analyses highlight how rates of changing ice cover concurrently vary with elevation. Due to the short growing season inherent to mountain lakes, the proportional increase in ice-free duration is 2.5 times larger in the GLV compared to the Northern Hemisphere average.
The current results are broadly consistent with previous analyses of alpine lake ecosystems. Our estimated trend in ice clearance (−0.25 d yr −1 ) is similar to another study of the GLV over a shorter time period (−0.21 d yr −1 ; Preston et al. 2016). Comparing across the Northern Hemisphere, Lake Lunz (608 m ASL) in Austria demonstrated a decrease in ice-cover duration of 0.35 d yr −1 from 1921 to 2015 (Kainz et al. 2017), whereas a lake in the Tatra Mountains, Poland (1935 m ASL) decreased ice cover at a rate of 0.7 d yr −1 from 1971 to 2007 (Choinski et al. 2010). Further, a study on the Tibetan plateau (mean: 4500 m ASL) found that ice cover decreased in some lakes from 2001 to 2010 by an average 0.65 d yr −1 (Kropacek et al. 2013). With respect to high-latitude lakes, Latifovic and Pouliot (2007) reported that boreal lakes in the Canadian Arctic exhibited both an earlier break-up (0.99 d yr −1 ) and later freeze-up (0.76 d yr −1 ) between 1985 and 2004 (Latifovic and Pouliot 2007). Overall, however, studies of changes in icecover duration for high-elevation systems-which can be variably defined among regions-remain relatively rare, and comparisons among studies can be highly sensitive to the sampling period and duration (Kainz et al. 2017), underscoring the value of the standardized comparisons among Northern Hemisphere lakes conducted here.
Mountain lakes are dynamic systems that exhibit and integrate numerous facets of a changing environment. The decreasing trend in maximum ice thickness in GL4 may be slowing given that Caine (2002) reported maximum ice thickness decreasing at a rate of 2 cm yr −1 (SE: 0.21) compared to our estimate of 0.88 cm yr −1 (SE: 0.17) in GL4. It is possible that permafrost melt, which was detected by a state shift of SO 4 − in 2000, may be causing the detected break point and subsequent weakening in decreasing ice thickness (Caine 2010). It is also possible another state shift could occur once this permafrost buffer expires. Other factors may also be acting upon lakes to produce varying results. Others have shown how lake morphometry, cloud cover, geography, and catchment characteristics (Brown and Duguay 2010), can influence ice phenology trends over large scales; such changes can also be nonlinear (Supplemental Material). In the GLV, it is possible that differing lake volumes, flushing rates, and topographic shading may be influencing differing rates of change. Comparatively, the lakes in the GLV are small and shallow compared to those of the Northern Hemisphere dataset, potentially amplifying ice-loss effects in the GLV. Also, decreased water clarity derived from chlorophyll (Chl) a or dissolved organic matter (DOM) has been linked to increases in lake temperature (Pilla et al. 2018), and this could have implications for ice phenology. In the GLV specifically, Preston et al. (2016) further showed that Chl a concentrations correlate positively with earlier ice out. It is also possible that elevational advancement of terrestrial vegetation (specifically, Salix glauca), which has been detected in the GLV based on aerial imagery between 1972 and 2008 (Bueno de Mesquita et al. 2018), could enhance inputs of DOM into the aquatic ecosystem, leading to decreased clarity. Finally, for lakes that are distributed within highly populated areas, which is common, we highlight the potential for human activity to further alter lake water quality through shifts in nutrients, ion concentrations, and dissolved organic carbon (Howell et al. 2012). Our findings suggest that, relative to other lakes with ice phenology data in the Northern Hemisphere, lake ice-loss is occurring at faster rates in alpine lakes in the Rocky Mountains. Such changes in ice cover have potential implications for within-lake processes and both the quality and quantity of downstream water resources. Mountain lakes are sensitive systems that provide water for half of the world's population (Woodwell 2004). As water resource demands increase under a growing global population (Vorosmarty et al. 2000), the quality and quantity of available water will be at the forefront of societal importance. We emphasize the need to further examine understudied mountain lakes and their role as a water resource globally.