The vulnerability of lakes to climate change along an altitudinal gradient

Studies of future 21st century climate warming in lakes along altitudinal gradients have been partially obscured by local atmospheric phenomena unresolved in climate models. Here we forced the physical lake model Simstrat with locally downscaled climate models under three future scenarios to investigate the impact on 29 Swiss lakes, varying in size along an altitudinal gradient. Results from the worst-case scenario project substantial change at the end of the century in duration of ice-cover at mid to high altitude (−2 to −107 days), stratification duration (winter −17 to −84 days, summer −2 to 73 days), while lower and especially mid altitude (present day mean annual air temperature from 9 °C to 3 °C) dimictic lakes risk shift to monomictic regimes (seven out of the eight lakes). Analysis further indicates that for many lakes shifts in mixing regime can be avoided by adhering to the most stringent scenario. Shifts in the mixing regimes of lakes in a changing climate are projected to vary with altitude, with mid-altitude lakes particularly affected, according to simulations of 29 Swiss lakes with a physical lake model and downscaled climate projections.

L akes are commonly cited as sentinels of climate change 1 . The physical response of lakes to changes in atmospheric forcing are traditionally assessed by quantifying modifications in their thermal structure. These modifications typically include evolution of lake surface and bottom temperatures, the duration of the summer and winter stratification and the duration of ice cover. The mixing regime of lakes 2 is a physical parameter describing the timing and frequency with which lake temperatures fully homogenize on an annual basis. When compared with other physical characteristics, these events with homogenization of the water column exert the largest overall influence on the functioning of lake ecosystems. A change in mixing regime thus profoundly alters a lake by enhancing or preventing vertical fluxes of nutrients and dissolved gases. This in turn can reshape food web dynamics in the lake [3][4][5][6][7][8][9] . Numerous studies have already reported specific effects of climate change on lake temperature [10][11][12] , ice cover [13][14][15][16][17] , stratification and mixing regime [18][19][20][21] and nonlinear seasonal interactions 22 at global and regional scales. These changes have been linked to trends in atmospheric climate variables (e.g. air temperature) and to individual lake characteristics (e.g. volume, area, transparency) either by statistical methods 12,16,20,22,23 or process-based numerical models 10,11,14,15,19,21,24,25 .
In addition to the well-studied latitudinal gradients in lake response to climatic change 19 , altitude-gradients in atmospheric conditions 26,27 may also influence the thermal structure of lakes 28 . Lake mixing regimes have classically been defined at global scales as a function of latitude and altitude 2 . Given the altitude-dependence of both climate trends 26,27 and lake mixing regimes 2 , we hypothesize that lake response to climate change varies as a function of altitude. This hypothesis carries major implications for the vulnerability of lakes in terms of their ecosystems and their role in the carbon cycle. The remoteness of high altitude lakes however limits the availability of long-term datasets that might elucidate these relations 26,29 . Given the lack of longterm data, modelling approaches can fill in gaps by estimating responses from past lake data and future projected climate scenarios. This study investigates altitude-dependencies on climate change by modelling the response of 29 lakes in Switzerland, spanning an altitudinal gradient from 193 to 1797 m a.s.l. (Fig. 1d and Table 1). The numerical simulations are performed with the deterministic one-dimensional physical lake model Simstrat 30,31 .
An important challenge when modelling lakes along an altitudinal gradient is that the altitude-dependence of climate forcing in Regional Climate Model (RCM) projections is represented by averages over individual grid cells. The complex topography of mountainous areas can strongly affect local atmospheric conditions. In the present study, we address this challenge by using recently developed RCM projections downscaled to the local scale for Switzerland 32 . We specifically focus on the period from 1981 to 2099 using RCM projections for the three Representative Concentration Pathways (RCP), RCP8.5, RCP4.5 and RCP2.6 ( Supplementary Fig. 1). The numbers at the end of the RCPs indicate the increased radiative forcing in W m −2 in the year 2100 compared to preindustrial levels. The worst-case scenario RCP8.5 implies continuously rising global greenhouse gas emissions, in RCP 4.5 emissions are peaking around 2050 and subsequently declining, while the most stringent scenario RCP2.6 was designed to limit global warming to 2°C in line with the Paris agreement and requires stringent measures for emission reduction and net negative emissions towards the end of the twenty-first century 33 . Differences in climate between the late twenty-first century and the base period are here expressed as the median, mean and variance for 2071-2099 compared to 1982-2010.
The results of our simulations indicate that climate change will cause substantial alterations in the thermal structure of the 29 investigated lakes. Mean annual lake surface temperatures are projected to consistently increase for all lakes. Projected changes in other thermal properties, such as lake bottom temperature and the duration of summer and winter stratification as well as ice cover, depend on individual lake properties such as lake volume and altitude. The projected changes in stratification and ice cover duration are largest for high altitude lakes. Nonetheless, these lakes will maintain their annual ice cover and a dimictic stratification regime throughout the twenty-first century. Low to mid altitude (<1500 m a.s.l.) small lakes (<0.5 km 3 ) were found to be especially sensitive to changes. These lakes are projected to lose ice cover and change from a dimictic to a monomictic regime during the twenty-first century. Results based on different emission scenarios indicate that changes in mixing regime and loss of ice cover can be counteracted, but not entirely avoided, with climate protection measures as projected by scenario RCP2.6. The results of the present study cannot be directly transferred to other regions because the altitudinal variation of both lake thermal properties and trends in climate variables are not necessarily similar elsewhere. However, we would expect similar sensitivities of lakes to climate change in other regions where lake mixing regimes range from primarily monomictic at lower altitudes to dimictic at high altitudes. The investigations in the present study are limited to the direct impacts of projected changes in climate conditions on the thermal structure of lakes. Further investigations are needed to quantify possible indirect impacts, especially those resulting from changes in the lake catchments 34 , including but not limited to glacier retreat and permafrost thaw 35 , alterations of snow melt dynamics 36,37 , river flow regimes and resulting changes in loads of suspended particles and nutrients 38,39 . On a hopeful note, this study shows that many lakes can potentially avoid shifts in mixing regimes and sustain ice cover throughout the twenty-first century if greenhouse gas concentration trajectories remain within envelopes envisioned under RCP2.6.

