A crisis of lake hypoxia in the Anthropocene: The long-term effects of climate and nutrients

doi:10.21203/rs.3.rs-3234938/v1

Download PDF

Article

A crisis of lake hypoxia in the Anthropocene: The long-term effects of climate and nutrients

https://doi.org/10.21203/rs.3.rs-3234938/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Climate change is altering thermal stratification in lakes worldwide. Reduction in winter mixing lead to prolonged oxygen depletion, lasting for years to centuries, potentially becoming permanent. Although there is convincing evidence of lake deoxygenation globally, its duration, timing, and impacts over decadal to centennial timescales remain uncertain. Here, we introduce a novel model-data assimilation approach using 150 years of limnological and paleolimnological data to evaluate the anthropogenic impact and future of deep dissolved oxygen in Lake Geneva. We find that climate change has influenced winter mixing, with divergent effects on bottom oxygen concentrations before and after eutrophication. Over centennial timescales, eutrophication, not climate warming, triggered unprecedented bottom-water hypoxia. However, by 2100, climate change will be the main driver of hypoxia in Lake Geneva and similar lakes, even with reduced phosphorus concentrations. With climate change locking in the effects of phosphorus loading on hypoxia, the significance of reducing loading remains intact.

Earth and environmental sciences/Climate sciences/Climate change/Climate-change impacts/Environmental health

Earth and environmental sciences/Climate sciences/Limnology

Earth and environmental sciences/Limnology

More than a century of research, starting at beginning of limnology as a discipline¹, has addressed oxygen dynamics in lakes using empirical, theoretical, and predictive approaches. Declining dissolved oxygen (DO), a major threat to freshwater and marine environments², is attributed to excessive nutrient inputs³ and anthropogenic climate change⁴. In the face of synchronous changes in climate and land use over the last century, understanding temporal dynamics is becoming increasingly important. In the Anthropocene, whereby human activity significantly modifies the planet’s climate⁵, abrupt, large and persistent shifts in the function and structure of ecosystems have made both reliable prediction and managed reversal of change exceptionally difficult⁶. Complex interactions among these regulators, which may act as triggers, contributors, or amplifiers of oxygen depletion, cause nonlinear long-lasting responses in ecosystems that are still largely difficult to predict⁷. Given that oxygen depletion has substantial implications for critical lake ecosystem services mediated by oxygen-sensitive biogeochemical processes, identifying the drivers and mechanisms of changes in oxygen dynamics is critical for understanding lake ecosystem functioning as a whole and mandatory for managing water quality.

Despite multiple concurrent observations on the oxygen decline in lakes worldwide, the interpretations of the drivers and the underlying mechanisms differ between long- and short-timescale approaches. On one hand, clear evidence of persistent oxygen depletion, based on a limnological synthesis of data from 393 temperate-zone lakes, indicates that the decline in the hypolimnion is associated with stronger thermal stratification and decrease of water clarity during the past 40 years⁸. On the other hand, paleolimnological synthesis of data from 365 lakes worldwide suggests that increased human activity and nutrient loading, but not climate, has led to the onset of hypoxia during the last 300 years⁹. To reconcile these apparent contradictions, comprehensive systematic long-term, temporally continuous data, are needed to fill a critical knowledge gap, and will be mandatory for unraveling the precise role of external drivers of ecosystem patterns. Here, we develop an explicit integration of (1) a process-based lake model¹⁰, with (2) in situ monitoring data for model calibration, and (3) paleolimnological proxies stored in lake sediments to provide the missing long-term perspective in continuous temporal information. This combination of techniques, applied over a period of 250 years (1850–2100), represents a novel approach that makes it possible to overcome current limitations in our ability to disentangle natural from anthropogenic, and endogenous from exogenous regulators of bottom lake hypoxia.

