Lake Geneva is a large deep perialpine lake located at the border between France and Switzerland at 372 m above sea level (46.45° N, 6.53° E). It is the largest lake in western Europe with a mean depth of 153 m and a maximum depth of 309 m. Its surface area is 582 km2, its volume is 89 km3, and its hydraulic retention time is around 11.4 years. It is thermally stratified from spring to early fall with a thermocline located at about 15 m depth in summer. Complete water mixing promoting deep reoxygenation occurs during very cold winters only once per five to ten years on average according to monitoring since the 1950s20. We selected Lake Geneva to apply our methodology because of its representativeness as it is one of the most studied and well-known lakes in the world, with a rich high-quality monitoring database spanning over 60 years. In addition, since the 1950s, it has experienced a history of anthropogenic eutrophication followed by re-oligotrophication similar to that of several other lakes over about the same time period11. And like other lakes, Lake Geneva faces the ongoing effects of climate change with significant warming trends and strengthening thermal stratification25.
In situ monitoring data
Monitoring station SHL2 (46.45° N, 6.59° E) is located at the deepest part of Lake Geneva, called the Grand Lac, representing more than 96% of the lake’s total water volume, where various physico-chemical and biological measurements have been conducted since 1957 by the Centre Alpin de Recherche sur les Réseaux Trophiques des Ecosystémes Limniques (CARRTEL). Physicochemical and biological data from probe sensors and analysis of water samples are freely available in the Observatory on Lakes database16. Discrete vertical measurements are taken from the surface to the bottom (at 0, 2.5, 5, 7.5, 10, 15, 20, 25, 30, 35, 40, 50, 100, 150, 200, 225, 250, 275, 280, 285, 290, 295, 300, 305, and 309 m depth) every two weeks, except from November to February, when sampling is performed monthly. The major tributary to Lake Geneva is the Rhône River with an average discharge of 186.4 m3 s–1 accounting for 75% of the lake’s total inflow volume. The Rhône River is monitored by the national long-term surveillance of Swiss rivers13 typically weekly since 1974 at the Porte de Scex station (46.35° N, 6.89° E), situated 5 km upstream of where the river enters the lake. Field data encompass the inflow discharge and several physicochemical variables. The outflow discharge of the Rhône River has been measured monthly since 1919.
Several field campaigns to collect sediment cores at the deepest point of the lake were conducted between 2004 and 2012. Two paleolimnological datasets obtained from well-dated sediment cores were used in the present study: (a) mean inferred-TP (total phosphorus) concentrations reconstructed from the absolute abundance of Daphnia remains using a transfer-function for 1850–2018 developed for large and deep perialpine lakes14. Daphnia has been successfully used as a quantitative indicator of lake trophic level due to its strong stoichiometric requirements for phosphorus30. Generalized additive models were used to smooth temporal data; and (b) the annual volumes of hypoxic waters reconstructed using a sedimentological proxy based on the formation of varves, i.e., annually laminated sediments17. This proxy relies on the assumption that when oxygen concentrations in the water-sediment interface fall below a critical point, determined by a combination of sufficient time and intensity, microbenthic life disappears, which prevents bioturbation and its related sediment mixing, thereby leading to the preservation of varves. However, because of the plurality of ecosystem components and their individual reactions, the term ‘hypoxia’ must refer to a process or low-oxygen stress rather than to some general application of a specific concentration threshold. Hence, the term ‘hypoxia’ is used in terms of the sedimentary consequence of a decrease in oxygen availability in bottom waters, combining duration and oxygen concentration.
Climatic data including air temperature at 2 m above the lake surface, wind speed at 10 m, and surface solar radiation were downloaded from a statistically bias-adjusted and downscaled climate model from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3b12). The data were previously verified against measurements from a local meteorological station located at Thonon-les-Bains (1987–2019 and 1971–2019 for air temperature and shortwave radiation, respectively, monitoring data from the INRAE CLIMATIK platform) and a scaling factor was implemented to correct the altitude bias25. Projections used historical climate from 1850 to 2014 and three future scenarios from 2015 to 2100: SSP1-RCP2.6, SSP3-RCP7.0, and SSP5-RCP8.5. The SSPs (Shared Socioeconomic Pathways) consider how societal choices will affect greenhouse emissions (SSP1 being the most sustainable scenario and SSP5 the least); and the RCPs (Representative Concentration Pathways) correspond to the range of the year 2100’s radiative forcing values, from 2.6 to 8.5 W m–2. These data were extracted at a daily resolution for the grid cell (55 km 55 km) containing the lake. The other meteorological variables for which local downscaling often lacks accuracy (relative humidity, cloud cover, and rainfall) were extracted from meteorological observations31 located at Nyon/Changins from 2000 to 2011 (MétéoSuisse Data Warehouse), from which daily means were calculated and replicated every year from 1850 to 2100 as we currently lack clear expectation for their fate under future scenarios. This approach of reducing model input variables to only those with a high confidence level was previously applied to Lake Geneva20 and was found to be well-adapted to long-term simulations by overcoming limitations on their scaling and by their negligible effects on the hydrodynamic processes25.
Historical reconstruction of long-term driver data
GLM-AED2 hydrodynamics model platform requires daily meteorological, hydrological, and nutrient load data as inputs. Due to monitoring data temporal insufficiency, a conceptually sound framework was developed to estimate long-term input forcing data able to capture seasonal variability over the simulation period. This framework accounts for the effects of climate change and eutrophication, the two most prevalent stressors of our study site, relying on (i) long-term climatic projections, (ii) long-term monitoring in Lake Geneva of total phosphorus (TP) and soluble reactive phosphorus (SRP) over 64 years (1957–2020), (iii) long-term monitoring in the Rhone River of discharge, temperature, electrical conductivity, dissolved oxygen, TP, SRP, ammonium, nitrate, total nitrogen, reactive silica, total and dissolved organic carbon over 47 years (1974–2020), and (iv) TP inferred from paleolimnological records from 1850 to 2018 extracted from sediment cores taken at the deepest point in the lake.
TP conditions in the lake were used to reconstruct the long-term trend of TP inputs to Lake Geneva in the absence of longer temporal trends in the lake tributaries. Firstly, TP records inferred from sediment cores at an annual scale from 1850 to 2018 were multiplied by a scaling factor of 2.64 to convert TP in the lake into TP in the tributaries; then the values were proportionally divided into SRP, dissolved and particulate organic phosphorus in the river, as required by the model. The scaling factors were calculated based on the ratios between median values of field measurements over the monitoring period (Fig. S1a). The contribution of phosphorus inputs (SRP, dissolved and particulate organic phosphorus) from the other tributaries was summed and aggregated into the Rhône River as a single inflow into the lake32 by applying a scaling factor of 1.33 (based on their contribution of 25% in the total discharge into the lake). Secondly, weekly measurements in the Rhone River from 1974 to 2020 were linearly interpolated to obtain a daily time series. This method has been used in other model applications33 but potentially underestimates the effect of storm events that may not be captured by routine monitoring. The mean daily concentration was normalized by the mean annual concentration over the period to create a synthetic seasonality
(Fig. S1b). Then, the daily variation was applied to the mean annual concentrations derived from paleolimnological records to reconstruct a synthetic long-term daily time-series of phosphorus input into Lake Geneva from 1850 to 2018 (Fig. S1c). Considering the long water residence time of 11.4 years in Lake Geneva, the phosphorus inputs were adjusted backward in time by 11 years. Finally, a correction factor of 0.64 was applied for achieving a better match of the TP peak magnitude. The adoption of correction factors to multiply input variables are commonly used (e.g., for wind factor) in the face of uncertainty in their reconstruction.
The absence of paleolimnological information for the remaining mass inputs into the lake (ammonium, nitrate, and reactive silica) prevented the reconstruction of their long-term trends. In this manner, the annual means from field measurements in the Rhone River were assumed as constant and a similar procedure to account for intra-seasonality based on weekly measurements was applied. Although this approach is a pragmatic compromise for detailed trajectories of nutrient inputs over the past, it provides a reasonable overview of major historical changes in the lake in accordance with negligible long-term change measured from 1974 to 2020 (Fig. S2). As we did for TP, a scaling factor of 1.33 was applied to account for the contribution of their loadings in the remaining tributaries. Long-term inflow and outflow discharges, inflow temperature, salinity, dissolved oxygen, and organic matter concentrations were estimated following common procedures found in the literature based on field measurements. A complete description of each of these methods is provided in the supplementary material.
Hydrodynamic and biogeochemical model description
Oxygen concentration and hypoxia dynamics were simulated in this study using the state-of-the-art Aquatic Eco-dynamics model (AED2), which has been recently applied to lakes worldwide34–38. It was selected because it includes the key processes of water quality dynamics, that depend on both the climatic environment and nutrient loads, being well-suited to studying impacts of eutrophication and climate change on lake internal processes. The AED2 model must be coupled to a host physical (hydrodynamic) driver model to connect the role of thermal stratification and vertical mixing on lake ecosystem dynamics. We used the one-dimensional (1D) General Lake Model (GLM) which has been shown to provide accurate reproduction of the mixing dynamics of lakes across a broad range of latitudes, climatic zones, and morphometric properties39. The one-dimensionality was mandatory for this long-term study given its low computational cost. Although other processes such as differential cooling or riverine underflow push towards the use of 3D models (or in some cases 2D models) in large deep lakes because hydro-meteorological conditions are heterogeneous over extensive water surface areas, 1D models can be sufficient to represent long-term dynamics and provide a good compromise between accuracy and computational cost, which facilitates model implementation, calibration and simulation, especially at high spatial resolution19. Notably, 1D models have proved to be a viable compromise of the representation of physics being able to reproduce the thermal regime with reasonable accuracy in perialpine lakes11. Particularly for Lake Geneva, 1D models have been found to be applicable beyond investigations solely of physical processes, but they are sufficient to attempt effects-oriented modelling for issues such as water quality in a changing climate40.
The hydrodynamic GLM model was previously applied to perform long-term hydrodynamic simulations in Lake Geneva from 1850 to 2100 and accurately reproduced thermal processes after calibration and validation with good agreement with observations (RMSE of
1.02 °C and 0.57 °C in the full profile and hypolimnion, respectively)25. In brief, GLM adopts a flexible Lagrangian scheme of homogeneous horizontal layers that dynamically changes the thickness of each layer as a function of the vertical density gradient in the water column. Within the lake domain, the model computes an energy balance that compares the available kinetic energy to the internal potential energy of the water column41. The vertical profile of temperature is provided by the GLM as input to the AED2 at every time step of the simulation for the computation of temperature-dependent chemical and biological processes; then the AED2 can feedback conditions to the hydrodynamic model by modifying the light extinction coefficient.
The AED2 model is a library of equations for the simulation of aquatic ecodynamics, i.e., water quality, aquatic biogeochemistry, biotic habitat, and aquatic ecosystem dynamics42. It consists of numerous modules that can be connected through specific variable dependencies to represent the most relevant aquatic biogeochemical processes, including the cycling of carbon, nitrogen and phosphorus, oxygen, organic matter, as well as different functional groups of phytoplankton and zooplankton. State variables defined within each module are subject to transport and mass conservation, thereby allowing the quantification of interactions that otherwise would be difficult to infer from observational data alone. Given its flexibility, the model can be applied across a range of scales (from 0D box models to 3D models) and contexts (wetlands, lakes, reservoirs, rivers, and estuaries). A detailed description of the model process equations, state variables, and parameters can be found in the application manuals for GLM41 and AED242.
The GLM-AED2 model (GLM: v.3.1.1, AED2: v.2.0) was run on an hourly time step from January 1st 1850 to December 31st 2020. A continuous run throughout the total simulation was performed to ensure the annual carryover of heat stored in the hypolimnion between full turnovers. The model was configured to simulate the dynamics of dissolved oxygen, inorganic nutrients (carbon, nitrogen, and phosphorus), silica, organic matter (particulate and dissolved), and phytoplankton. The dominant biogeochemical processes in the water column were represented: fluxes between the water-sediment interface, phosphorus internal loading, mineralization of dissolved organic matter, hydrolysis of particulate organic matter, sedimentation of particulate organic matter, nitrification, and denitrification. In addition, the model computed phytoplankton primary production, respiration, vertical movement, excretion, and mortality according to nutrient, temperature, and light limitation. The phytoplankton composition was limited a single group representing the total biomass. The configuration of the GLM model that had been previously calibrated and validated25 was adopted and all AED2 parameters were set at their default values42. As the initial conditions of the state variables in 1850 remained unknown, the initial values of nutrients, organic matter, and chlorophyll-a were set as zero, while the initial vertical profile of DO was set as its mean value from field measurements (1957–2020). The statistical software R 4.1.2 was used to simulate and analyze modelling results.
Model reliability and limitations
The modelling procedure followed well-defined standard techniques, i.e., calibration based on field measurements over a representative time period and validation against an independent dataset43 in accordance with other aquatic ecosystem modeling studies33,36,44, based on a comparison against field measurements collected in the deepest point of the lake over 64 years (1957–202016; Table S1, Table S2, Fig. S4). We used time periods prior to and after the calibration period for validation to stress test the model by applying it at periods with different ecological characteristics. A detailed description of the calibration and validation procedures, as well as the assessment of model performance are provided in the Supplementary material.
As the model performance assessment and acceptability are highly dependent on the model's purpose, we applied traditional goodness-of-fit metrics for evaluating the ability of the model to capture short-term processes in addition to a complimentary analysis to quantitatively assess the model performance in replicating biogeochemical dynamics at long-term scales, the focus of the present study. Some metrics errors were less than established appropriate levels of confidence (RMSE lower than 3 mg L–1 for DO and RRMSE lower than 50% for TOC and lower than 100% for Chl-a45), while other error metrics (e.g., coefficient of determination and Nash-Sutcliffe efficiency) were relatively high compared to other studies that performed short-term simulations. The reason is that many of the former model studies have covered shorter time scales10 and systems that have not undergone such high nutrient variability (TP from 10 to 90 µg L–1), which somewhat eases the ability of the model to capture ecosystem dynamics.
The performance assessment of the model’s ability to represent bottom oxygen and the hypoxia regime over different temporal scales confirms the reliability of the methodological framework in providing robust modelled long-term trends of bottom oxygen in terms of general behavior, amplitude, and timing (Fig. S5, Fig. S6). Discrepancies between model results and observed data were attributed mostly to the inability of the hydrodynamic module to reproduce the timing of deep mixing by using ISIMIP climatic projections, instead of meteorological time series from a nearby station, as DO dynamics at the bottom of the lake are strongly associated with the occurrence of stochastic deep mixing events. Nevertheless, the coupled GLM-AED2 model was able to capture important long-term transitions in bottom DO of Lake Geneva at a decadal scale. For instance, the model corrected predicted reoccurring long-term patterns of dissolved oxygen in the bottom of the lake, such as the decadal DO replenishment after complete mixing events followed by consumption under well-established thermal stratification. Finally, while the model is not suitable for fully representing certain short-term processes accurately, it is nevertheless a robust tool to adequately describe long-term evolution of the hypoxia regime in the lake over 250 years, which was an especially critical gap for lake research. If we want to successfully respond in an adaptive manner, our understanding of future risks needs to embrace these uncertainties46.
Characterising long-term hypoxia
After calibration and validation of the model at a fine temporal scale against field data, modelled TP and oxygen were qualitatively compared to paleolimnological data to ascertain the accuracy of the long-term trend over the period 1860–2018 in terms of amplitude, timing, and behavior. The first decade of simulations (1850–1860) was taken as a “warm-up period”, meaning we did not analyze model outputs over this period. A quantitative performance assessment was not possible because the concentrations in the sediments cannot be directly translated to the concentrations in water. The paleolimnological proxy of hypoxia volume reflects a combination of duration and intensity of low-oxygen in bottom-water, rather than some general concentration threshold. As a result, we characterized hypoxia intensity by calculating daily minimum values from modelled oxygen concentrations throughout the water column summarized as annual averages each year, matching the same frequency as the paleolimnological data. Hypoxia duration was calculated as the number of days per year for which hypolimnetic minimum concentrations were below specific thresholds. We tested the often-used thresholds of 1–4 mg L−1, and found a better match with the paleolimnological proxy at 1 mg L−1, but the other thresholds presented similar patterns.
For analysing the effects of climate on lake hydrodynamics, thereby affecting DO replenishment in the hypolimnion, the occurrence of winter complete mixing was defined as the event of temperature and oxygen concentration difference less than 2 ºC and 2 mg L–1, respectively, between the surface and the bottom of the lake for at least 5 consecutive days. The temperature threshold is in the range of often used values from 1 ºC40 to 3 ºC20 found in previous studies in Lake Geneva.
Future projections of lake hypoxia
We used climatic projections for predicting changes in solar radiation, air temperature, and wind speed in three different future scenarios (SSP1–RCP2.6, SSP3–RCP7.0, and SSP5–RCP8.5). Long-term future projections of phosphorus inputs into the lake from 2020 to 2100 were based on several prospective scenarios to represent different management strategies. Phosphorus external load in 2019 gradually increases or decreases until 2030, after which the average annual value was assumed constant up to 2100. The minimum and maximum phosphorus external loads were between 2 and 100 µg L–1, respectively, to cover the entire range based on the historical condition. The step sizes between the minimum and maximum phosphorus external load varied, with a step size of 1 µg L–1 between 2 and 25 µg L–1 and a step size of 5 µg L–1 between 25 and 100 uµ L–1 (Fig. S8). Although the annual phosphorus load remained constant from 2030, the TP in the lake stabilized as a constant annual value only after 2050 due to the long inertia of such a large
(582 km2) and deep (309 m depth) lake (Fig. S9). As a result, future projections were analysed from 2050 to 2100. All the remaining input data adopted for the historical modelling were assumed constant and the initial conditions of the state variables in 1st January 2020 were set as simulated values in 31st December 2019.
In addition to the model intrinsic limitations and the simplifications adopted in the present study, future projections are based on the assumption that fixed parameter values employed by the model that induce stationarity in process rates are able to sufficiently simulate the fluctuating ecosystem states in the future as they were for the strong changes on TP and oxygen conditions in the past (1860–2020). In addition, simulated future conditions do not consider the effects of quagga mussel populations that are expected to increase in European sub-alpine lakes in the next decades, with potential strong impacts on DO dynamics29.
The areal hypolimnetic mineralization rate (AHM, [g O2 m−2 d−1]) over the whole hypolimnion describes the processes of (1) organic matter mineralization in the water column,
(2) sediment oxygen uptake, and (3) flux of reduced compounds from the sediment. It was computed from modelled mean DO concentrations between 15 m and 309 m depth at a daily basis as the difference between the modelled maximum DO concentration modelled (February to April) and the minimum DO at the end of summer stratification (late October/November) divided by the number of days of summer stratification23.