Earlier ice breakup induces changepoint responses in duration and variability of spring mixing and summer stratification in dimictic lakes

Reductions in ice cover duration and earlier ice breakup are two of the most prevalent responses to climate warming in lakes in recent decades. In dimictic lakes, the subsequent periods of spring mixing and summer stratification are both likely to change in response to these phenological changes in ice cover. Here, we used a modeling approach to simulate the effect of changes in latitude on long‐term trends in duration of ice cover, spring mixing, and summer stratification by “moving” a well‐studied lake across a range of latitudes in North America (35.2°N to 65.7°N). We found a changepoint relationship between the timing of ice breakup vs. spring mixing duration on 09 May. When ice breakup occurred before 09 May, which routinely occurred at latitudes < 47°N, spring mixing was longer and more variable; when ice breakup occurred after 09 May at latitudes > 47°N, spring mixing averaged 1 day with low variability. In contrast, the duration of summer stratification showed a relatively slower rate of increase when ice breakup occurred before 09 May (< 47°N) compared to a 109% faster rate of increase when ice breakup was after 09 May (> 47°N). Projected earlier ice breakup can result in important nonlinear changes in the relative duration of spring mixing and summer stratification, which can lead to mixing regime shifts that influence the severity of oxygen depletion differentially across latitudes.

As the climate continues to warm, one of the most widespread physical responses of lake ecosystems has been longterm decreases in ice cover (Magnuson et al. 2000;Benson et al. 2012). These decreases in ice cover have vast ecological, economic, and cultural implications (Knoll et al. 2019;Sharma et al. 2019). Both the duration and the timing of ice cover phenology have been assessed for trends over time, and important nonlinear patterns have emerged. For example, in lakes < 61 N, the timing of ice breakup and thereby the duration of ice cover showed high and increasing variability with decreasing latitude compared to more northern lakes (Weyhenmeyer et al. 2011). In lakes at 61 N, the mean timing of ice breakup was on 10 May (day of year 130), and variability in timing of ice breakup increased when ice breakup was earlier and remained relatively constant when ice breakup was later. Continued climate warming may lead to greater variability in timing of ice breakup and total ice cover duration even at latitudes north of this 61 N changepoint due to shorter and more variable periods below the critical 0 C freezing threshold. Currently, lakes in mid-latitude or temperate regions are more likely to experience reduced to intermittent ice cover with even modest warming rates, whereas more northern lakes will require greater warming before experiencing intermittent ice cover Warne et al. 2020). However, northern latitudes are experiencing amplified, rapid rates of climate warming (IPCC 2013), making them vulnerable to very rapid changes in ice cover phenology and the associated ecological consequences if future warming continues.
With changes in timing and duration of ice cover, the phenology of subsequent periods of spring mixing and summer stratification will also necessarily be altered. For example, the duration of spring mixing has been extended as ice cover duration decreased in Lake Langtjern, Norway (Couture et al. 2015). In other studies of individual and multiple lakes, the duration of summer stratification has become notably extended (Austin and Colman 2007;Fang and Stefan 2009;Woolway et al. 2017;Ayala et al. 2019;Woolway and Merchant 2019). As with timing of ice breakup, timing of the onset of summer stratification is strongly influenced by seasonal climatic patterns in air temperature, solar radiation, and wind (Boehrer and Schultze 2008), and can also vary substantially in the same lakes in warm vs. cool years (Read et al. 2014). A space-for-time analysis of 14 lakes with a wide range in lake morphometry from throughout the northern hemisphere found that, across lakes, both spring mixing and summer stratification were extended when ice cover duration was shorter, with both showing nonlinear patterns (Pilla 2021). Here, when ice cover duration was relatively short, typically in mid-latitude lakes, spring mixing duration was extended, but showed no response in lakes with longer ice cover. On the other hand, lakes with longer periods of ice cover duration showed notably extended periods of summer stratification relative to spring mixing. The relative changes in these phenological periods have important implications for many ecosystem processes, and for periods of oxygen depletion vs. replenishment in particular. For example, longer periods of spring mixing can enhance oxygen replenishment throughout the water column (Couture et al. 2015), whereas longer summer stratified periods can result in greater oxygen depletion (Fang and Stefan 2009;Foley et al. 2012;Rösner et al. 2012). Hence, ice cover seems to play a role in mediating the relative responses of periods of spring mixing vs. summer stratification and the consequent responses of water column oxygen dynamics. We focus here on spring mixing and summer stratification as responses to the previous period of ice cover. Though autumn mixing is also likely to be influenced by phenological changes in ice cover, it is not considered here as it is generally more closely linked to the subsequent onset of ice cover.
Latitude and the corresponding seasonal patterns in air temperature and incident solar radiation are key variables for determining the duration and timing of ice cover (Weyhenmeyer et al. 2011) and the subsequent responses of spring mixing and summer stratification across lakes (Pilla 2021). Single-lake studies can control for lake-specific characteristics that influence lake phenology, such as morphometry, but they are not able to assess the range in climate variation over broad latitudinal gradients. On the other hand, space-for-time studies can capture the influence of these geographic variables on lake phenology across a broad latitudinal range (i.e., Weyhenmeyer et al. 2011;Pilla 2021), but the differences in the other lake-specific characteristics that influence lake phenology cannot be systematically controlled for in such studies, especially when sample size is low. In order to control for lake-specific characteristics while also exploring the influence of latitude on lake phenology, we implemented a modeling approach wherein we model a single lake's phenological responses as if it were located across different latitudes. This approach to systematically model the effect of latitude on lake phenology while also controlling for other lake characteristics, including morphometry, may be able to more clearly elucidate potential nonlinear dynamics between the duration and timing of ice cover phenology vs. spring mixing and summer stratification.
We modeled the responses of Lake Giles, a temperate dimictic lake with seasonal ice cover in northeastern Pennsylvania, USA (Williamson 2020; Fig. 1), which has been extensively studied for more than three decades regarding long-term trends in lake browning (Williamson et al. 2015), thermal structure , dissolved oxygen , and zooplankton responses . Using this approach to "move" Lake Giles to different latitudes to assess lake phenological responses and relationships while controlling for other lake-specific characteristics, we address the following research questions: (1) Given its established influence on ice cover duration, how does latitude influence the duration of the phenological periods of spring mixing and summer stratification? (2) As shorter ice cover periods and earlier ice breakup become more common, how do the duration and timing of the subsequent periods of spring mixing and summer stratification respond across a latitudinal gradient?