The novel methodological framework proposed here (Fig. 1) uses a state-of-the-art one-dimensional vertical process-based model (Aquatic EcoDynamics coupled to General Lake Model, GLM-AED2) to perform continuous simulations of lake biogeochemistry, with a focus on the dynamics of bottom oxygen concentration, driven by both climatic forcing and nutrient loading in Lake Geneva (France, Switzerland). This water body, one of the most thoroughly monitored and well-known lake in the world, with a rich high-quality database of monitoring data for over
60 years, has experienced a distinct history of changing phosphorus input over the last century similar to that of several other lakes worldwide¹¹ (see “Methods”). Model forcing data, applied to iconic Lake Geneva, consist of (a) climatic models of air temperature, shortwave radiation, and wind speed derived from the ISIMIP framework¹² comprising a historical period (1850–2014) and future simulations (2015–2100) of three projections based on socioeconomic pathways and forcing levels of the representative greenhouse gas concentration pathways (SSP1–RCP2.6, SSP3–RCP7.0, and SSP5–RCP8.5¹²); (b) field measurements of discharge, temperature, and nutrients collected in the main tributary of the lake (Rhone River) over 47 years (1974–2020¹³);
(c) paleolimnological proxies, i.e., preserved Daphnia subfossils (Cladoceran) extracted from well-dated sediment cores that allow the development of a Daphnia-based transfer function to predict volume-averaged total phosphorus (TP) concentration in the lake water from 1850 to 2018 based on their strong relationship, also referred to as Daphnia-inferred TP¹⁴; and d) lake-specific hypsography¹⁵. The historical reconstruction of long-term data considered the temporal evolution of driver variables with the strongest influence on biogeochemical in-lake processes while coarser assumptions were made for the other variables necessary for calculating energy and mass balances (see “Methods”). The model’s ability to reproduce the bottom oxygen and hypoxia regime over different temporal scales was assessed by comparison against field measurements from probe sensors and analysis of water samples collected in the deepest point of the lake every two weeks¹⁶. Hypoxia regime was further validated using paleolimnological record of water hypoxia, expressed as annual volume of hypoxic waters, derived form a varve index, i.e. annually laminated sediments¹⁷. Both assessments confirmed the reliability of our methodological framework in providing robust modelled long-term trends of bottom oxygen in terms of general behavior, amplitude, and timing (see “Methods”).

Simulated daily oxygen concentrations throughout the water column from 1850 to 2100 are summarised as annual values to characterise the hypoxia regime. We derived metrics for intensity (annual minimum oxygen concentration throughout the water column) and duration (number of days over the year with a minimum oxygen concentration below 1 mg L^–1; see “Methods”). We also investigated the primary drivers of hypoxia by disentangling the effects of climatic forcing and nutrient-enrichment in order to isolate the influence of each, both of which have been considered important for the oxygen content in the bottom of the lake. We did this by running the model with detrended phosphorus input, taking it constant at its natural background over the simulation period, while all the remaining forcings were unchanged. Finally, we ran
117 future projections, representing the evolution of the hypoxia regime in the lake, subject to three different scenarios (the low, medium, and high greenhouse gas emission), by assuming a wide range of phosphorus inputs into the lake from a low-input value until the maximum historical input value (see “Methods”).

Our daily modelled results capture the remarkable phosphorus history in the lake for both high- and low-nutrient conditions (Fig. 2a). Baseline conditions of total phosphorus, before the 1950s, revealed the range of inherent natural variability (4.1 ± 4.6 µg L^–1). Starting in the mid-20^th century, a dramatic TP increase up to 66.5 µg L⁻¹ in model simulation following the demographic and economic development in the region resulted in eutrophic conditions (TP > 24 µg L^–1)¹⁸ from 1952 and reaching a TP peak in the 1970s. From then, remediation measures including the construction of wastewater treatment plants and a ban on phosphates in detergents were taken to control phosphorus point sources and diffuse inputs. These resulted in a consistent TP decline termed “re-oligotrophication”. The amplitude of the simulated annual averages (0.4–66.5 µg L^–1) was comparable to that of the paleolimnological proxy (0–63.6 µg L^–1), although both are a smooth representation of the long-term sharp changes in TP over time as revealed by field measurements.

The magnitude and duration of hypoxia were characterized by modelled daily minimum oxygen concentrations throughout the water column summarised as annual averages each year and by the number of days per year for which hypolimnetic concentrations were below 1 mg L⁻¹. Model simulations provide important insight into reference oxygen conditions by hindcasting the range of its inherent natural variability before the monitoring period (i.e., 15 ± 51 days yr^–1 of hypoxia and 5.5 ± 2.0 mg L^–1of minimum DO throughout the water column), something that is not captured by the paleolimnological proxy. The reconstructed hypoxia dynamics demonstrated that Lake Geneva was well-oxygenated before the second half of the 20^th century exhibiting only short-duration, episodic events of hypoxia, i.e., not longer than one year. The lake exhibited resilience to the establishment of hypoxia until the 1970s. An early signal of oxygen dynamics change is provided by the paleolimnological proxy in the 1960s. It was only in 1973, when median TP was as high as 67 µg L^–1, that instead of a smooth continuous trend following nutrient enrichment, the hypoxia regime experienced a sharp shift in duration and magnitude unprecedented in the previous 112 years (1860–1972; Fig. 2b, c). Indeed, between 1973 and 2020, the lake was hypoxic year-round in 21 years while DO averaged lower than 1 mg L^–1 in 34 of those years. In contrast, from 1860 to 1972, the lake did not experience a single year with those conditions.