Results and discussion
Altitudinal gradients in RCM forcing and heat fluxes. An analysis of the downscaled dataset confirmed altitude-dependent trends associated with several climate variables (Fig. 1). For RCP8.5, averaged air temperature trends across all RCM simulations are~20% larger at 1700 m a.s.l (0.54°C decade −1 , Fig. 1c) than those at 450 m a.s.l. (0.44°C decade −1 , Fig. 1a). Conversely, surface downward solar radiation (305-2800 nm) is projected to decrease at high altitudes (−0.84 W m −2 decade −1 at 1700 m a.s.l. for RCP8.5, Fig. 1f) but not at low altitudes (Fig. 1d). Changes in shortwave solar radiation can in principle result from changes in cloud cover, in atmospheric water content, or in aerosol concentrations. A recent analysis showed that global climate models (GCMs) generally project significantly increasing surface downward shortwave radiation over Europe due to projected reduced anthropogenic aerosols 40 . However, most EURO-CORDEX RCMs do not take into account the future evolution of these aerosols. The projected spatially variable shortwave radiation in these models, with decreasing radiation at high altitudes and latitudes (Fig. 1i, j) and increasing radiation elsewhere, is therefore likely related to the spatial variability in projected trends of cloud cover, cloud types and atmospheric water content in the RCM simulations. In the alpine region of Switzerland, this results in a significant negative altitudinal gradient of the solar radiation trend of −0.63 W m −2 decade −1 km −1 . For parts of the pre-alps, Jura Mountains and most mid to low altitudes (<~1000 m a.s.l.) in south-west Europe, the RCMs show a positive trend for solar radiation (Fig. 1h-j). The projected wind speed reductions of 0.003 (low altitude) to 0.008 (high altitude) m s −1 decade −1 for RCP8.5 are in line with historical observations for the alpine region but notably smaller than observed recent atmospheric  Table 1 for details) modelled using the Simstrat physical lake model (also in h). Decadal model median trend in surface downward solar radiation (RCP8.5) from 17 regional climate models for h Switzerland and i Europe, together with j regional European altitude vs. latitude dependent trends (10W to 40E, 35N to 70N). Maps created with QGIS v3.4, geographical data obtained for topography from www.gadm.org version 2.8, country boundaries www.naturalearthdata.com version 3.1.0, and lakes www.diva-gis.org/datadown. stilling in central and northern Europe 24 . Trends are also projected for some other atmospheric variables used in the forcing of the Simstrat model ( Supplementary Fig. 2), but their effects on lake temperatures are comparably small as discussed in the following paragraph.
Lake surface temperatures vary in equilibrium with atmospheric forcing, typically increasing by~0.8°C for a 1.0°C change in air temperature, by~0.4°C for a 10 W m −2 change in surface downward solar radiation, by 0.1°C for a 0.1 m s −1 decrease in wind speed and by 0.1°C for a 1% increase in relative humidity 41 . These responses are, however, modulated by lake characteristics and mixing conditions. The projected forcing trends ( Fig. 1a-f) therefore imply that for RCP8.5, air temperature trends would cause high altitude lakes to warm~0.1°C decade −1 faster than low altitude lakes. The projected decrease in surface downward solar radiation, however, would reduce this warming rate bỹ 0.035°C decade −1 . Projected trends in wind speed (Supplementary Fig. 2a-c) suggest a surface warming at high altitude of 0.008°C decade −1 and relative humidity trends (Supplementary Fig. 2g-i) a cooling by~0.02 (mid to low altitude) to~0.008 (high altitude)°C decade −1 , respectively. Precipitation is only used in Simstrat for the calculation of the snow cover on ice, which then affects ice growth and melting, and the projected precipitation trends will not have a relevant impact on this process. In summary, the air temperature and solar radiation trends have the largest impact on lake temperatures, and the trends of the other atmospheric variables are therefore not discussed further in this manuscript.
Changes in lake temperature are driven by heat fluxes between air and water, and these fluxes depend on meteorological conditions and lake surface temperatures 41,42 . Average annual heat fluxes among Swiss lakes display an altitude-dependent change in the response to climate change, strongest for RCP8.5 and weakest for RCP2.6 ( Fig. 2). The effect is strongest in mid-to high-altitude lakes (from~800 m a.s.l.), where surface downward solar radiation (H S ) and uptake (H A ) as well as emission (H W ) of infrared longwave radiation change substantially. The heat-flux trends primarily arise due to the reduced duration of lake ice cover (Fig. 3e). As ice cover recedes, lakes absorb more heat in the form of incoming longwave radiation (H A ) and surface downward solar radiation (H S ). These changes also mean increased heat loss by latent heat flux (H E ) and outgoing longwave radiation (H W ). Altitudinal gradients in air temperature and surface downward solar radiation forcing for different RCMs, combined with the large heat flux modification due to loss of lake ice, create a complex altitude-dependence in air-water heat flux variation. Finally, despite the negative trend in solar radiation at high altitudes in the atmosphere (Fig. 1), the net annual absorbed shortwave solar radiation is projected to increase for high altitude lakes due to the loss of ice (Fig. 2). Using the relationship above, the increased absorption of surface downward solar radiation due to loss of ice cover under RCP8.5 at 1700 m a.s.l. (~1.98 W m −2 decade −1 ) leads to an increase in lake surface temperature bỹ 0.079°C decade −1 .
Climate change impact on the thermal structure of lakes. The variation in the projected impacts of climate change on lake properties including thermal structure, ice and stratification results from: (i) the projected greenhouse gas concentration trajectories (RCP climate assumptions); (ii) the variation of projected atmospheric changes by the climate model chains for a given RCP; and (iii) on individual lake characteristics and how they modulate lake response to forcing scenarios. For lake surface temperatures in the late twenty-first century, the RCPs are clearly the most important driver of change (Figs. 3a and 4a and Supplementary Fig. 3a-c). The average projected lake surface temperature warming compared to the reference period (average of 29 lakes in Figs. 3a and 4a) is 3.3°C (standard deviation, SD: 0.28°C) for RCP8.5, 1.7°C (SD: 0.16°C) for RCP4.5 and 0.9°C (SD: 0.11°C) for RCP2.6. The standard deviations given here are a measure for the individual lake response. In comparison, the variation induced by the climate models is measured with the standard deviations of all model runs for each individual lake. These are on average 0.57°C for RCP8.5, 0.35°C for RCP4.5 and 0.27°C for RCP2.6, indicating that the variance induced by the climate models is about twice as large as the variation between individual lakes. A slight altitudinal difference exists ( Fig. 3a and Supplementary Fig. 3c) with RCP8.5 causing an average lake surface temperature increase of 3.26°C (SD: 0.31°C) for high altitude lakes (>1600 m a.s.l., #27-29), 3.58°C (SD; 0.33°C) for mid altitude lakes (800-1600 m a.s.l., #21-26), and 3.29°C (SD: 0.23°C) for low altitude lakes (400-800 m a.s.l., #1-20).
Lake bottom temperatures are projected to warm at slower rates, estimated (average of 29 lakes in Figs. 3b and 4b) as 1.6°C (SD: 0.87°C) for RCP8.5, 0.93°C (SD: 0.59°C) for RCP4.5 and 0.48°C (SD: 0.37°C) for RCP2.6. For bottom temperatures, the variation between lakes is larger than that induced by the climate models for individual lakes, which amounts to 0.39°C for RCP8.5, 0.32°C for RCP4.5 and 0.25°C for RCP2.6. Remarkably, all emission scenarios result in lower bottom-water warming rates in smaller lakes (<0.15 km 3 and <8.5 km 2 #10, 14, 18, 19, 21-29; with RCP8.5 modelled median changes from 0.06 to 2.10°C, Fig. 4b). The reason for this is that despite climate change, the cooling phase in fall and winter will remain strong enough to extract the amount of heat accumulated in summer, and cool the entire lake down toward the temperature of maximum density of 4°C. The lake bottom temperature is thus reset every winter independent of the exact meteorological conditions. This effect is strongest at high altitude lakes, resulting in diminished bottom temperature warming rate at high altitudes. For RCP8.5 the average projected deep water temperature increase is 0.18°C (SD: 0.10°C) for high altitude lakes, 1.02°C (SD; 0.71°C) for mid altitude lakes and 2.03°C (SD: 0.6°C) for low altitude lakes ( Fig. 3b and Supplementary Fig. 3f).
Future stratification duration in lakes. Our model results suggest that differential warming of surface and deep waters can influence the strength and duration of a given lake's stratified period. Results show substantial changes in the durations of summer and winter stratification (Fig. 3c, d). . The corresponding values for RCP2.6 are 13 days (SD: 2 days) for high-altitude lakes, 16 days (SD; 9 days) for mid-altitude lakes and 7 days (SD: 9 days) for low-altitude lakes. Estimates of the duration of winter stratification show similar altitude-dependence, where highaltitude lakes exhibit the largest decrease (Fig. 3d). RCP8.5 results in an average decrease in inverse stratification of 73 days (SD: 14 days) at high altitudes, 22 days (SD; 19 days) at mid altitudes and 14 days (SD: 15 days) at low altitudes. The corresponding values for RCP2.6 are 21 days (SD: 2 days), 12 days (SD; 9 days) and 6 days (SD: 1 day). Projected changes in summer stratification also depend on lake size (Fig. 4c), with smaller lakes (<0.5 km 3 ) changing faster than larger lakes (>0.5 km 3 ), which is in line with the slower bottom temperature changes in these lakes. There is no clear size dependence of the projected trends in inverse stratification (Fig. 4d).
Climate impact on ice cover duration. During the twentieth century, ice cover for larger and deeper lakes on the Swiss plateau (between 400 and 700 m a.s.l.) drastically diminished to a point of near disappearance 43 . Model estimates for the twenty-first century indicate that persistent ice cover only develops on small, mid-to high-altitude lakes. These lakes experience shorter icecover duration (from #14 to 29; Fig. 3e). The ice cover duration (average of 12 lakes, Figs. 3e and 4e) is projected to decrease by 54 days (SD: 35 days) for RCP8.5, 34 days (SD: 18 days) for RCP4.5 and 18 days (SD: 7 days) for RCP2.6 in the late twentyfirst century compared to the reference period.
As for inverse stratification, the RCP8.5 scenario has the most severe impacts on ice cover for lakes at high altitudes. The average decrease of ice cover duration under RCP8.5 is 85 days (SD: 24 days) at high altitudes, 44 days (SD; 39 days) at mid altitudes Even though the rate of change is faster at higher altitudes, lower altitude lakes risk complete loss of ice cover. Of the three lakes that are projected to completely lose ice cover during late twenty-first century (#14 with RCP8.5, #22 with all RCPs, and #24 with RCP4.5 and RCP8.5), one resides at midaltitude and two at low altitude. Lake volume cannot be linked to ice loss (Fig. 4e), since only the smaller lakes in this study are icecovered in winter. Loss of ice leads to seasonal increased warming mainly in spring, with consequences for stratification duration and mixing regimes as discussed below.
Changes in lake mixing regimes. Changes in the thermal structure of lakes also modify mixing regimes. For example, a dimictic lake that presently experiences seasonal ice-cover will first lose ice cover under climate warming. This subsequently prevents winter stratification and shifts the lake over to a monomictic state. Depending on its inherent tendency to undergo complete mixing, which is a function of its morphology and exposure to wind, a lake may either remain a warm monomictic lake or further shift to an oligomictic or meromictic state with additional warming. Table 1 and Fig. 5 show this shift in mixing regime along an altitudinal gradient for four typical lakes. High altitude lakes such as Lake St. Moritz (#27, 1768 m a.s.l., Fig. 5b) are projected to remain dimictic under all climate scenarios although with reduced duration of winter stratification and ice cover. This loss of ice will cause a greater increase in early summer lake surface temperatures and duration of summer stratification relative to that estimated from air temperature trends alone. Mid altitude lakes, such as Lac de Joux (#23, 1004 m a.s.l., Fig. 5c) or Klöntalersee (#21, 848 m a.s.l., Fig. 5d), are projected to partially or completely shift from dimictic to monomictic under some climate scenarios. This shift clearly depends on greenhouse gas concentration trajectories. The model indicates seven out of the eight dimictic lakes will shift from predominantly dimictic to primarily monomictic regimes during late twenty-first century under RCP8.5. Only three lakes make this projected transition under RCP2.6 and five make it under RCP4.5. The larger, low altitude lakes that already exhibit monomictic regimes will remain in their present states ( Table 1). The present study focuses on shifts from dimictic to monomictic mixing regimes. Due to lack of modelled instances, we did not include a detailed analysis of shifts among deep lakes  Table 1. ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-021-00106-w from monomictic to oligomictic states or the frequency of deep mixing events among oligomictic lakes. Our results generally show that mixing regimes will shift at mid altitudes and that this shift will be larger under scenarios of higher greenhouse gas concentration trajectories (Fig. 5a). Our study suggests that the global distribution of lake thermal regions will not only shift across latitudes 10 but also across altitudinal gradients. This shift is likely to obscure historical definitions of lake mixing regimes based on altitude and latitude 2 especially under RCP8.5. New mixing regime criteria should also consider the effects of water clarity and bathymetry recorded by this and other studies 44,45 .
By regulating heat storage, mixing regimes are important modulators of climate impacts on lakes 18 . A shift from a dimictic to a warm monomictic regime increases the heat storage potential of lakes which do not cool below 4°C. A reduced frequency and/ or intensity of deep water renewal can cause hypoxia or anoxia 46 and thereby increase deep water phosphorus content 47 . Reduced oxygen availability can severely reduce the habitat for higher life forms. Along with temperature increases, these may especially affect cold water fish, who face elevated temperatures from above and diminished oxygen concentrations below 6,48 . Shifts in mixing regime with climate change are one of the most fundamental changes that can take place in lakes. The present study outlines an important altitude dependence of these shifts, which appear most pronounced for mid altitude lakes.

Methods
Climate forcing. In this study, we use RCM simulations from EURO-CORDEX 49 . The RCMs cover the European region and are forced along their boundaries using information from GCMs ( Supplementary Fig. 1) for the period 1981-2099 for three different future climate scenarios. Within the RCM domain, the climate is allowed to evolve freely and includes variables for formation and breakdown of terrestrial snow and ice.
The climate data 50 used to model lake responses in this study were developed under the Swiss climate scenarios framework, CH2018 (ref. 32 ). In this framework, RCM simulations were statistically downscaled from the model grid to local (station) scales and aligned to atmospheric observations through non-parametric empirical quantile mapping to remove RCM biases 32,51-54 . To do this, a day-ofyear dependent correction function was created by relating the statistical distribution of the observed station data to the overlying model grid cell during the calibration period 1981-2010. The correction function was then applied to the entire simulated period 51 , implicitly removing the RCM biases 52,53 and linking the grid scale with the local scale 54 .
Here, we use daily data for mean air temp (tas), precipitation (pr), global radiation (rsds), relative humidity (hurs) and wind speed (sfcwind) covering the Fig. 4 Median changes in projected lake thermal properties of lakes from the reference period  to the late twenty-first century (2071-2099) under three RCP scenarios. Changes in a lake surface temperature (at 1 m depth) and b bottom temperatures (1 m above lake bottom), c duration of summer stratification, d duration of winter stratification and e duration of ice cover for scenario RCP8.5 (orange), RCP4.5 (green) and RCP2.6 (blue). Lakes are ordered by volume as denoted by the circular arrows from the smallest (#1; 0.003 km 3 ) to the largest (#29; 89 km 3 ). For lake details see Table 1.  Table 2 lists station details, and the pairing of meteorological stations with individual lakes is given in Table 1. Using downscaled RCM data at meteorological stations ensures consistency and representativeness for the lake model, which was calibrated with observations from the stations. For four lakes (# 23, 24, 28, 29), some of the downscaled atmospheric variables were not available at the closest station. In these cases, data were combined from two nearby stations (multiple stations, Table 1). To maintain consistency, measured atmospheric variables used in the calibration were combined in the same manner. Air temperature was adjusted for station altitude differences as −0.56°C for every 100 m increase in altitude (known gradient for Switzerland 28 ).
The ensemble median Φ is provided on the EUR-11 grid and combines N model simulations originally on the EUR-11 (Ψ) grid and M model simulations originally on the EUR-44 (ϕ) grid. N and M vary depending on climate scenario (see Supplementary Fig. 1). Overbars indicate regridded data, while k and i indicate whether the data was originally derived from the EUR-11 or EUR-44 grid, respectively. The term ik indicates the use of information from both resolutions. First, the EUR-11 simulations were regridded to the EUR-44 grid ( Ψ n i ) using a bilinear distance weighting. The median (Q) was calculated for all simulations on the EUR-44 grid, and then regridded to EUR-11 resolution as indicated by the first term on the right-hand side of Eq. (1). Second, in order to retain the high-resolution features, Ψ n i was regridded back to the original EUR-11 grid. This resulted in low resolution versions ( Ψ ik ) of the high-resolution simulations (Ψ k ) for the EUR-11 grid. The median was then calculated from the anomalies in the EUR-11 simulations relative to the low-resolution version (Eq. (2)). These were added to the combined median estimate from EUR-44 and EUR-11. The physical lake model Simstrat requires air temperature, wind speed, precipitation, vapour pressure, downward shortwave solar and either downward longwave radiation or fractional cloud cover as forcing variables. These variables were obtained from the CH2018 framework downscaled scenarios, yet two variables were not available, downward longwave radiation and fractional cloud cover. Therefore, shortwave solar radiation was used to estimate fractional cloud cover, c, which was then used by Simstrat to calculate downward longwave radiation. The fractional cloud cover was estimated as where P sf is the potential solar fraction defined as P sf = I s /C ssr . In this expression, I s is observed (or downscaled RCM) shortwave solar irradiance and C ssr is clear sky solar radiation. Equation (3) linearly interpolates cloud cover between 1.0 (P sf equal to the clearness index for complete cloud cover, k cld ) and 0.0 (P sf equal to the clear sky index, k clr ). These two indices depend on location 55 and were obtained during the calibration period (lake dependent; Supplementary  Table 3) using PEST (model-independent parameter estimation software; http://pesthomepage.org). Parameters selected for calibration and the rationale for their selection have been described previously 31 . The remaining model parameters are set to standard values. To match the temporal resolution of climate projections, calibrations used daily averaged measurements from nearby meteorological stations (Table 1). To check the sensitivity of the daily temporal resolution of the forcing data, Simstrat was additionally calibrated using hourly instead of daily averaged forcing. This gave similar model performance for both time increments (Supplementary Fig. 5).
Among available one-dimensional lake models, Simstrat has been shown to generate accurate results for both surface and bottom waters in deep lakes 60 . The original Simstrat model has a tendency to overestimate the internal wave energy in winter and thus the occurrence of mixing events in deep lakes. This can be avoided by filtering out high-frequency wind events (removing wind events shorter thañ ¼ the first internal wave mode period) 61 or by using seasonally varying parameterization for α 7 . In order to limit the number of parameters and manage computational resources, the model used seasonally varying α only for lakes deeper than 150 m, an approach which adequately improved model results (Supplementary Table 3). The point where α_w (winter) changes to α_s (summer) can be assigned to a specific day of the year 7 , but this transition is expected to change with changing climate. The model therefore switches α_s to α_w once thermal stratification disappears, i.e., if the maximum Brunt-Väisälä frequency (N 2 ) during a time step drops below a set limit. N 2 is equal to −g/ρ*dρ/dz where g is gravitational acceleration, ρ is water density, and z is the vertical axis pointing upward. The N 2 limit depends on the grid resolution of the model (0.5 m). N 2 was set to 2 × 10 −4 s −2 based on temporal evaluation of maximum values for N 2 during the reference period (1982-2010). Using real measured atmospheric forcing of Lake Geneva between 1981 and 2011 (31 years), the seasonal α method resulted in 16 modelled overturn events compared with 15 observed events (detected in monthly CTD profiles).
Calculations. According to the conventions of the CH2018 climate analyses, a 30year reference time period from 1981 to 2010 was used and compared to the late twenty-first century period from 2070 to 2099. However, since the first simulation year could include spin-up effects, we dropped it from the analysis and used 1982 to 2010 as the reference period instead. To be consistent, the first year was also dropped for the late twenty-first century, and averages were calculated from 2071 to 2099. The approximate thirty years' classical time frame aims at removing random effects due to interannual variability 33 . All calculations were performed in MATLAB R 2017b and statistical properties such as trends, means, medians and variances were calculated first for each individual RCM and then for each RCP. For example, in Fig. 3a, the median surface temperature change from the reference period to late twenty-first century was calculated as follows. First, annual means (temperatures) and annual sums (stratification, ice) for each RCM simulation were obtained. Second, differences between the two periods were calculated for each RCM simulation. Third, medians, means and variances across the models (7 models in RCP2.6 and RCP4.5, and 17 for RCP8.5) of this difference were calculated (Figs. 3 and 4 and Supplementary  Fig. 6).
A lake is considered to be stratified if the temperature difference (ΔT) between the surface water (10 m mean) and the bottom water (lowest point) exceeds 1°C 62 . The surface mean temperature was averaged over the top 10 m to remove shortterm temperature fluctuations. Lakes investigated in this study are all deep (Table 1), and water exchange between surface and bottom is limited when a significant density gradient develops. In summer ΔT > 1°C results in stable stratification, while in winter ΔT < −1 describes an inverse stable stratification. Multiple stratified periods in both summer and winter, separated with short gaps (<15 days) due to strong short-term mixing events, were merged prior to estimating the total annual stratification duration.
Future and past mixing regimes were classified as either warm monomictic (overturn once in winter), dimictic (overturn twice per year, in autumn and in spring) or oligomictic (irregular overturns, not occurring every year). A lake overturn is defined to have occurred if surface temperature (10 m average) becomes colder or equal to bottom temperature after a previous summer stratification (monomictic, oligomictic and dimictic lakes), or if the surface temperature becomes warmer or equal to bottom temperature following an inverse stratification (dimictic lakes).
Stratification regimes were obtained annually (Fig. 5a), and as periodic means during the reference period (1982-2010) and for the late twenty-first century (2071-2099) ( Table 1 and Fig. 5b-e). A lake was first classified for each simulated year as dimictic, monomictic or meromictic (no mixing). The final classification was then defined as dimictic or monomictic, if that was the regime occurring in most years for each climate scenario, or as oligomictic, if it was classified as meromictic in most years.