Meteorological input data
Daily meteorological data were collected from an on-lake weather station from Lake Lacawac, which is 16.6 km from Lake Giles. These data have been used effectively as a surrogate for meteorology at Lake Giles in previous studies (Williamson et al. 2015;Pilla et al. 2018;Pilla and Couture 2021). Daily averages from 1997 through 2020 were used for precipitation (mm d À1 ), wind speed at 2 m (m s À1 ), relative humidity (%), and air pressure (hPa). For daily air temperature and radiation, we used the global Gaussian-gridded data set from National Centers for Environmental Prediction/Department of Energy (NCEP/DOE) 2 Reanalysis data provided by the National Oceanic and Atmospheric Administration/Oceanic and Atmospheric Research/Earth System Research Laboratories Physical Sciences Laboratory (NOAA/OAR/ESRL PSL), Boulder, Colorado, USA, from their website at https://psl. noaa.gov/. We found the coordinate nearest to the location of Lake Giles to represent the local daily air temperature and scenarios across Lake thermal regions (gray circles) and the default location of Lake Giles (gray triangle) that was used for model calibration. The latitudinal gradient spanned multiple lake thermal regions (Maberly et al. 2020) with variability in the occurrence of seasonal ice cover and presence of subsequent periods of spring mixing and stable summer stratification.
radiation from 1997 through 2020. For the model scenarios across a latitudinal gradient, we used the latitudinal coordinates between 35.2 N and 65.7 N all at the same longitude of 97.5 W that spans the central United States and Canada, for a total of 17 modeling scenarios (Fig. 1). This longitude ensures the range of latitudes for the simulations were extracted over land in the center of the continent with less direct influence of coastal weather patterns. Latitudes tested outside of this range either had no occurrence of ice cover (< 35.2 N) or continuous year-round ice cover (> 65.7 N). Besides air temperature and radiation that systematically change over latitude, the other meteorological input variables of precipitation, wind speed, relative humidity, and air pressure were held constant for all modeling scenarios in MyLake using the weather station data from the default location of Lake Giles.