Our results confirm the role of historical anthropogenic phosphorus over-enrichment as the primary environmental forcing that triggered the onset of hypoxia which is, in turn, a product of microbial respiration and consumption of sinking organic matter. By running the model with detrended phosphorus input, i.e., at its background concentration, while keeping the trends in historical climatic forcing, we assess the hypoxia regime as influenced solely by climate, with no nutrient-enrichment influence. Under this condition, very low TP concentrations in the lake
(TP < 0.8 µg L^–1) were obtained in the model (Fig. 3a). Minimum DO in the lake lower than
1 mg L^–1 occurred at a lower frequency (in 22 years; Fig. 3b) and, more importantly, annual minimum DO of 2.30 mg L^–1 on average from 1973 to 2009 (Fig. 3c), the period of persistent hypoxia over the historical condition (with 34 years of hypoxia and minimum DO of 0.60 mg L^–1 on average), revealed that the hypoxia onset must be the result of excessive catchment phosphorus load. However, both the intensity and duration of hypoxia can be explained by the combined influence of climatic forcing and high phosphorus loads since DO experienced lower concentrations starting in the 1970s accompanied by longer periods of hypoxia even in the absence of eutrophication.

The resistance mechanism preventing the spread of hypoxia in the 1950s was mixing events sufficiently effective to compensate for high-nutrient levels which maintained hypolimnion oxygenation from the 1950s to 1970. In seasonally stratifying lakes, the dominant process driving the replenishment of dissolved oxygen in the hypolimnion is winter deep mixing⁷. Indeed, from 1860 to 1970, complete mixing (see “Methods”) occurred during 67 winters (dotted lines in Fig. 3), being the dominant source of DO mass discharges from the surface to the hypolimnion. Since 1970, deep mixing has become intermittent mixing (on average 1 event dec^–1 in 1971–2020) according to our model results and limnological data, that severely limited the replenishment of bottom oxygen. Despite some mismatches between the timing of modeled and observed mixing events, the frequency is comparable (1.9 event dec^–1 in model results against 1.2 event dec^–1 in monitoring data for 1971–2020)^19,20. As mixing events occur at lower frequency, unprecedented hypoxia in the hypolimnion extended over longer intervals, generally exceeding one complete year. The worldwide decreased frequency of mixing is an on-going process in Lake Geneva which started in the 1970s due to warming, predicted in the future according to numerical simulations and different IPCC emission scenarios²¹, and lower wind speed over almost all continental areas in the northern mid-latitude since 1979²².

Having demonstrated the long-term trajectory of historical hypoxia dynamics in the lake’s hypolimnion, over the past 160 years (1860–2020), we investigated projected changes under future climatic forcing scenarios (SSP1–RCP2.6, SSP3–RCP7.0, and SSP5–RCP8.5) and prospective scenarios of phosphorus external load by varying phosphorus inputs to explore the outcome of different management strategies from 2020 to 2100 (see “Methods”). We used minimum and maximum phosphorus concentrations (resulting from loading differences) between 2 and
100 µg L^–1, respectively, to cover the entire range of loads based on the historical condition. According to the modelled dynamics, the severity of change in the hypoxia regime by the end of the twenty-first century will be primarily driven by climatic forcing for the three scenarios, despite extreme variation in phosphorus concentration in the lake (Fig. S10). We project that the lake will experience hypoxia lasting longer than 15 d yr ^–1and minimum DO lower than 5.5 mg L^–1 (average duration and concentration for pre-1950). As expected, poorer conditions were found in the high-emission scenario (SPP 5–RCP 8.5) and were explained primarily by the magnitude of change in the climatic drivers investigated (air temperature, solar radiation, and wind speed).

