Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.
Advertisement
Nature Climate Change (2024)
Metrics details
Oxygen depletion constitutes a major threat to lake ecosystems and the services they provide. Most of the world’s lakes are located >45° N, where accelerated climate warming and elevated carbon loads might severely increase the risk of hypoxia, but this has not been systematically examined. Here analysis of 2.6 million water quality observations from 8,288 lakes shows that between 1960 and 2022, most northern lakes experienced rapid deoxygenation strongly linked to climate-driven prolongation of summer stratification. Oxygen levels deteriorated most in small lakes (<10 ha) owing to their greater volumetric oxygen demand and surface warming rates, while the largest lakes gained oxygen under minimal stratification changes and improved aeration at spring overturns. Seasonal oxygen consumption rates declined, despite widespread browning. Proliferating anoxia enhanced seasonal internal loading of C, P and N but depleted P long-term, indicating that deoxygenation can exhaust redox-sensitive fractions of sediment nutrient reservoirs.
The availability of dissolved oxygen (DO) is fundamental to the health and functioning of lake ecosystems but an alarming number of lakes experience rapid deoxygenation today1,2,3. Oxygen is an essential driver of the biogeochemical processes that regulate nutrient cycling, drinking water quality, biodiversity and aquatic food webs4. Aerobic respiration can deplete oxygen concentrations at depth when seasonal stratification prevents aeration of the water column5. This detrimental condition, termed hypoxia (DO < 2 mg l−1), diminishes the cold-water habitat of aquatic life, including commercially important fish6 and can cause mass mortality if persistent7. Anoxia (DO < 0.5 mg l−1) promotes the reductive dissolution of mineral-bound carbon and phosphorus8,9 which fuels harmful algal blooms10 and limits microbial oxidation of methane, a potent greenhouse gas11. The negative environmental consequences of deoxygenation ultimately compromise human livelihoods and health12.
Historically, oxygen depletion in lakes has been associated with excessive inputs of nutrients (eutrophication) from urban and agricultural sources, which increase oxygen demand via enhanced primary production and decomposition of sinking biomass1,13. Elevated loading of terrestrial organic matter (browning) linked to soil recovery from acidification, climate-driven shifts in precipitation patterns and land use change14,15,16, could similarly contribute to deoxygenation. Climate warming has the potential to cause further oxygen loss via prolonged summer stratification17,18 and enhanced aerobic respiration rates19. Seasonal hypoxia tends to develop more rapidly in stratified lakes with high nutrient concentrations20 and a large sediment surface area relative to the overlying water volume (Ased/Vwater), as benthic respiration is often the dominant oxygen sink5,21. This implies that oxygen levels in specific lake classes, such as small, eutrophic or humic lakes, may be more sensitive to climate forcing but whether predictors of seasonal hypoxia extend to long-term deoxygenation is unclear.
Most of the world’s lakes are located in boreal and arctic regions22, where the climate is warming at two to four times the global rate23. Rising air temperatures, atmospheric stilling and earlier ice loss have prolonged summer stratification in lakes across the region24,25,26 and, combined with widespread browning15,27, put pristine freshwater ecosystems at increasing risk of hypoxia28. At the same time, declining phosphorus and nitrogen stocks in northern landscapes associated with reduced terrestrial loading and atmospheric deposition29,30,31,32 could limit fresh organic substrate fuelling oxygen demand. Moreover, overall shortening of the ice-cover season33 may prolong spring mixing and enhance aeration before summer stratification34,35. The rapid accumulation of counteractive anthropogenic pressures means that the sign and magnitude of long-term bottom-water DO trends in northern lakes remains unknown but a continental-scale assessment, as was successfully done in 429 temperate lakes3,18, is lacking. Our inability to resolve the global extent and causes of deoxygenation hampers international commitments to address this dimension of climate impacts and protect freshwater quality.
In this study, we conducted a large-scale systematic analysis of long-term water quality monitoring data from northern lakes, comprising 2.6 million observations from 8,288 lakes in Fennoscandia and North America from 1960 to 202236. A map of the sampled lakes is shown in Extended Data Fig. 1. We compiled, harmonized and quality-checked measurements of DO, temperature, total phosphorus (TP), total nitrogen (TN), colour (a measure of chromophoric dissolved organic matter, CDOM) and lake morphometry from governmental databases and non-profit organizations via data portals and data requests. The goals of our analyses were to: (1) identify long-term changes in the bottom-water oxygen content of northern lakes, (2) assess whether standard predictors of seasonal hypoxia (surface area, maximum depth and trophic state) apply to long-term DO trends and (3) evaluate potential mechanisms driving long-term (de)oxygenation, including persistent shifts in stratification phenology, oxygen demand and spring aeration.
Between 1960 and 2022, ice-free mean oxygen concentrations in lake bottom waters (≥0.7 × maximum depth) have decreased substantially across monitored northern lakes (Fig. 1). In ‘trend lakes’ with ≥15 year time series, the occurrence of bottom-water anoxia has increased from 39% to 61% between the first and last 5 years of measurement (n = 912). Oxygen trends scaled inversely with lake surface area from −0.40 to +0.07 mg l−1 decade−1 in lakes sized ≤10 ha and >104 ha, respectively (Fig. 2a). This scaling pattern was robust across continents (Fennoscandia and North America) and biomes (tundra, boreal forests, north-temperate forests and grasslands)37 (Extended Data Fig. 2). Because small lakes and ponds are disproportionally abundant22, our results imply that most northern lakes have experienced major deoxygenation. In our dataset, small lakes are under-represented compared to their real relative abundance. Adjusting for this bias via weighted resampling, we estimated that bottom-water DO has decreased below the critical hypoxic threshold (2 mg l−1) in over half of northern lakes (Fig. 1a, black line).
a–d, Bottom-water oxygen concentrations classified by lake surface area (a), maximum depth (b), trophic state (c) and humic content (d) across 7,394–8,288 lakes, depending on data availability. Annual class-medians (dots) were computed from ice-free season lake–year mean DO corrected for bias toward mid-season sampling (Methods). GAM smooths were weighted by inverse variance to account for uncertainty at low sample counts. Shaded areas represent 95% confidence intervals (CIs). The black curve in a shows a weighted population-median estimate based on the size distribution of Swedish lakes (Methods). In a and b, legend labels define class boundaries (for example in a, green marks Asurf = 102–103 ha). Trophic state categories are oligotrophic (O), mesotrophic (M) and eutrophic (E). Humic content categories are clear (C), humic (H) and polyhumic (PH). Dashed and dotted lines mark hypoxia (DO < 2 mg l−1) and anoxia (DO < 0.5 mg l−1), respectively.
a–c, Box plots of long-term trends (≥15 years) of bottom-water DO (a), stratification strength (Δρ) (b) and sediment surface area to water volume ratio (Ased/Vwater) (c) by logarithmic lake surface area class. d, GAM residuals of DO trends as a function of Ased/Vwater, Δρ trends and sampling day-of-year trends. Boxes bound the interquartile range (IQR) (25th to 75th percentile) and whiskers extend to the smaller quantity of data extremes or medians ± 1.5× IQR. Dots and lines mark means and medians, respectively. In a, b and d, open symbols signify means whose 95% CIs include 0. Outliers are shown by grey dots. Distinct letters (top) mark statistically different groups (Padj < 0.05, Conover–Iman tests). Bottom labels show the sample size (number of trend lakes) in each surface area class, which differed per variable depending on data availability.
The areal scaling of long-term bottom-water oxygen trends (Fig. 1a) divided into four distinct size groups (Conover–Iman test, Padj < 0.05) at 102, 103 and 104 ha (trend lakes, Fig. 2a), with significant downward trends in lakes <103 ha and non-significant mean increases in the largest lakes (>104 ha). Among lakes classified by their maximum depth (zmax, 20th percentiles), those of intermediate zmax (4–31 m) showed substantial oxygen declines, whereas bottom waters of shallow (zmax < 4 m) and deep (zmax ≥ 31 m) lakes showed no clear long-term trend and generally remained oxygenated throughout the ice-free season (Fig. 1b). In addition, long-term oxygen loss was markedly more pronounced in eutrophic and polyhumic (highly coloured) lakes compared to less productive and clearer water bodies (Fig. 1c,d). These covariates are consistent with empirical models of seasonal oxygen loss as a function of basin morphometry and nutrient concentrations20 and demonstrate that similar emergent mechanisms modulate long-term oxygen trends.
To identify the dominant drivers of the observed oxygen trends we conducted random forest (RF) analyses38 of ice-free mean bottom-water DO saturation (DOsat) and its long-term trend. Separating the integration timescale is necessary because cross-scale interactions among drivers identified in seasonal and comparative studies may shift their importance on multiyear scales39. Predictor variables included stratification strength (Δρ), sampling day-of-year, temperature, colour, TP, TN and their respective long-term trends and Ased/Vwater. We used surface-water colour and nutrients because in the hypolimnion, anoxia-induced internal loading confounds causal relations with DO8. The RF models explained 76% and 28% of the variance of DOsat and DOsat trends, respectively. Predictor variable importance scores ranked stratification strength as the primary covariate of seasonal and long-term DO depletion, followed by Ased/Vwater (Extended Data Fig. 3). Partial dependence plots show that bottom-water DOsat tended to decrease with increasing Δρ, Ased/Vwater and surface-water TP and TN (Fig. 3a). DOsat trends were more negative (deoxygenation) in lakes experiencing strong positive Δρ trends and surface-water browning and eutrophication (Fig. 3b). To assess changes in oxygen consumption, we computed the seasonal bottom-water oxygen depletion rate (∂DO/∂t) from annual sets with more than four observations of decreasing DO concentrations at each profiling location under stratified conditions (Δρ ≥ 0.1 kg m−3) (n = 2,528). Assuming negligible transport of DO through the thermocline40, this rate closely approximates biogenic volumetric oxygen demand5. We found that ∂DO/∂t trended negative in 64% of a small but representative set of trend lakes (−1.99 ± 1.31 µg−1 l−1 d−1 decade−1, n = 53), suggesting that DO consumption rates have not systematically increased in northern lakes and therefore do not generally explain the trend towards hypoxia.
a,b, RF partial dependence plots of ice-free mean bottom-water DOsat (a) and its long-term trend (b). a shows the six most important predictor variables identified in the RF model; Δρ, Ased/Vwater, sampling day-of-year and surface-water colour, TP and TN. b includes the top three predictors of the DOsat trend (Δρ trend, Ased/Vwater, sampling day-of-year trend) as well as trends in surface-water colour, TP and TN. Bars summarize statistics of the independent variables and follow the legend of Fig. 2. Numbers rank the permutation importance z-score relative to that of the most important feature. The y axes are the same scale for each subpanel; grey dashed lines mark the intercepts. The complete importance ranking of covariates is shown in Extended Data Fig. 3.
We evaluated long-term trends in bottom-water oxygen during the spring overturn and the start of summer stratification (first 3 weeks after ice-off). Size-class-mean concentrations varied between 5.6 and 10.6 mg l−1 (48%–84% saturation) and increased with lake size (Fig. 4). Contrary to summer observations, initial oxygen trended positive in all but the smallest lakes. To examine how these partly opposing trends are linked, we estimated decadal changes in seasonal DO cycles by fitting a hierarchical generalized additive model (HGAM)41, using a normalized time variable—day-of-year scaled to mean ice-on and ice-off dates—to account for geographical differences in the length of the ice-free season (Methods). The HGAM revealed that long-term trends toward oxygenation in spring can offset deoxygenation from earlier onset of stratification (Fig. 5) and that the balance between the two processes is, in part, a function of lake morphometry. Enhanced spring aeration elevated summer DO in some of the largest lakes with negligible shifts in stratification phenology (≥104 ha) but the net offsetting effect diminished with increasing thermal stability trends toward mid-sized lakes (Figs. 2b and 5).
Annual median bottom-water oxygen concentrations in the first 3 weeks after ice-off by logarithmic lake surface area class. GAM smooths were weighted by inverse variance to account for uncertainty at low sample counts. Shaded areas represent 95% CIs. Legend labels define area class boundaries (for example, green marks Asurf = 102–103 ha). Dashed and dotted lines mark hypoxia (DO < 2 mg l−1) and anoxia (DO < 0.5 mg l−1), respectively.
HGAM estimates of mean bottom-water DO as a function of (scaled) day-of-year, year, lake surface area (Asurf), maximum depth (zmax) and relative sampling depth (zsample/zmax) (Methods). Predictions are shown for years 1970 and 2020 at a depth of 0.85 × zmax in order-of-magnitude lake surface area classes and class-median zmax (that is, for Asurf = 10 ha, the median zmax of lakes 5–50 ha). Day-of-year was scaled to lake-mean ice-on and ice-off derived from the ERA5-Land ‘lake ice depth’ product61 to account for geographic differences in ice-free season length. Vertical dashed lines mark local oxygen maxima and approximate the timing of vernal and autumnal lake mixing. Shaded areas cover 95% Bayesian CIs.
We examined whether any one mechanism dominates the morphometric scaling of DO trends (Fig. 2a). First, trends in stratification strength (Fig. 2b) and duration (Fig. 5) increased with decreasing lake size. Second, the predictive power of Ased/Vwater in the RF analyses on both seasonal and interannual timescales (Fig. 3) suggests that the response to climate perturbations (for example, an additional day of stratification per year) is proportional to areal extent of the sediment oxygen sink (Ased/Vwater), which scales inversely with lake surface area (Fig. 2c). Third, DO at the start of the ice-free season trended positive in large lakes (Figs. 4 and 5). We constructed a generalized additive model (GAM) of DO trends with terms for Δρ trends (1), ln(Ased/Vwater) (2) and ice-free mean sampling day-of-year trends (3), plus an interaction term (1) × (2) for lakes where both DO and Δρ trends could be computed (n = 787 lakes). Term (3) accounted for progressively later sampling in some lakes which could have biased DO trends low (Fig. 3b) but was not significant (P = 0.112). Figure 2d presents summary statistics of the model residuals. The model explained class-mean differences between DO trends (Kruskal–Wallis test on model residuals, P = 0.19, d.f. = 4) (Fig. 2d), while significant differences remained in the residuals of GAMs excluding either term (1) (P = 0.0071, d.f. = 4) or term (2) (P < 0.0001, d.f. = 4). Spring DO trends were omitted from the GAM for lack of trend lakes (n = 84) but accounted for the non-significant, positive mean DO trend in lakes >104 ha (Fig. 2a) remaining in the model residuals (Fig. 2d). Thus, the lake-size dependency of bottom-water DO trends emerged from a synergy between biotic and physical drivers.
Seasonal anoxia is known to enhance internal loading of CDOM and P (via reductive desorption) and N (as NH4 via nitrification) in the hypolimnion of stratified lakes but the long-term effects of the observed proliferation of anoxia (Fig. 1) remain unknown. We computed seasonal rates of change of bottom-water colour, TP, TN and NH4 under distinct oxygen regimes for every lake–year with four or more observations in the stratified period (Δρ ≥ 0.1 kg m−3). As expected, we found significant accumulation of all four substances under anoxic conditions (Fig. 6). Rate magnitudes decreased with increasing DO and were negligible in highly oxic (>8 mg l−1) waters. The seasonal increase in TN under anoxic conditions was accounted for by NH4. We next computed bottom-water colour, TP, TN and NH4 trends (≥15 years) under the same oxygen regimes to assess their long-term impacts. Overall, mean bottom-water colour trends were positive, TP trends were negative and TN and NH4 trends tended to zero (Fig. 6). The magnitudes of colour trends increased with decreasing DO and were approximately four times higher in anoxic compared to oxic hypolimnia. In contrast, and perhaps counterintuitively, we found that recurring anoxia was associated with dispersed trends in bottom-water TN and NH4 and significant downward trends in TP.
a–d, Estimates of ice-free seasonal rates (within years, orange) and long-term trends (across years, cyan) of bottom-water colour (a), TN (b), TP (c) and NH4 (d) within distinct oxygen regimes: anoxic (≤0.5 mg l−1), hypoxic (0.5–2 mg l−1), low oxic (2–5 mg l−1), mid oxic (5–8 mg l−1), high oxic (>8 mg l−1) and ‘All’ (no oxygen criterion). Long-term surface-water trends are shown for comparison as open symbols. In the y axis labels, f stands for frequency (d−1 and yr−1). The category ‘All’ includes all bottom-water samples. Error bars mark 95% CIs of the means. Statistics associated with this figure are summarized in Extended Data Table 1.
Long-term deoxygenation of northern lakes was strongly coupled to increased density stratification (Δρ) (Fig. 3), which in turn is driven by climate warming18,24,25. While long-term patterns in air temperature and wind speed were near-identical across lake-size classes (Extended Data Fig. 4), smaller lakes experienced the greatest prolongation of summer stratification (HGAM, Fig. 5) and Δρ trends (Fig. 2b). This suggests that lake morphometry mediates thermal responses to climate warming. This notion is corroborated by model and field studies; in small dimictic lakes with limited fetch, reduced wind mixing confines heating to a shallow surface layer, increasing the thermal gradient42, whereas in lakes with a large surface area, wind mixing distributes heat over greater depths43 and interannual variation of stratification duration is governed by wind speed more than air temperature44.
Our results further indicate that prolonged thermal stratification, rather than fewer mixing events, is the main driver of long-term oxygen loss. Local oxygen maxima in the modelled seasonal DO cycles (Fig. 5) revealed that spring and autumn mixing occurred significantly earlier and later in the year, respectively, consistent with reports of widespread shortening of lake ice-cover periods across the northern hemisphere33. Our findings are also in agreement with previous studies in temperate lakes which report coupled increases in stratification duration and bottom-water hypoxia17,18,34,45 and, more rarely, negligible DO trends where the stratification period has not lengthened but shifted46. Here, dividing long-term DO trends by seasonal ∂DO/∂t in lakes where both estimates were available (n = 251), we estimate that an extra (median) 1.7 days of DO drawdown per decade accounted for the DO trends. This is comparable to model estimates of historic (1960–2005) prolongation of summer stratification in northern hemisphere lakes larger than 100 ha (ref. 24), about 2 d decade−1.
In addition to duration of stratification, the degree of aeration at the start of the stratified period in spring can severely impact the trajectory of oxygen loss over summer5,35. We found that initial (spring) oxygen concentrations decreased predictably with lake surface area (Figs. 4 and 5). This pattern may in part be a remnant of winter oxygen levels, which, as in summer, were elevated in larger lakes47. Additionally, during the brief spring overturn in smaller lakes, their smaller fetch can limit efficient air–water gas exchange (aeration)48. The time series of initial DO (Fig. 4) and the HGAM (Fig. 5) further showed that, in spring, DO trended positive in lakes larger than ~103 ha. Such improved bottom-water aeration could indicate that earlier ice-off has caused a lengthening of the spring mixing period35 but the drivers of this widespread and consistent pattern remain unclear. We hypothesize that in large lakes with substantial thermal inertia, excess warming in winter26 has had a greater effect on the timing of ice-off than on stratification onset.
Previous studies have linked hypolimnetic oxygen loss to cultural eutrophication of surface waters and increased availability of organic matter1,10,49. In agreement with these reports, we found that long-term bottom-water oxygen trends scaled inversely with trends in surface-water TP and colour (Fig. 3b). Here, oligotrophication of northern lakes29,30,31 may have offset some of the climate-driven deoxygenation, which would be consistent with the observed diminishing of seasonal oxygen depletion rates (∂DO/∂t). Continental-scale patterns in surface-water browning14,15,27 may have exacerbated long-term bottom-water oxygen loss via enhanced heat absorption at the water surface strengthening thermal stratification50. Browning may also have increased lacustrine degradation of terrestrially derived organic matter51, although this seems inconsistent with observed ∂DO/∂t declines. Overall, biogeochemical predictors ranked consistently lower in importance than physical quantities Δρ and Ased/Vwater (Extended Data Fig. 3). Compared to climate warming, browning and nutrient trends appear to have had a limited effect on oxygen loss in northern lakes at a monitoring timescale of ∼25 years (trend lake median). Another way to interpret this result is that biologically mediated changes in lake metabolism tend to occur at a much slower pace1.
Conversely, in an increasing number of lakes experiencing annual anoxia (Fig. 1), biogeochemical cycling of carbon and nutrients may be affected long-term. The positive CDOM trends in anoxic waters (Fig. 6) imply that anoxic liberation of buried organic matter, previously shown in sediment incubations and in situ in lakes on seasonal timescales8,52,53, is not only widespread but has also increased over the last six decades. As a result, bottom-water colour trends were significantly greater than landscape-driven browning rates measured at the surface (trend lakes, pairwise one-sided Wilcoxon signed rank test, P ≪ 0.0001, n = 678), with a factor two difference between mean trends (Fig. 6). Notably, a recent study found that dissoved organic matter (DOM) ‘deprotected’ from iron complexes is highly reactive and that up to 21% can escape coprecipitation with iron upon re-aeration54. If retained in the lake, this labile DOM pool could, together with reduced substances that also accumulate in anoxic hypolimnia (CH4, NH4, S(−II), Mn(II), Fe(II))55, deplete the limited supply of epilimnetic DO mixed down during the autumn overturn56 and reinforce anoxic conditions over longer timescales. Latent oxygen sinks may be particularly important in smaller lakes with higher Ased/Vwater (Fig. 2c) and greater build-up of sediment-derived substances55. Here, the HGAM (Fig. 3) indicates that bottom-water oxygen concentrations at ice-on have declined in lakes smaller than 100 ha, but more work is needed to identify the cause(s) of this trend.
As with CDOM, we might expect prolonged stratification (and anoxia) to elevate TP concentrations long-term through increased internal loading. While bottom-water TP trends were generally more positive than surface-water TP trends (pairwise one-sided Wilcoxon signed rank test, P = 0.003, n = 820), this pattern reversed under persistently recurring anoxia (≥15 years) (Fig. 6c). We hypothesize that the strong downward TP trend in anoxic bottom waters was caused by the gradual depletion of iron-bound P reservoirs in the thin, periodically oxygenated layer of the surface sediment. One mechanism that could drive such depletion is incomplete coprecipitation of P with Fe upon re-aeration57, if the residual P is exported out of the lake or sequestered into less redox-sensitive complexes. Prolonged anoxia evidently did not exhaust the redox-sensitive organic matter fraction in surface sediments (Fig. 6a), probably because its concentration tends to be several orders of magnitude higher than that of iron-bound P58,59.
Northern lakes, which have until recently been less affected by human development, are now at severe risk of prolonged hypoxia with consequent ecological deterioration. Recent long-term studies in north-temperate lakes show the rapid diminishing of the oxythermal habitat of commercially important fish under concurrent browning and warming60 and the potential for mass kills during heat waves7. On the basis of our results, we expect a similar, ongoing impact in boreal and arctic regions, where cold-water fisheries fulfil important cultural and economic functions. Moreover, the proliferation of anoxia could have major biogeochemical implications, including increased remobilization of terrestrially derived DOM and a decrease in the amount of terrestrial carbon that is ultimately buried in lake sediments53,54. Deoxygenation has particularly affected small lakes, globally the most abundant lake type and is a widespread problem with important socio-economic consequences that demands urgent attention.
In this section, we describe how we sourced the lake water quality and morphometry data, quality-controlled the observations, computed derived quantities, such as volumetric oxygen demand and analysed the data with established statistical methods in R v.4.2.0 (ref. 62). Variables were selected on the basis of known (proxy) relations with oxygen and availability across the different datasets. The goal of the statistical analysis was to identify environmental covariates of bottom-water oxygen concentrations and their long-term trends. This approach enabled us to disentangle the impact of climate warming on deoxygenation from that of anthropogenic changes in lake metabolism associated with eutrophication and browning. A detailed description of these methods follows.
In this study, we focused on seasonally ice-covered lakes in the northern regions of North America (>44° N) and Northern Europe (>55° N), where most of the world’s lakes are located22 and the climate is warming significantly faster compared to the global average23. We compiled lake water quality measurements from a wide variety of governmental and not-for-profit organizations in the United States (Alaska), Canada (British Columbia, Québec, Ontario, Saskatchewan, Alberta, Manitoba, Nova Scotia, New Brunswick, Prince Edward Island and the Northwest Territories), Finland, Sweden and Norway (Supplementary Data 1). Data were accessed directly through online data portals (7,855 lakes) or obtained via data requests (433 lakes). Lake morphometry information was similarly collected from public repositories, scientific literature and data requests (Supplementary Data 2) and matched with the water quality data through common water body identifiers or manual colocation. Lakes without maximum depth or surface area information were excluded from the study. The synthesized dataset comprises 2.6 million measurements in 8,288 lakes36.
We included freshwater lakes and ponds in the dataset but we excluded reservoirs, non-permanent lakes and actively managed water bodies. We classified managed water bodies as those undergoing physical measures against deoxygenation, such as aeration, dredging or hypolimnetic withdrawal, as well as chemical treatments to combat eutrophication which are likely to affect oxygen consumption rates, such as the addition of aluminium sulfate. Because measures to improve water quality are often costly and disruptive they are generally well documented. We identified 197 managed lakes and ponds in our dataset with help from the data providers, through news media, reports and published research and via regional authorities and local conservation organizations. Basins of large lakes that were sampled separately were treated as distinct water bodies. Each lake was assigned a new identifier to replace originals which could overlap between datasets and to merge observations from water bodies sampled by multiple organizations, for example across borders.
Physical and chemical variables were selected on the basis of known (proxy) relations with bottom-water oxygen concentrations and availability in the different datasets. In our data synthesis, we included the following variables: water temperature, true water colour, TP, TN, ammonia and chlorophyll. TP and TN were computed as the sum of particulate and dissolved fractions (that is, unfiltered samples). When TN measurements were unavailable, such as in some profiles of the US WQP dataset63, we computed TN from total Kjeldahl nitrogen (ammonia + organic N), nitrate and nitrite. In analyses involving long-term TN trends, we confirmed that TN was computed consistently through time. We also included sampling coordinates, sampling date, sampling depth and lake altitude, surface area, mean depth and maximum depth in the dataset.
We used standard equations to compute freshwater density64 and DO saturation with respect to atmospheric concentrations65. Strength of thermal stratification was estimated as the difference between mean bottom- and surface-water density (Δρ, kg m−3). There exists no universal definition of stratification in lakes but density thresholds are frequently applied66 and we defined stratified profiles as those with Δρ exceeding 0.1 kg m−3(ref. 24). Salinity measurements were unavailable in most lakes and were not included in the computation of freshwater density. This means that in chemically stratified naturally hypersaline lakes, which can be found in select areas of the northern Great Plains in western Canada67, our RF estimate of the importance of Δρ (Extended Data Fig. 3a) is probably conservative. Trophic state classes were assigned on the basis of surface-water chlorophyll or, if unavailable, TP concentrations as oligotrophic (chlorophyll (Chl) ≤ 2.6 µg l−1, TP ≤ 12 µg l−1), mesotrophic (Chl = 2.6–7.3 µg l−1, TP = 12–24 µg l−1) and eutrophic (Chl > 7.3 µg l−1, TP > 24 µg l−1)68. Humic state was divided into three categories on the basis of surface-water colour, measured on a platinum–cobalt (Pt) colour scale; clear (<30 mg Pt l−1), humic (30–90 mg Pt l−1) and polyhumic (>90 mg Pt l−1)69. The Ased/Vwater ratio below 0.7 × zmax was estimated using hypsographs computed with the approx.bathy function in the rLakeAnalyzer R package70.
Quality control of the synthesis dataset was conducted as follows. First, we removed observations without a sampling date or (plausible) sampling depth (0–lake maximum depth). However, bathymetric surveys based on transects can be incomplete, so in water bodies where the sampling depth consistently (in more than five profiles) exceeded the reported zmax we adjusted zmax accordingly. We standardized units of measurement across the dataset and only kept measurements with unambiguous units, for example NH4, TN or TP in molar quantities (µM) or concentrations ‘as N’ and ‘as P’ and clearly defined fractions (‘total’ or ‘dissolved’)71. We subsequently removed water quality measurements outside a predefined plausibility range, generally excluding the top 0.1% quantile: DO (0–20 mg l−1), water temperature (0–40 °C), TP (0–3 mg l−1), TN (0–60 mg l−1), water colour (0–800 mg Pt l−1) and NH4 (0–8 mg l−1). The range filter was applied to the raw data before the computation of derived quantities (see the following sections). For all analyses, except the HGAM of seasonal DO cycles (Fig. 5), we removed profiles outside the ice-free period (see section on ‘Definition of the ice-free period’). In figures comparing class-mean values, trends or rates, we removed extreme outliers (>3 s.d. from the mean) which comprised 1.3%–2.2% of values.
We defined the bottom water of each lake as the layer below 0.7 × maximum depth and surface water as the layer above the smaller quantity of 0.5 × maximum depth and 1 m depth. A more traditional approach defines the surface and bottom layer as the epi- and hypo-limnion, respectively, but this would exclude lakes and lake-years without a clearly defined seasonal thermocline. By including periods when lakes are fully mixed, our definition does not mask the ongoing shift from polymictic to dimictic mixing regimes in many small and shallow lakes as a result of surface warming24.
We estimated ice phenology for each lake from the ‘lake ice depth’ product in the ERA5-Land reanalysis dataset61. The product is based on the ECMWF weather model coupled with the 1D hydrodynamic lake model FLake72. FLake has been shown to accurately predict ice-on and ice-off dates in north-temperate lakes with a high accuracy73. In each lake, the ice-free season was defined as the portion of the period between 1 June and 31 October when the ice thickness was 0. We also computed long-term (1960–2022) mean ice-on and ice-off day-of-year for each lake as temporal reference points for the HGAM.
We estimated the magnitude of the interannual trends of bottom-water oxygen concentrations and environmental variables. First, we computed ice-free, depth-averaged annual mean values for each variable and lake–year. In lakes where persistent anoxia biased the oxygen trend towards 0, rather than omitting the lakes from the trend analysis entirely, we removed each middle value of three consecutive years with ice-free mean DO values below 0.5 mg l−1. This procedure preserved initial deoxygenation preceding long-term anoxia in lakes with multidecadal DO time series. We then selected the subset of lakes with time series exceeding 15 years, comprising 460–951 trend lakes depending on the variable. Omitting anoxic years reduced the number of DO trend lakes from 912 to 805. For each lake and variable, we quantified long-term changes with the Theil–Sen slope estimator74, a non-parametric technique for trend estimation that is robust to outliers and non-normality.
We estimated ice-free season (within-year) seasonal oxygen depletion rate (∂DO/∂t) as the rate of change of bottom-water DO with the Theil–Sen estimator74 for time series with a minimum of four measurements per lake–year. We selected, for each lake and year, the part of the time series with decreasing DO concentrations and excluded non-stratified profiles (Δρ < 0.1 kg m−3). In this way, we excluded periods of intermittent and autumn mixing from the computation of ∂DO/∂t. The computed seasonal oxygen depletion rate reflects a balance of oxygen sources (diffusion and photosynthesis) and sinks (biochemical oxygen consumption). During summer stratification, transport of DO through the thermocline is typically orders of magnitude lower than hypolimnetic oxygen consumption40. In addition, photosynthesis is less likely to play a significant role as 87.4% of bottom-water oxygen samples were collected from the aphotic zone (<1% of surface photosynthetically active radiation at >2.7 × Secchi depth4). We therefore interpreted ∂DO/∂t to reflect primarily the volumetric oxygen demand of aquatic organisms. Even so, we cannot completely exclude an influence of oxygen sources and, as such, oxygen consumption rates inferred from ∂DO/∂t should be considered conservative.
To examine functional relations between oxygen and biochemical variables, we computed seasonal rates and long-term trends in bottom-water colour, TP, TN and NH4 under distinct DO regimes: anoxic (≤0.5 mg l−1), hypoxic (0.5–2 mg l−1), low oxic (2–5 mg l−1), mid oxic (5–8 mg l−1) and high oxic (>8 mg l−1), plus an additional class that included all DO concentrations. We excluded non-stratified profiles (Δρ < 0.1 kg m−3) to minimize the influence of turbulent mixing. For each variable, we estimated seasonal rates by first computing mean values for each sampling date and within each available oxygen regime and then estimating a rate for each lake–year and oxygen class with more than four observations. Long-term trends were similarly computed from mean values across years and oxygen classes for each lake with more than 15 years of data. Both rates and trends were computed with the Theil–Sen estimator74.
In water quality datasets, small lakes and ponds are typically under-represented relative to their natural abundance. Lakes with surface areas <10 ha comprise over three-quarters of global lakes75 but only 26% of our dataset. We therefore computed size-class abundance weighted medians of annual lake-mean oxygen concentrations (black line in Fig. 1). Abundance fractions were derived from the Swedish national lake registry, which lists the size class of over 100,000 water bodies (Svensk Vattenarkiv (SVAR) v.2016_3, retrieved 23 September 2022 from the Swedish Hydrological and Meteorological Institute). The census is based on a compilation of maps and is considered to be a complete record of Swedish lakes ≥1 ha (ref. 75).
At the level of individual lakes, we corrected for temporal sampling biases using GAM residuals. The number of samples per month peaks in August when seasonal oxygen minima occur most frequently, which could bias ice-free mean DO estimates low. The procedure is similar to detrending of a time series using a linear regression model. Here, we fitted an HGAM of the bottom-water DO concentration with a smooth for day-of-year (lake–year means), random intercepts for each lake and a Tweedie conditional distribution to account for the disproportionate number of zeros in the response variable. We fitted the GAM with the R package mgcv41. The day-of-year term was significant (P ≪ 0.0001). Model residuals were subsequently used in pooled time-series analyses (Fig. 1). We also detected a mean tendency towards sampling later in the year (1.6 ± 0.4 d decade−1, n = 951), which we adjusted for using GAM residuals before conducting the comparative analysis of DO and Δρ trends in Fig. 2. In the RF analyses of bottom-water oxygen saturation (trends), we included day-of-year (trends) as features.
We evaluated whether the variation in Δρ with lake morphometry (Fig. 2b) could mask a geographic climate forcing pattern. Overlapping time series of ERA5 reanalysis 2 m air temperature and 10 m wind speed across lake area and depth classes rule out a bias in climate forcing (Extended Data Fig. 4) and instead confirm that lake morphometry mediates lake thermal responses to climate warming.
We conducted RF analysis38 to identify the most important environmental covariates of ice-free mean bottom-water water DO saturation (1) and its long-term trend (2). We chose to model DOsat rather than concentration to distinguish temperature effects on solubility and biochemical conversion rates. RFs are based on an ensemble of regression trees, where each tree is trained on a bootstrap sample of data and a random subset of the independent variables splits each node. Model performance is evaluated as the error of its predictions for a subset of data not used to train the model (33% of the dataset). The importance of each variable corresponds to the decrease in model performance upon shuffling of its values (permutation). RF models were fitted with the Ranger package in R76. In RF (1) we included predictor variables that represent known drivers of DO5,20; stratification strength, surface-water TN, TP and colour, bottom-water temperature, Ased/Vwater and day-of-year. In RF (2) we added long-term trends of all variables, except Ased/Vwater. We used surface-water colour and nutrients as predictor variables because, in the hypolimnion of stratified lakes, anoxia-induced internal loading may confound causal relations with DOsat. Because RF cannot easily distinguish the importance of highly correlated variables (collinearity), we computed Kendall’s τ for each independent variable pair but found none that exceeded the recommended threshold of |τ| > 0.7 (ref. 77). We used the tuneMtryFast function to determine the optimal number of splitting variables and set up the RF with 1,500 trees to ensure stable outcomes. We created partial dependence plots with the R package ‘pdp’78 to visualize the marginal effect of each predictor variable. We re-fit model (1) replacing bottom trends with surface trends of TN, TP and colour to test hypotheses about the causal nature of the correlations (see main text).
To assess statistical significance of central tendency differences across classes (Fig. 2) we conducted one-way analysis of variance with the Kruskal–Wallis test79 and post hoc Conover–Iman tests for stochastic dominance with Benjamini–Hochberg adjustment for multiple comparisons80,81. The null hypothesis (no difference between groups) was rejected at Padj < 0.05. To test the hypothesis that long-term TP, TN and colour trends are greater in deep waters than in surface waters we performed pairwise one-sided Wilcoxon signed rank tests. Non-parametric tests were chosen because conditional distributions were not necessarily normal and because they tend to be more robust against outliers.
We fitted an HGAM to evaluate how the seasonal cycle of bottom-water DO varies between years and across lake morphometric gradients. GAMs are a class of generalized linear models, in which the response variable depends linearly on unknown smooth functions of the predictor variable(s)41,82. Hierarchical or mixed models contain fixed and random effects to explicitly model variance–covariance structures of non-independent, clustered data (for example, oxygen observations in lakes)83. The hierarchical model allowed us to incorporate information from the full dataset, rather than a subset of lakes with long-term, high-resolution time series. The HGAM was fitted with the mgcv package41 in R v.4.2 compiled with OpenBLAS v.0.3.15, on a high-performance computing cluster hosted by the National Academic Infrastructure for Supercomputing in Sweden (NAISS). As fixed effects we included smooth terms for day-of-year (DoY), year, relative sampling depth (zsample/zmax), log(surface area) and log(maximum depth) and their two- and three-way interactions. We modelled the effects as smooth functions that were represented in the model via cyclic (DoY) and cubic (all other variables) regression spline basis expansions of the covariates. We scaled day-of-year in each lake to long-term mean ice-on and ice-off dates from the ERA5 reanalysis dataset to account for geographical differences in ice-free season length. Including random slopes and intercepts for each lake resulted in a model that was too large to fit, so we opted for lake clusters within which we could reasonably expect highly synchronous oxygen cycles. To this end, the lakes were divided into 500 groups by standardized (zero mean and unit variance) latitude, longitude, annual mean air temperature, surface area and maximum depth via k-means clustering. The number of groups was limited by the memory capacity of the computing cluster. We modelled the conditional error distribution as a Tweedie location scale family84, in which the additional distribution parameters (power (p) and scale (φ)) are also modelled as sums of smooth functions of covariates, to account for distinct conditional error distributions, for example, in the mixing periods versus at the height of summer stratification.
Seasonal cycles were predicted for years 1970 and 2020 and a relative depth of 0.85 × zmax for four order-of-magnitude lake surface area classes with corresponding class-median maximum depths, using the ‘predict’ function in mgcv. Confidence intervals of predictions were computed from standard errors based on the Bayesian posterior covariance matrix of the parameters. We omitted more uncertain predictions for sparsely sampled data regions, such as small ponds (Asurf = 1 ha) and early sampling dates (1960s). Like most linear regression models, the HGAM provides expected means at specified values of the independent variables. Predictions therefore represent the expected mean DO concentrations of all lakes in a defined variable space (for example, 10 ha lakes in 1970) and should not be interpreted as ‘typical’ patterns in specific lake types or individual lakes. We estimated the long-term change in the duration of stratification by tracking the timing of the vernal and autumnal oxygen maxima. These maxima represent model expectations of the time when most lakes of the specified morphology reach peak DO concentrations. Mean peak locations and 95% credible intervals (CIs) were computed from predictions based on repeated draws from the posterior distribution of model parameters (n = 1,000) using the fitted_samples function in the gratia R package85.
Water quality monitoring data used in this study are available via Zenodo at https://doi.org/10.5281/zenodo.11243331 (ref. 36). We did not receive permission to republish data for Lake identifier 8051. Original data sources, including open access data repositories and online data request forms, are listed in Supplementary Data 1. Morphometric properties of the study lakes were collected from open access data repositories, literature sources and water quality datasets listed in Supplementary Data 2. ECMWF ERA5 Reanalysis climate data are available through the Copernicus Climate Data Store: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview.
R code to perform the analyses and generate the figures is available via Zenodo at https://doi.org/10.5281/zenodo.11243331 (ref. 36).
Jenny, J.-P. et al. Urban point sources of nutrients were the leading cause for the historical spread of hypoxia across European lakes. Proc. Natl Acad. Sci. USA 113, 12655–12660 (2016).
Article CAS Google Scholar
Yuan, L. L. & Jones, J. R. Modeling hypolimnetic dissolved oxygen depletion using monitoring data. Can. J. Fish. Aquat. Sci. 77, 814–823 (2020).
Article Google Scholar
Jane, S. F. et al. Widespread deoxygenation of temperate lakes. Nature 594, 66–70 (2021).
Article CAS Google Scholar
Wetzel, R. G. Limnology (Elsevier, 2001).
Livingstone, D. M. & Imboden, D. M. The prediction of hypolimnetic oxygen profiles: a plea for a deductive approach. Can. J. Fish. Aquat. Sci. 53, 924–932 (1996).
Article Google Scholar
Ficke, A. D., Myrick, C. A. & Hansen, L. J. Potential impacts of global climate change on freshwater fisheries. Rev. Fish. Biol. Fish. 17, 581–613 (2007).
Article Google Scholar
Till, A., Rypel, A. L., Bray, A. & Fey, S. B. Fish die-offs are concurrent with thermal extremes in north temperate lakes. Nat. Clim. Change 9, 637–641 (2019).
Article Google Scholar
Carey, C. C. et al. Anoxia decreases the magnitude of the carbon, nitrogen and phosphorus sink in freshwaters. Glob. Change Biol. 28, 4861–4881 (2022).
Article Google Scholar
Hupfer, M. & Lewandowski, J. Oxygen controls the phosphorus release from lake sediments—a long-lasting paradigm in limnology. Int. Rev. Hydrobiol. 93, 415–432 (2008).
Article CAS Google Scholar
Michalak, A. M. et al. Record-setting algal bloom in Lake Erie caused by agricultural and meteorological trends consistent with expected future conditions. Proc. Natl Acad. Sci. USA 110, 6448–6452 (2013).
Article CAS Google Scholar
Pajala, G. et al. The effects of water column dissolved oxygen concentrations on lake methane emissions—results from a whole‐lake oxygenation experiment. J. Geophys. Res. Biogeosci. 128, e2022JG007185 (2023).
Brooks, B. W. et al. Are harmful algal blooms becoming the greatest inland water quality threat to public health and aquatic ecosystems? Environ. Toxicol. Chem. 35, 6–13 (2016).
Article CAS Google Scholar
Cornett, R. J. & Rigler, F. H. The areal hypolimnetic oxygen deficit: an empirical test of the model. Limnol. Oceanogr. 25, 672–679 (1980).
Article CAS Google Scholar
Monteith, D. T. et al. Dissolved organic carbon trends resulting from changes in atmospheric deposition chemistry. Nature 450, 537–540 (2007).
Article CAS Google Scholar
De Wit, H. A. et al. Cleaner air reveals growing influence of climate on dissolved organic carbon trends in northern headwaters. Environ. Res. Lett. 16, 1–13 (2021).
CAS Google Scholar
Škerlep, M., Steiner, E., Axelsson, A. & Kritzberg, E. S. Afforestation driving long‐term surface water browning. Glob. Change Biol. 26, 1390–1399 (2020).
Article Google Scholar
Foley, B., Jones, I. D., Maberly, S. C. & Rippey, B. Long-term changes in oxygen depletion in a small temperate lake: effects of climate change and eutrophication. Freshw. Biol. 57, 278–289 (2012).
Article CAS Google Scholar
Jane, S. F. et al. Longer duration of seasonal stratification contributes to widespread increases in lake hypoxia and anoxia. Glob. Change Biol. 29, 1009–1023 (2023).
Article CAS Google Scholar
Gudasz, C. et al. Temperature-controlled organic carbon mineralization in lake sediments. Nature 466, 478–481 (2010).
Article CAS Google Scholar
Nürnberg, G. K. Quantification of anoxia and hypoxia in water bodies. In Encyclopedia of Water: Science, Technology, and Society (ed. Maurice, P.) 525–534 (Wiley, 2019).
Cornett, R. J. & Rigler, F. H. Decomposition of seston in the hyolimnion. Can. J. Fish. Aquat. Sci. 44, 146–151 (1987).
Article Google Scholar
Messager, M. L., Lehner, B., Grill, G., Nedeva, I. & Schmitt, O. Estimating the volume and age of water stored in global lakes using a geo-statistical approach. Nat. Commun. 7, 13603 (2016).
Article CAS Google Scholar
Rantanen, M. et al. The Arctic has warmed nearly four times faster than the globe since 1979. Commun. Earth Environ. 3, 168 (2022).
Article Google Scholar
Woolway, R. I. et al. Phenological shifts in lake stratification under climate change. Nat. Commun. 12, 2318 (2021).
Woolway, R. I. et al. Northern hemisphere atmospheric stilling accelerates lake thermal responses to a warming world. Geophys. Res. Lett. 46, 11983–11992 (2019).
Article Google Scholar
Li, X., Peng, S., Xi, Y., Woolway, R. I. & Liu, G. Earlier ice loss accelerates lake warming in the Northern Hemisphere. Nat. Commun. 13, 5156 (2022).
Article CAS Google Scholar
Fölster, J., Johnson, R. K., Futter, M. N. & Wilander, A. The Swedish monitoring of surface waters: 50 years of adaptive monitoring. Ambio 43, 3–18 (2014).
Article Google Scholar
Couture, R., Wit, H. A., Tominaga, K., Kiuru, P. & Markelov, I. Oxygen dynamics in a boreal lake responds to long‐term changes in climate, ice phenology and DOC inputs. J. Geophys. Res. Biogeosci. 120, 2441–2456 (2015).
Article CAS Google Scholar
Crossman, J. et al. Can recovery from disturbance explain observed declines in total phosphorus in Precambrian Shield catchments? Can. J. Fish. Aquat. Sci. 73, 1202–1212 (2016).
Article CAS Google Scholar
Rekolainen, S., Mitikka, S., Vuorenmaa, J. & Johansson, M. Rapid decline of dissolved nitrogen in Finnish lakes. J. Hydrol. 304, 94–102 (2005).
Article CAS Google Scholar
Huser, B. J., Futter, M. N., Wang, R. & Fölster, J. Persistent and widespread long-term phosphorus declines in Boreal lakes in Sweden. Sci. Total Environ. 613–614, 240–249 (2018).
Article Google Scholar
Isles, P. D. F. et al. Widespread synchrony in phosphorus concentrations in northern lakes linked to winter temperature and summer precipitation. Limnol. Oceanogr. Lett. 8, 639–648 (2023).
Article CAS Google Scholar
Sharma, S. et al. Widespread loss of lake ice around the Northern Hemisphere in a warming world. Nat. Clim. Change 9, 227–231 (2019).
Article Google Scholar
Mammarella, I. et al. Effects of similar weather patterns on the thermal stratificationmixing regimes and hypolimnetic oxygen depletion in two boreal lakes with different water transparency. Boreal Environ. Res. 23, 237–247 (2018).
Google Scholar
Dugan, H. A. A comparison of ecological memory of lake ice‐off in eight north‐temperate lakes. J. Geophys. Res. Biogeosci. 126, e2020JG006232 (2021).
Jansen, J. Lake water quality monitoring data from 8287 northern lakes [Data set]. Zenodo https://doi.org/10.5281/zenodo.11243331 (2024).
Dinerstein, E. et al. An ecoregion-based approach to protecting half the terrestrial realm. Bioscience 67, 534–545 (2017).
Article Google Scholar
Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001).
Article Google Scholar
Soranno, P. A. et al. Cross‐scale interactions: quantifying multi‐scaled cause–effect relationships in macrosystems. Front. Ecol. Environ. 12, 65–73 (2014).
Article Google Scholar
Nakhaei, N., Ackerman, J. D., Bouffard, D., Rao, Y. R. & Boegman, L. Empirical modeling of hypolimnion and sediment oxygen demand in temperate Canadian lakes. Inl. Waters 11, 351–367 (2021).
Article CAS Google Scholar
Wood, S. N. Generalized Additive Models (Chapman and Hall/CRC, 2017).
MacIntyre, S. & Melack, J. M. in Encyclopedia of Inland Waters (ed. Likens, G. E.) 603–612 (Elsevier, 2009). https://doi.org/10.1016/B978-012370626-3.00040-5
Winslow, L. A., Read, J. S., Hansen, G. J. A. & Hanson, P. C. Small lakes show muted climate change signal in deepwater temperatures. Geophys. Res. Lett. 42, 355–361 (2015).
Article Google Scholar
Magee, M. R. & Wu, C. H. Response of water temperatures and stratification to changing climate in three lakes with different morphometry. Hydrol. Earth Syst. Sci. 21, 6253–6274 (2017).
Article Google Scholar
Jankowski, T., Livingstone, D. M., Bührer, H., Forster, R. & Niederhauser, P. Consequences of the 2003 European heat wave for lake temperature profiles, thermal stability and hypolimnetic oxygen depletion: implications for a warmer world. Limnol. Oceanogr. 51, 815–819 (2006).
Article Google Scholar
Nelligan, C., Jeziorski, A., Rühland, K. M., Paterson, A. M. & Smol, J. P. Long-term trends in hypolimnetic volumes and dissolved oxygen concentrations in Boreal Shield lakes of south-central Ontario, Canada. Can. J. Fish. Aquat. Sci. 76, 2315–2325 (2019).
Article CAS Google Scholar
Bengtsson, L. & Ali-Maher, O. The dependence of the consumption of dissolved oxygen on lake morphology in ice covered lakes. Hydrol. Res. 51, 381–391 (2020).
Article CAS Google Scholar
Vachon, D. & Prairie, Y. T. The ecosystem size and shape dependence of gas transfer velocity versus wind speed relationships in lakes. Can. J. Fish. Aquat. Sci. 70, 1757–1764 (2013).
Article Google Scholar
Vollenweider, R. A. Scientific Fundamentals of the Eutrophication of Lakes and Flowing Waters, with Particular Reference to Nitrogen and Phosphorus as Factors in Eutrophication (Organization for Economic Co-operation and Development, 1968).
Brothers, S. et al. A feedback loop links brownification and anoxia in a temperate, shallow lake. Limnol. Oceanogr. 59, 1388–1398 (2014).
Article CAS Google Scholar
Lapierre, J.-F., Guillemette, F., Berggren, M. & del Giorgio, P. A. Increases in terrestrially derived carbon stimulate organic carbon processing and CO2 emissions in boreal aquatic ecosystems. Nat. Commun. 4, 2972 (2013).
Article Google Scholar
Skoog, A. C. & Arias-Esquivel, V. A. The effect of induced anoxia and reoxygenation on benthic fluxes of organic carbon, phosphate, iron and manganese. Sci. Total Environ. 407, 6085–6092 (2009).
Article CAS Google Scholar
Peter, S., Agstam, O. & Sobek, S. Widespread release of dissolved organic carbon from anoxic boreal lake sediments. Inl. Waters 7, 151–163 (2017).
Article CAS Google Scholar
Lau, M. P., Hutchins, R. H. S., Tank, S. E. & A. del Giorgio, P. The chemical succession in anoxic lake waters as source of molecular diversity of organic matter. Sci. Rep. 14, 3831 (2024).
Article CAS Google Scholar
Müller, B., Bryant, L. D., Matzinger, A. & Wüest, A. Hypolimnetic oxygen depletion in eutrophic lakes. Environ. Sci. Technol. 46, 120824131307006 (2012).
Article Google Scholar
Ghane, A. & Boegman, L. The dissolved oxygen budget of a small Canadian Shield lake during winter. Limnol. Oceanogr. 68, 265–283 (2023).
Article Google Scholar
Kleeberg, A., Herzog, C. & Hupfer, M. Redox sensitivity of iron in phosphorus binding does not impede lake restoration. Water Res. 47, 1491–1502 (2013).
Article CAS Google Scholar
Tammeorg, O., Nürnberg, G. K., Nõges, P. & Niemistö, J. The role of humic substances in sediment phosphorus release in northern lakes. Sci. Total Environ. 833, 155257 (2022).
Article CAS Google Scholar
Peter, S. & Sobek, S. High variability in iron-bound organic carbon among five boreal lake sediments. Biogeochemistry 139, 19–29 (2018).
Article CAS Google Scholar
Jane, S. F. et al. Concurrent warming and browning eliminate cold-water fish habitat in many temperate lakes. Proc. Natl Acad. Sci. USA 121, 2017 (2024).
Article Google Scholar
Muñoz Sabater, J. ERA5-Land Hourly data from 1950 to Present (Copernicus Climate Change Service, 2019); https://cds.climate.copernicus.eu/cdsapp#!/dataset/10.24381/cds.e2161bac
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2022).
Read, E. K. et al. Water quality data for national‐scale aquatic research: the Water Quality Portal. Water Resour. Res. 53, 1735–1745 (2017).
Article Google Scholar
Chen, C.-T. & Millero, F. J. The use and misuse of pure water PVT properties for lake waters. Nature 266, 707–708 (1977).
Article CAS Google Scholar
Benson, B. B. & Krause, D. The concentration and isotopic fractionation of oxygen dissolved in freshwater and seawater in equilibrium with the atmosphere. Limnol. Oceanogr. 29, 620–632 (1984).
Article CAS Google Scholar
Gray, E., Mackay, E. B., Elliott, J. A., Folkard, A. M. & Jones, I. D. Wide-spread inconsistency in estimation of lake mixed depth impacts interpretation of limnological processes. Water Res. 168, 115136 (2020).
Article CAS Google Scholar
Last, W. M. & Ginn, F. M. Saline systems of the Great Plains of western Canada: an overview of the limnogeology and paleolimnology. Saline Syst. 1, 10 (2005).
Article Google Scholar
Carlson, R. E. A trophic state index for lakes. Limnol. Oceanogr. 22, 361–369 (1977).
Article CAS Google Scholar
Lyche Solheim, A. et al. Water Framework Directive Intercalibration Technical Report—Northern Lake Phytoplankton Ecological Assessment Methods (ed. Poikane, S.) (EU Publications Office, 2014); https://doi.org/10.2788/70684
Winslow, L. et al. rLakeAnalyzer: Lake Physics Tools v.1.11.4.1 (2019).
Shaughnessy, A. R., Wen, T., Niu, X. & Brantley, S. L. Three principles to use in streamlining water quality research through data uniformity. Environ. Sci. Technol. 53, 13549–13550 (2019).
Article CAS Google Scholar
Mironov, D. V. Parameterization of Lakes in Numerical Weather Prediction. Description of a Lake Model (DWD, 2008).
Bernhardt, J., Engelhardt, C., Kirillin, G. & Matschullat, J. Lake ice phenology in Berlin-Brandenburg from 1947–2007: observations and model hindcasts. Clim. Change 112, 791–817 (2012).
Article Google Scholar
Sen, P. K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 63, 1379–1389 (1968).
Article Google Scholar
Cael, B. B. & Seekell, D. A. The size-distribution of Earth’s lakes. Sci. Rep. 6, 29633 (2016).
Article CAS Google Scholar
Wright, M. N. & Ziegler, A. ranger: a fast implementation of random forests for high dimensional data in C++ and R. J. Stat. Softw. 77, 1752–1758 (2017).
Article Google Scholar
Dormann, C. F. et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).
Article Google Scholar
Greenwell, B. M. pdp: an R package for constructing partial dependence plots. R J. 9, 421–436 (2017).
Article Google Scholar
Kruskal, W. H. & Wallis, W. A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621 (1952).
Article Google Scholar
Conover, W. J. & Iman, R. L. Rank transformations as a bridge between parametric and nonparametric statistics. Am. Stat. 35, 124–129 (1981).
Article Google Scholar
Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300 (1995).
Article Google Scholar
Pedersen, E. J., Miller, D. L., Simpson, G. L. & Ross, N. Hierarchical generalized additive models in ecology: an introduction with mgcv. PeerJ 7, e6876 (2019).
Article Google Scholar
Zuur, A. F., Ieno, E. N., Walker, N., Saveliev, A. A. & Smith, G. M. Mixed Effects Models and Extensions in Ecology with R (Springer, 2009).
Wood, S. N. & Fasiolo, M. A generalized Fellner–Schall method for smoothing parameter optimization with application to Tweedie location, scale and shape models. Biometrics 73, 1071–1081 (2017).
Article Google Scholar
Simpson, G. L. gratia: Graceful ggplot-based graphics and other functions for GAMs fitted using mgcv. R package version 0.8.1 (2023).
Massicotte, P. & South, A. rnaturalearth: World Map Data from Natural Earth. R package version 1.0.1.9000 https://github.com/ropensci/rnaturalearth (2024).
Muñoz Sabater, J. ERA5-Land monthly averaged data from 1950 to present (Copernicus Climate Change Service, 2019); https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means
Download references
This work was funded by grant 2020-06460 from the Swedish Research Council to J.J. We wish to thank all staff and volunteers who have, over many decades, contributed to lake water quality monitoring. We are grateful to the Finnish Environment Institute, the Finnish Centres for Economic Development, Transport and the Environment, the Swedish University of Agricultural Sciences, the Swedish Meteorological and Hydrological Institute, the Norwegian Environment Agency, the Norwegian Water Resources and Energy Directorate, Environment and Climate Change Canada, Fisheries and Oceans Canada, the Ontario Ministry of the Environment, Conservation and Parks, Dorset Environmental Science Centre, IISD Experimental Lakes Area, Turkey Lake Watershed Programme, the British Columbia Ministry of Environment and Climate Change Strategy, the Nova Scotia Department of Environment and Climate Change, the New Brunswick Departments of Environment and Local Government and Natural Resources and Energy Development, the New Brunswick Lakes Monitoring Program, the Government of Alberta, the Alberta Geological Survey, the Alberta Lake Management Society, Lac La Biche County, Prince Edward Island Department of Environment, Water and Climate Change, Parks Canada (Riding Mountain National Park, Wood Buffalo National Park), Province of Manitoba Environment and Climate Change, Manitoba Hydro, the Saskatchewan Water Security Agency, the Saskatchewan Ministry of Environment, Ministère de l’Environnement, de la Lutte contre les changements climatiques, de la Faune et des Parcs (Québec), the Northwest Territories Ministry of Environment and Climate Change, the US National Water Quality Monitoring Council, the Alaska Departments of Environmental Conservation and Fish and Game, the US Long-term Ecological Research Network and the Toolik Lake GIS Office at the University of Alaska Fairbanks. We thank J. Fölster, S. Mitikka, M. Paterson, S. Moras, C. Rickard, M. Nicholson, J. Kennedy, S. Hase, M. Rahman, L. Roy, S. Bourget, R. Staples, A. Bethe, A. Giblin and G. Kling for help with data processing and quality control. The National Academic Infrastructure for Supercomputing in Sweden (NAISS) provided computing resources. G.A.W. acknowledges financial support for this study from the Swedish Research Council (VR grant no. 2020-03222 and FORMAS grant no. 2020-01091). L.H.H. was supported by the Finnish Ministry of Agriculture and Forestry under the programme Catch the Carbon (project SysteemiHiili (VN/28536/2020)). P.A.d.G. was supported by the CarBBAS (Carbon Biogeochemistry in Boreal Aquatic Systems) programme funded by NSERC and Hydro-Québec.
Open access funding provided by Uppsala University.
Department of Ecology and Genetics/Limnology, Uppsala University, Uppsala, Sweden
Joachim Jansen & Gesa A. Weyhenmeyer
Université du Québec à Montréal, Montreal, Canada
Joachim Jansen, Paul A. del Giorgio & Yves T. Prairie
Groupe de Recherche Interuniversitaire en Limnologie, Montreal, Canada
Joachim Jansen, Paul A. del Giorgio & Yves T. Prairie
Department of Animal and Veterinary Sciences and Climate, Aarhus University, Aarhus, Denmark
Gavin L. Simpson
Finnish Environment Institute, Helsinki, Finland
Laura H. Härkönen
Ontario Ministry of the Environment, Conservation and Parks, Dorset Environmental Science Centre, Dorset, Ontario, Canada
Andrew M. Paterson
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
You can also search for this author in PubMed Google Scholar
J.J. conceptualized the study, acquired funding, compiled the dataset, carried out the analyses and wrote the manuscript. L.H.H. and A.M.P. contributed data to the study. G.A.W., Y.T.P. and P.A.d.G. provided additional funding and periodic feedback on the analyses. G.L.S. guided the statistical analyses. All authors contributed critical feedback to draft revisions of the manuscript.
Correspondence to Joachim Jansen.
The authors declare no competing interests.
Nature Climate Change thanks Mohammed Hamdan, Yingxun Du, Claudia Dresti and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Panels show Fennoscandia (a) and northern North America (b). Lake locations are shown in blue. Orange dots denote ‘trend’ lakes with more than 15 years of bottom-water DO observations. Coloured regions mark terrestrial biomes37. The biome ‘Temperate Forest and Grasslands’ includes ‘Temperate Conifer Forests’, ‘Temperate Broadleaf and Mixed Forests’ and ‘Temperate Grasslands, Savannas and Shrublands’. World map data is from the ‘rnaturalearth’ R package86. The ESRI shapefile of the terrestrial biomes37 can be accessed at https://storage.googleapis.com/teow2016/Ecoregions2017.zip.
a,b. Box plots of long-term dissolved oxygen trends (≥15 years) by lake surface area class in distinct regions (a) and terrestrial biomes (b)37. Boxes bound the interquartile range (25th-75th percentile, IQR) and whiskers extend to the smaller quantity of data extremes or medians ± 1.5 × IQR. Dots and lines mark means and medians, respectively. Outliers are shown by grey dots. Bottom labels show the sample size (number of trend lakes) in each surface area class. The biome ‘Temperate Forest and Grasslands’ includes ‘Temperate Conifer Forests’, ‘Temperate Broadleaf and Mixed Forests’ and ‘Temperate Grasslands, Savannas and Shrublands’. The ESRI shapefile of the terrestrial biomes37 can be accessed at https://storage.googleapis.com/teow2016/Ecoregions2017.zip.
Ranked feature importance of environmental covariates of DO saturation (a) and the long-term DO saturation trend (b). Permutation importance z-scores reflect increases in model root-mean-square-error upon shuffling variable values. Bracketed letters indicate whether the variable was measured in bottom (b) or surface waters (s).
ERA5-Land Reanalysis Monthly87 2 m air temperature (a,b) and 10 m wind speed (c,d) across classes of surface area (a,c) and maximum depth (b,d). Legend labels define area class boundaries (for example green marks Asurf = 102–103 ha). Lines connect class-medians of ice-free (June-October) lake-year means, computed for each lake over 1960–2022 (nlakes = 8288). Black lines show linear regressions with 95% confidence interval bands. The figure shows that climate forcing is not biased toward specific morphometric lake classes.
Provenance of the lake water quality monitoring dataset.
Provenance of the lake morphometry dataset.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and permissions
Jansen, J., Simpson, G.L., Weyhenmeyer, G.A. et al. Climate-driven deoxygenation of northern lakes. Nat. Clim. Chang. (2024). https://doi.org/10.1038/s41558-024-02058-3
Download citation
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41558-024-02058-3
Anyone you share the following link with will be able to read this content:
Sorry, a shareable link is not currently available for this article.
Provided by the Springer Nature SharedIt content-sharing initiative
Advertisement
© 2024 Springer Nature Limited
Sign up for the Nature Briefing: Microbiology newsletter — what matters in microbiology research, free to your inbox weekly.