Model calibration and scenarios
We used the open-source MyLake Version 1.2, which is a one-dimensional, process-based model that is renowned for its lake ice module (Saloranta and Andersen 2007). We calibrated seven parameters (Table 1): light extinction coefficient of water for photosynthetically active radiation (PAR) wavelengths (λ = 400-700 nm), light extinction coefficient of water for non-PAR wavelengths (λ < 400 nm or λ > 700 nm), melting snow albedo, melting ice albedo, wind sheltering coefficient, openwater diffusion parameter, and ice-covered diffusion parameter (Saloranta and Andersen 2007). We calibrated the model against the default location of Lake Giles using empirical highfrequency daily-averaged surface (0.5 m) and bottom (22 m) water temperature data measured from Precision Measurement Engineering miniDOT sensors during 11 August 2017 through 05 May 2020. Due to previously known model issues with MyLake capturing ice breakup during the spring of 2018 (Pilla and Couture 2021), we removed data during March 2018 from the calibration routine in order to minimize model bias of temperature data for this month, resulting in a total of n = 1735 daily water temperature measurements used during calibration. The Nelder-Mead simplex calibration routine (Box 1965;Nelder and Mead 1965) from the "nloptr" package in R (Johnson 2018) optimized the root mean squared error of the surface and bottom water temperatures, with a final value of 1.83 C, which is within range of that found for another recent modeling application of MyLake at Lake Giles using highfrequency water temperature data (Pilla and Couture 2021).
The calibrated parameter values (Table 1) were held constant in the subsequent model scenario runs, and used with latitude-specific time series of air temperature and radiation (NCEP/DOE 2 Reanalysis). These modeling scenarios were used to analyze the direct responses to changes in meteorological input data as a function of latitude, as if Lake Giles were replicated and "moved" to a different latitude, with all other lake-specific characteristics remaining constant. We ran the model from 05 August 1997 through 31 May 2020 and used data starting in the winter of 1998 to allow for a 5-month spin-up time. The distributions of air temperature and radiation data for these 17 simulations became less similar to the distribution of these two meteorological variables from the calibration data set as latitude increased (Fig. 2). In general and as expected, northern latitudes had cooler air temperatures and lower radiation compared to the calibration data distribution, but a maximum of 36% and 26% of the air temperature and radiation data in the most northern latitude were not represented in the calibration data range (Fig. 2).
We simulated three additional sets of scenarios to test the influence of lake size and depth on the phenological relationships. Across the same set of latitudes and using the same calibrated parameters and meteorological data as described above, we simulated (1) a small, shallow lake with the same surface area as Lake Giles but with one-third the maximum depth (8 m); (2) a large, deep lake with 10 times the surface area of Lake Giles but with the same maximum depth (24 m); and (3) a large, shallow lake with 10 times the surface area of Lake Giles and one-third the maximum depth (8 m).