Our model projects that changes in the annual duration of bottom hypoxia and the accompanying oxygen concentrations will not be as strong as they were in the past since the shift to hypoxia (1973–2020). Indeed, the deepest layers (270–300 m) are expected to experience an annual average daily minimum oxygen concentration of 2.0–3.0 mg L^–1 (volume-weighted values), and hypoxia lasting 134–226 d yr^–1 by the end of the century. This is in contrast to the contemporary average concentration of 1.2 mg L^–1 lasting 246 d yr^–1 (Table 1). Although the impacts of future climate warming will not be as strong as the impacts of eutrophication in the past, the projected changes should not be underestimated as they still represent a loss of 43–62% of average daily minimum DO concentration, an additional 104–196 days of hypoxia in comparison to the natural condition (pre-1950).

Table 1 | Modelled mean hypoxia duration and annual minimum oxygen concentration throughout the water column.

Indicators	Natural	Contemporary	Future projections (2050–2100)
	(pre-1950)	(1973–2020)	SSP1–RCP2.6	SSP3–RCP7.0	SSP5–RCP8.5
Hypoxia duration (d yr^–1)	15 (±51)	246 (±150)	134 (±82)	166 (±84)	226 (±72)
Minimum DO (mg L^–1)	5.5 (±2.0)	1.2 (±2.0)	3.0 (±1.9)	2.6 (±1.6)	2.0 (±1.3)

We computed the areal hypolimnetic mineralization rate (AHM) over the whole hypolimnion, that describes the processes of water column mineralization, sediment oxygen uptake, and flux of reduced compounds from the sediment, a strong indication of the nutrient effects in the bottom of the lake. Our results show AHM oscillation from 0.3 to 1.4 g O₂ m⁻² d⁻¹ over the historical period in the range of previous studies²³ with exceptionally high values
(1.5–2.3 g O₂ m⁻² d⁻¹) during the peak of eutrophication (1967–1996). The logarithmic relation between the number of years without complete mixing and AHM reveals a stronger effect of climate in the future (2020–2100), as a narrower range of AHM would be required to maintain well-oxygenation in the deep hypolimnion under less frequent mixing events (Fig. 4). We therefore conclude that oxygen replenishment during complete mixing events will be the governing driver controlling the severity of hypoxia in deep hypolimnion. This analysis demonstrates that over the 19^th and 20^th centuries, hypoxia trajectories have been predominantly driven by the combined effects of climatic forcing and phosphorus concentrations, whereas in the 21^st, it will be primarily under climatic control affecting physical processes that regulate thermal dynamics, and less by nutrient content that influences organic matter and its decomposition. Given the projected worldwide increase in stratification intensity and duration of deep lakes as a response to climate warming^24–26, prolonged isolation of the hypolimnion from the atmosphere is expected. As a result, lakes will be increasingly susceptible to hypoxia because of less frequent mixing with more severe consequences expected for the systems with a eutrophic past. Our findings strongly support the need to reduce local phosphorus inputs in stratified lakes to make them more resistant to less frequent mixing events under climate change, thereby revealing a twofold protecting benefit in addition to reducing harmful algal blooms’ susceptibility²⁷. Furthermore, conducting global scale lake models without taking into consideration changes in nutrient dynamics is likely to result in flawed predictions by disregarding the role of internal loads and the direct effects of management actions on these processes.

We elucidated shifts in the oxygen regime by comparing oxygen dynamics during
(1) natural conditions, i.e., pre-1950, (2) nutrient-altered conditions, i.e., 1950 to 2020, and
(3) climate-altered conditions, i.e., 2050 to 2100, while taking in-lake TP concentrations to be
< 24 µg L^–1. Longer periods of hypoxia (> 240 d yr^–1) and lower annual minimum DO concentration (< 1.3 mg L^–1) throughout the water column for 1950–2020 in comparison to
1860–1950 resulting from nutrient-enrichment in the lake (Fig. S11a) projected to remain as long-lasting regime shifts extending through the end of the century (Fig. S11b) reveal a transition from nutrient- to climate-altered conditions, even in the presence of TP < 24 µg L^–1. Minimum DO concentrations higher than 7.3 mg L^–1 still occurred during nutrient-altered condition, but were absent during climate-altered condition (Fig. S11b), revealing lake resilience to phosphorus input, but not to climate warming. Nutrient over-enrichment acted as a trigger to the establishment of hypoxia so that the ecosystem shifted to a new physico-chemical state at which point climate warming becomes as a new key control of lake oxygen dynamics (Fig. S11c). The fact that we were unable to detect any recovery from hypoxia by 2100, the end of our model run, may mean that climate forcing causes the eutrophication-driven transition to hypoxia to delay recovery so that it takes on a centennial scale²⁸.

The model projections presented here should be interpreted with some caution as the confidence in future forcing remains of uncertain reliability. For instance, changes in food web composition supported by biological invasion (e.g. quagga mussel populations first detected in Lake Geneva in the 2010s that are expected to increase in European alpine lakes in the next decades²⁹) would drastically affect organic matter and nutrient cycling and indirectly critical characteristics of the lake such as transparency, i.e. key parameters in our model, that could induce diverging trajectories from the modeled ones. Nevertheless, the combination of modelling, paleolimnological proxies, and field measurements provide strong insights into potential future changes in the hypoxia regime. Our overall results indicate a path for lake management and restoration programs which must consider that reducing local phosphorus inputs will likely not be sufficient to prevent hypolimnetic deoxygenation, thus demanding the implementation of other strategies or engineering approaches to prevent hypoxia. A non-trivial challenge for 21st-century environmental management will be to systematically counteract hypoxia. From a global perspective, lakes have undergone abrupt changes since the 1950s due to multiple kinds of environmental forcing worldwide⁶. As Lake Geneva shares similar trophic history with other lakes in Europe and North America, our results suggest a common trajectory for lake recovery and challenge for lake management.

Lake Geneva

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 km², its volume is 89 km³, 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 1950s²⁰. 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 period¹¹. And like other lakes, Lake Geneva faces the ongoing effects of climate change with significant warming trends and strengthening thermal stratification²⁵.

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 database¹⁶. 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 Rhône River with an average discharge of 186.4 m³ s^–1 accounting for 75% of the lake’s total inflow volume. The Rhône River is monitored by the national long-term surveillance of Swiss rivers¹³ 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 Rhône River has been measured monthly since 1919.

Paleolimnological data

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 lakes¹⁴. Daphnia has been successfully used as a quantitative indicator of lake trophic level due to its strong stoichiometric requirements for phosphorus³⁰. 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 sediments¹⁷. 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

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 (ISIMIP3b¹²). 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 bias²⁵. 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 observations³¹ 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 Geneva²⁰ 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 processes²⁵.

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 Rhône River as a single inflow into the lake³² 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 applications³³ 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 worldwide^34–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 properties³⁹. 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 resolution¹⁹. 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 lakes¹¹. 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 climate⁴⁰.

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)²⁵. 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 column⁴¹. 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 dynamics⁴². 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 GLM⁴¹ and AED2⁴².

Model configuration

The GLM-AED2 model (GLM: v.3.1.1, AED2: v.2.0) was run on an hourly time step from January 1^st 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 validated²⁵ was adopted and all AED2 parameters were set at their default values⁴². 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 dataset⁴³ in accordance with other aquatic ecosystem modeling studies^33,36,44, based on a comparison against field measurements collected in the deepest point of the lake over 64 years (1957–2020¹⁶; 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-a⁴⁵), 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 scales¹⁰ 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 uncertainties⁴⁶.

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⁻¹, and found a better match with the paleolimnological proxy at 1 mg L⁻¹, 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 ºC⁴⁰ to 3 ºC²⁰ 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^–1and a step size of 5 µg L^–1between 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 km²) 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 1^st January 2020 were set as simulated values in 31^st 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 dynamics²⁹.

The areal hypolimnetic mineralization rate (AHM, [g O₂ m⁻² d⁻¹]) 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 stratification²³.