Metrics of phenology and statistical analysis
We calculated metrics of phenology from the model output across all years for all scenarios, using the model output variables of ice cover (binary) and vertical temperature profile output ( C at 1-m intervals). Ice freeze date was the first day starting on October 1 or later of each water year that had ice cover present. Ice breakup date was the day following the last day of ice cover that occurred before the end of the water year on September 30. The duration of ice cover was the total number of days with ice cover present between the ice freeze and ice breakup dates. If ice cover did not occur in a given year, it was not included in subsequent statistical models described below. Onset of summer stratification was the first day after ice breakup where a temperature difference of ≥ 1 C per 1 m was reached (Wetzel 2001). The end of summer stratification was the last day before ice freeze when ≥ 1 C per 1 m occurred. The duration of summer stratification was the total number of days when ≥ 1 C per 1 m occurred in the vertical temperature profile output. Finally, the duration of spring mixing was the number of days between ice breakup and onset of summer stratification. We transformed spring mixing duration and summer stratification duration to proportion of days of the year to use a binomial response distribution with their relationships with timing of ice breakup. We compared three candidate models for the relationship between timing of ice breakup vs. spring mixing duration using the "mcp" package in R (Lindeløv 2020): (1) binomial model with no changepoint, (2) segmented changepoint model with binomial distribution on both sides of the changepoint, and (3) segmented changepoint model with binomial distribution before the changepoint and intercept-only after the changepoint. After convergence, models were compared using efficient approximate leave-one-out (LOO) cross-validation, using the "loo" package in R (Vehtari et al. 2020). In brief, this LOO crossvalidation compares theoretical expected log pointwise predictive density (ELPD) across candidate models, and the candidate model with the lowest ELPD is generally considered the best model, analogous to Akaike information criterion, for example. Comparable models are those with a low ELPD difference to standard error difference ratio, where absolute values of this less than $ 1.96 (i.e., 95% probability as with zscores) can be considered competitive models (Lindeløv 2020). We conducted a similar model comparison for the relationship between timing of ice breakup vs. summer stratification duration with three candidate models: (1) binomial model with no changepoint, (2) uninformed segmented changepoint model with binomial distribution on both sides of the changepoint, and (3) informed segmented changepoint model using a uniform prior of the changepoint result from the above model and with binomial distribution on both sides of the changepoint. We used LOO model comparison described above to determine the best model fit for this relationship.
For each of the two best models, we determined the coefficients and residual deviation for relationships before and after  the changepoint. We modeled the relationship and changepoint between timing of ice breakup vs. spring mixing duration for the additional sets of scenarios to assess how differences in lakes morphology influenced this relationship. We assumed the same model structure as the best fit for that of ice breakup vs. spring mixing duration in the primary simulation, but allowed the changepoint to vary with an uninformative prior. All analyses were conducted in R version 4.0.2 (R Core Team 2020), and figures were created using the "ggplot2" (Wickham 2016) and "ggpubr" (Kassambara 2019) packages in R.

Results
As expected, the duration of lake phenological periods varied with latitude ( Fig. 3) with longer and less variable ice cover duration at higher latitudes (Fig. 3a). Spring mixing duration was highly variable at low latitudes, but was typically only a few days at mid-to high-latitudes (Fig. 3b). The duration of summer stratification decreased with increasing latitude (Fig. 2c).
We found support for changepoint relationships between the timing of ice breakup vs. both spring mixing duration and summer stratification duration (Table 2). For the response of spring mixing duration, the one changepoint model with binomial distribution before and intercept-only after the changepoint had the lowest ELPD and standard error. Though the one changepoint model with a binomial distribution both before and after the changepoint was a competing model (difference ratio = 0.1; Table 2), this model had additional parameters and was thereby slightly more complex. From the less complex candidate model of the timing of ice breakup vs. spring mixing duration, the estimated changepoint was 09 May with a credible interval from 17 April to 04 June (Fig. 4a). When ice breakup was before 09 May, spring mixing duration was highly variable and was extended by 2.59% per day earlier ice breakup. When ice breakup was after 09 May, spring mixing showed no change and averaged only 1 day or less, with very low variability. This changepoint represents a 100% decrease in the rate of response of spring mixing duration and an 80% decrease in residual deviance before vs. after 09 May (Table 3).
For the response of summer stratification duration, the best model was that with one changepoint with the prior from the best model for the response of spring mixing duration (uniform prior of 09 May; Table 3). Here, the duration of summer stratification across latitude was extended as ice breakup became earlier, but at different rates around the 09 May changepoint (Fig. 4b). When ice breakup occurred before 09 May, summer stratification was extended by 1.08% per day earlier ice breakup, and when ice breakup occurred after 09 May, summer stratification was extended by 2.25% per day earlier ice breakup. This changepoint represents a 109% increase in the rate of response of summer stratification duration and a 15% increase in residual deviance after the changepoint compared to before (Table 3).
Based on the changepoint of 09 May, two groupings of latitude were established to predict the response in spring mixing duration and summer stratification duration with earlier timing of ice breakup ( Fig. 4c; Table 4). Per each day of earlier ice breakup, spring mixing duration will be an average of 2.59% longer and more variable in lakes at lower latitudes (< 47 N) where ice breakup tends to occur before 09 May, but will not change in lakes at higher latitudes (> 47 N). Summer stratification duration will become longer across all latitudes, but with a slow average rate of response of 1.08% per day of earlier ice breakup in low latitudes (< 47 N) and compared to a faster rate of response of 2.25% longer in higher latitudes (> 47 N; Table 4).
In comparing modeled lakes with differing morphology, we found generally similar relationships between timing of ice breakup vs. spring mixing duration (Fig. 5). All relationships indicated earlier ice breakup was associated with extended and more variable spring mixing duration. Large lakes showed the strongest increase in spring mixing duration with earlier ice breakup, and had up to three times longer average spring mixing periods compared to small lakes. Small, shallow lakes had the lowest average and lowest variability in spring mixing duration of all four lake types, while large, shallow lakes had the highest average and variability in spring mixing duration. Both small and large shallow lakes had similar changepoints to 09 May (22 May and 05 May, respectively), while large, deep lakes had a notably earlier estimated changepoint of 21 March (Fig. 5).