Forel, F. A. Le Léman: Monographie Limnologique. Tome 2. Méchanique, Chimie, Thermique, Optique, Acoustique. (1895).
Wilson, J., Ucharm, G. & Beman, J. M. Climatic, physical, and biogeochemical changes drive rapid oxygen loss and recovery in a marine ecosystem. Sci Rep 9, (2019).
Jenny, J. P. et al. Scientists’ Warning to Humanity: Rapid degradation of the world’s large lakes. J Great Lakes Res 46, 686–702 (2020).
Grant, L. et al. Attribution of global lake systems change to anthropogenic forcing. Nat Geosci 14, 849–854 (2021).
McCarthy, F. M. G. et al. The varved succession of Crawford Lake, Milton, Ontario, Canada as a candidate Global boundary Stratotype Section and Point for the Anthropocene series. Anthropocene Review (2023) doi:10.1177/20530196221149281.
Huang, S., Zhang, K., Lin, Q., Liu, J. B. & Shen, J. Abrupt ecological shifts of lakes during the Anthropocene. Earth-Science Reviews vol. 227 Preprint at https://doi.org/10.1016/j.earscirev.2022.103981 (2022).
Deyle, E. R. et al. A hybrid empirical and parametric approach for managing ecosystem complexity: Water quality in Lake Geneva under nonstationary futures. (2022) doi:10.1073/pnas.
Jane, S. F. et al. Widespread deoxygenation of temperate lakes. Nature 594, 66–70 (2021).
Jenny, J. P. et al. Global spread of hypoxia in freshwater ecosystems during the last three centuries is caused by rising local human pressure. Glob Chang Biol 22, 1481–1489 (2016).
Soares, L. M. V. & Calijuri, M. do C. Deterministic modelling of freshwater lakes and reservoirs: Current trends and recent progress. Environmental Modelling and Software vol. 144 Preprint at https://doi.org/10.1016/j.envsoft.2021.105143 (2021).
Dresti, C., Fenocchi, A. & Copetti, D. Modelling physical and ecological processes in medium-to-large deep European perialpine lakes: A review. J Limnol 80, (2021).
Golub, M. et al. A framework for ensemble modelling of climate change impacts on lakes worldwide: The ISIMIP Lake Sector. Geosci Model Dev 15, 4597–4623 (2022).
Eawag & FOEN. NADUF - National long-term surveillance of Swiss rivers (2022-1) (Version 2022-1) [Data set]. (2022) doi:https://doi.org/10.25678/00069G.
Berthon, V., Alric, B., Rimet, F. & Perga, M. E. Sensitivity and responses of diatoms to climate warming in lakes heavily influenced by humans. Freshw Biol 59, 1755–1767 (2014).
Stockwell, J., Brentrup, J., Bruel, R., Doubek, J. & Anneville, O. Lake Geneva hypsometry data. Global Evaluation of the Impacts of Storms on freshwater Habitat and structure of phytoplankton Assemblages (GEISHA) Working Group, FRB Centre for the Synthesis and Analysis of Biodiversity and USGS John Wesley Powell Center. (2021).
Rimet, F. et al. The Observatory on LAkes (OLA) database: Sixty years of environmental data accessible to the public. doi:10.4081/jlimnol.2020.1944ï.
Jenny, J. P. et al. Inherited hypoxia: A new challenge for reoligotrophicated lakes under global warming. Global Biogeochem Cycles 28, 1413–1423 (2014).
Carlson, R. E. A trophic state index for lakes. Limnol Oceanogr 22, 361–369 (1977).
Gaillard, R., Perroud, M., Goyette, S. & Kasparian, J. Multi-column modelling of Lake Geneva for climate applications. Sci Rep 12, (2022).
Schwefel, R., Gaudard, A., Wüest, A. & Bouffard, D. Effects of climate change on deepwater oxygen and winter mixing in a deep lake (Lake Geneva): Comparing observational findings and modeling. Water Resour Res 52, 8811–8826 (2016).
Woolway, R. I. & Merchant, C. J. Worldwide alteration of lake mixing regimes in response to climate change. Nat Geosci 12, 271–276 (2019).
Vautard, R., Cattiaux, J., Yiou, P., Thépaut, J. N. & Ciais, P. Northern Hemisphere atmospheric stilling partly attributed to an increase in surface roughness. Nat Geosci 3, 756–761 (2010).
Müller, B. et al. Oxygen consumption in seasonally stratified lakes decreases only below a marginal phosphorus threshold. Sci Rep 9, (2019).
Woolway, R. I. et al. Phenological shifts in lake stratification under climate change. Nat Commun 12, (2021).
Desgué-Itier, O. et al. Past and future climate change effects on the thermal regime and oxygen solubility of four peri-alpine lakes. Hydrol Earth Syst Sci 27, 837–859 (2023).
Shatwell, T., Thiery, W. & Kirillin, G. Future projections of temperature and mixing regime of European temperate lakes. Hydrol Earth Syst Sci 23, 1533–1551 (2019).
Kakouei, K. et al. Phytoplankton and cyanobacteria abundances in mid-21st century lakes depend strongly on future land use and climate projections. Glob Chang Biol 27, 6409–6422 (2021).
Haas, M. et al. Roman-driven cultural eutrophication of Lake Murten, Switzerland. Earth Planet Sci Lett 505, 110–117 (2019).
Kraemer, B. M. et al. Title: An abundant future for quagga mussels in deep European lakes 1. doi:10.1101/2023.05.31.543086.

Methods References

Sterner, R. W. & Elser, J. J. Ecological Stoichiometry: The Biology of Elements from Molecules to the Biosphere. (Princeton University Press, 2002).
MétéoSuisse Data Warehouse. Data management. https://www.meteoswiss.admin.ch (2023).
Bueche, T., Hamilton, D. P. & Vetter, M. Using the General Lake Model (GLM) to simulate water temperatures and ice cover of a medium-sized lake: a case study of Lake Ammersee, Germany. Environ Earth Sci 76, 461 (2017).
Ladwig, R. et al. Lake thermal structure drives interannual variability in summer anoxia dynamics in a eutrophic lake over 37 years. Hydrol Earth Syst Sci 25, 1009–1032 (2021).
Farrell, K. J. et al. Ecosystem-scale nutrient cycling responses to increasing air temperatures vary with lake trophic state. Ecol Modell 430, (2020).
Weng, W. et al. Coupling Natural and Human Models in the Context of a Lake Ecosystem: Lake Mendota, Wisconsin, USA. Ecological Economics 169, 106556 (2020).
Fenocchi, A. et al. Applicability of a one-dimensional coupled ecological-hydrodynamic numerical model to future projections in a very deep large lake (Lake Maggiore, Northern Italy/Southern Switzerland). Ecol Modell 392, 38–51 (2019).
Ladwig, R. et al. Lake thermal structure drives inter-annual variability in summer anoxia dynamics in a eutrophic lake over 37 years. Hydrology and Earth System Sciences Discussions 1–45 (2020) doi:10.5194/hess-2020-349.
Bucak, T. et al. Modeling the effects of climatic and land use changes on phytoplankton and water quality of the largest Turkish freshwater lake: Lake Beyşehir. Science of the Total Environment 621, 802–816 (2018).
Bruce, L. C. et al. A multi-lake comparative analysis of the General Lake Model (GLM): Stress-testing across a global observatory network. Environmental Modelling and Software 102, 274–291 (2018).
Perroud, M., Goyette, S., Martynov, A., Beniston, M. & Anneville, O. Simulation of multiannual thermal profiles in deep Lake Geneva: A comparison of one-dimensional lake models. Limnol Oceanogr 54, 1574–1594 (2009).
Hipsey, M. R. et al. A General Lake Model (GLM 3.0) for linking with high-frequency sensor data from the Global Lake Ecological Observatory Network (GLEON). Geosci Model Dev 12, 473–523 (2019).
Hipsey, M. R. Modelling Aquatic Eco-Dynamics: Overview of the AED modular simulation platform. Preprint at (2022).
Hipsey, M. R. et al. A system of metrics for the assessment and improvement of aquatic ecosystem models. Environmental Modelling and Software 128, (2020).
Mi, C. et al. Ensemble warming projections in Germany’s largest drinking water reservoir and potential adaptation strategies. Science of the Total Environment 748, 141366 (2020).
Soulignac, F. et al. Contribution of 3D coupled hydrodynamic-ecological modeling to assess the representativeness of a sampling protocol for lake water quality assessment. Knowl Manag Aquat Ecosyst 42 (2019) doi:10.1051/kmae/2019034.
John, A., Nathan, R., Horne, A., Stewardson, M. & Angus Webb, J. How to incorporate climate change into modelling environmental water outcomes: A review. Journal of Water and Climate Change vol. 11 327–340 Preprint at https://doi.org/10.2166/wcc.2020.263 (2020).

There is NO Competing Interest.

Supplementarymaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

A crisis of lake hypoxia in the Anthropocene: The long-term effects of climate and nutrients

Status:

Version 1

Abstract

Figures

Background

General approach

Historical lake hypoxia

Drivers of historical hypoxia

Hypoxia under future climate change

Methods

References

Additional Declarations

Supplementary Files

Status:

Version 1