Discussion
Modeling of phenological periods in small, deep Lake Giles across latitudes showed that latitude was an important driver of ice cover duration, timing of ice breakup, spring mixing duration, and summer stratification duration, with important nonlinear changepoint relationships that contrasted for different phenological periods. Earlier timing of ice breakup led to longer duration of both spring mixing and summer stratification across latitudes, but at differential rates based on a key changepoint of ice breakup at 09 May, consistent with that of 10 May reported in Weyhenmeyer et al. (2011). This date represents a time frame following the vernal equinox when there are rapidly increasing rates of change in incoming radiation and air temperature with fewer days below 0 C. These two meteorological variables are important drivers of ice dynamics and lake thermal structure, and they vary systematically with latitude. This changepoint also corresponded with 47 N latitude. Here, earlier ice breakup before 09 May that typically occurred in lakes < 47 N resulted in longer and more variable spring mixing and somewhat longer summer stratified periods. Before 09 May, lower incoming radiation and shorter day length are less likely to allow for rapid formation of thermal stratification immediately following ice breakup, thus the spring mixing period is extended. In contrast, ice breakup after 09 May in lakes > 47 N resulted in fairly consistent spring mixing of only 1 day or less, but with rapid increases in the duration of summer stratification due to rapid onset of thermal stratification following ice breakup. Hence, higher latitude lakes, which are already most likely to experience rapid changes in ice cover phenology , could experience more severe ecological consequences due to much longer periods of summer stratification with minimal change in spring mixing.
A similar changepoint (10 May) has been reported between the average timing of ice breakup vs. its variability, which corresponded to a latitudinal break of 61 N (Weyhenmeyer et al. 2011). Just as the variability in timing of ice breakup increased in lakes when ice breakup occurred before 10 May (Weyhenmeyer et al. 2011), we found here that the variability in spring mixing also increased when ice breakup occurred before our similar changepoint of 09 May. Though a similar changepoint timing, the latitudinal break of 61 N that corresponded with Weyhenmeyer et al.'s (2011) 10 May changepoint in timing of ice breakup was notably higher latitude compared to the latitudinal break of 47 N that corresponded to our 09 May changepoint in timing of ice breakup. This difference in latitude is likely driven by the differences in water bodies included in the respective studies. Weyhenmeyer et al. (2011) included 1449 lakes and rivers with a wide range in morphometry, whereas we controlled for morphometry by focusing solely on modeling small and deep Lake Giles. Differences in lake morphometry alone influence ice and mixing dynamics, where large lakes with larger fetch can have earlier ice breakup and extended spring mixing compared to smaller, more sheltered lakes (Kirillin et al. 2012). We found that when we applied the same modeling scenarios to different lake morphotypes, large lakes indeed had higher average and more variable duration of spring mixing. Small, deep lakes and large, shallow lakes had a similar changepoint timing and credible intervals, suggesting a coherence in lakes of similar volume. In this case, both of these lake morphotypes need to heat a similar volume of water for ice breakup and to induce onset of summer stratification (Brown and Duguay 2010;Kirillin et al. 2012;Warne et al. 2020). Small, shallow lakes had a less precise changepoint, as they are most likely to have shorter periods of spring mixing due to limited fetch and thereby minimal wind-driven mixing. Hence, a precise changepoint in the timing of ice breakup for small, shallow lakes may be less clear given the overall low variability in spring mixing duration across the range in timing of ice breakup. Large, deep lakes had a notably earlier changepoint for timing of ice breakup compared to the other three morphotypes. This may indicate the importance of interannual differences in wind mixing that can strongly influence mixing duration in lakes with a large fetch, which would not be directly related to timing of ice breakup. In large, deep lakes, intermittent spring mixing is also common depending on winter and spring weather conditions, which can have long-lasting influences on deep-water oxygen dynamics (Flaim et al. 2020). These links between phenological periods of ice cover and subsequent spring mixing show important nonlinear relationships with generally consistent changepoints across studies, but with important influences of lake morphometry.
Both the duration and the variability in spring mixing varied with latitude, which influence the efficacy of vertical water column mixing. For example, in many small lakes including Lake Giles, a spring mixing period of just a few days was sufficient for full vertical mixing, but other lakes can experience incomplete vertical mixing when the spring mixing period is less than 3 days (Pilla 2021). In our study, small, deep lakes like Lake Giles at higher latitudes > 47 N averaged just 1 day of spring mixing. At these higher latitudes, ice breakup is more likely to occur after 09 May, when solar radiation and air temperatures are higher and can promote the rapid onset of thermal stratification and thereby result in a shorter period of spring mixing. Because of their short average spring mixing duration, high latitude lakes are more likely to experience incomplete vertical mixing before the rapid onset of summer stratification. On the other hand, lakes at lower latitudes < 47 N had an average duration of spring mixing of 4.9 d with a standard deviation of 7.4 d. In small, sheltered lakes, much of the interannual variability in spring mixing duration and thereby the efficacy of full water column mixing can be linked to weather patterns immediately following ice breakup (Boehrer and Schultze 2008). For example, cool air temperatures and high winds following ice breakup can extend spring mixing and promote full water column mixing (Adrian et al. 2009;Foley et al. 2012). In contrast, warm air temperatures and low wind speeds following ice breakup promote the rapid onset of thermal stratification, thereby limiting the spring mixing period (Adrian et al. 2009;Foley et al. 2012). Hence, in years with warm air temperatures and low wind speeds following ice breakup, even if it occurs before 09 May, onset of stratification may be rapid and lead to short spring mixing periods and incomplete water column mixing. The expected lengthening of the summer stratified period resulting from earlier ice breakup across all latitudes, though at different rates, is likely to result in more severe deep-water oxygen depletion. The relationship between longer summer stratification duration and increased oxygen depletion has been reported in several studies (Fang and Stefan 2009;Foley et al. 2012;Rösner et al. 2012). The consequences of this response may be especially pronounced at higher latitudes due to the strong increase in summer stratification duration paired with minimal change in the spring mixing. Here, earlier ice breakup would result in minimal change in spring mixing duration and thereby no extension of vertical mixing for oxygen replenishment, which could lead to incomplete vertical mixing of oxygen to deeper waters. This potential for incomplete oxygen replenishment can carry over to the longer periods of summer stratification that are tightly linked with worsened oxygen depletion and increased prevalence of hypoxia or anoxia (Fang and Stefan 2009;Foley et al. 2012;Rösner et al. 2012). The ecological consequences of longer hypoxia and anoxia are vast, ranging from loss of oxy-thermal habitat for fishes (Fang et al. 2004;Guzzo and Blanchfield 2017) to increased anaerobic production of greenhouse gasses (Marotta et al. 2014).
While small dimictic lakes like Lake Giles are globally dominant (Downing et al. 2006), the generalizability to lake with other morphologies or mixing types requires further investigation. For example, ice breakup in Lake Langtjern (Norway, 60.372 N) had become earlier over the past four decades, but rarely occurred before 09 May (Couture et al. 2015). Based on our results, ice breakup at or after 09 May would be expected to result in very short spring mixing. However, the primary  Fig. 4a), as well as three simulated lakes of differing morphology: large and deep, small and shallow, and large and shallow. The significant changepoint for each lake type for spring mixing duration vs. timing of ice breakup is marked with a vertical dashed line, and gray band indicates the 95% credible interval. Note that y-axes are listed in days, but relevant statistical analyses were conducted with the response variables as proportion of the year. response to earlier ice breakup in this lake was extended periods of spring mixing, with important implications for deep-water oxygen ventilation in spring (Couture et al. 2015). Lake Langtjern is slightly smaller in both surface area and maximum depth compared to Lake Giles, but with two inlet streams and an elongated basin shape and long fetch. Hence, lake-specific differences in hydrology and basin shape that play an important role in water column mixing can modify the expected response of spring mixing duration to earlier ice cover. Large lakes in particular may have more variable responses due to the two to three times larger range in spring mixing duration compared to small lakes. For instance, changes in ice cover phenology have been linked to changes in mixing patterns in large Lake Tovel (Italy,46.261 N), where preceding periods of autumn mixing seem to be a more important driver of deep-water oxygen than spring mixing (Flaim et al. 2020). Further analysis of the hydrological and watershed controls of the nonlinear responses of year-round lake phenology would be a valuable addition to our understanding and prediction of changing lake ice cover and mixing regimes.
Changes in ice cover and mixing dynamics can lead to regime shifts in lakes, especially those at lower latitudes < 47 N. This latitude of 47 N is just above the boundary between northern temperate and northern warm thermal regions in the western hemisphere ( Fig. 1; Maberly et al. 2020). These two thermal regions have a distinct difference in lake ice coverage, with 92% of higher latitude northern temperate lakes experiencing ice cover, compared to only 9% of lower latitude northern warm lakes (Maberly et al. 2020). Many lakes in North America below 47 N already experience intermittent ice cover . With just moderate climate warming, lakes up to 45 N to 50 N will begin experiencing intermittent or total loss of ice cover. Furthermore, many dimictic lakes in temperate latitudes may experience alteration of mixing regimes due to climate warming and changes in ice cover (Woolway and Merchant 2019). Up to 80% of lakes are likely to experience a mixing regime shift by the end of the century, including 30% transitioning from dimictic systems to either monomictic or polymictic (Woolway and Merchant 2019). Combined, lakes < 47 N are most susceptible to total loss of ice cover and resulting alteration of mixing regimes with continued climate warming. These potential regime shifts from continuous to intermittent ice cover or from intermittent to no ice cover both fundamentally alter the lake ecosystem structure, including changing the phenology of underwater light regimes as well as frequency of vertical mixing of oxygen and nutrients that can in turn influence lake productivity (O'Reilly et al. 2003) and habitat availability (Fang et al. 2004;Guzzo and Blanchfield 2017).
The responses of spring mixing and summer stratified periods to earlier ice breakup can lead to a range of biological responses and ecological consequences. Longer summer stratified periods have amplified trends of warming surface waters (Austin and Colman 2007) and led to greater deepwater oxygen depletion (Fang and Stefan 2009;Foley et al. 2012;Rösner et al. 2012), resulting in reductions in oxythermal habitat for fishes (Fang et al. 2004;Guzzo and Blanchfield 2017). Earlier onset of summer stratification can result in mismatches in predator-prey interactions as each responds to changing phenology at different rates. For example, earlier onset of stratification in Lake Washington (Washington, USA) resulted in an approximately equal advancement in the spring phytoplankton bloom, but no change in zooplankton blooms of Leptodiaptomus and Daphnia (Winder and Schindler 2004a). This decoupling of food resources and zooplankton grazers can result in critical imbalances in lake food webs, as Daphnia in particular are important links to higher trophic levels (Winder and Schindler 2004b). Higher trophic levels may experience substantial trophic mismatches due to temperature-dependent spawning (Straile et al. 2015) or seasonal shifts in feeding behaviors (Wagner et al. 2012). More resilient and adaptable species will have a lower risk for trophic mismatch or resulting population decline, and highly tolerant species a lower risk with changes in oxy-thermal habitat due to changing lake phenology in response to climate warming and reduced ice cover. These ecological implications are important to consider as lakes are already experiencing rapid climate warming (IPCC 2013) and reductions in ice cover ) that drive phenological responses and potential regime shifts. The latitude and morphometry dependence of these differential responses of lake phenology therefore have critical implications for ecosystem structure and function, such as severe deep-water oxygen depletion, trophic mismatches, and enhanced greenhouse gas production.