Anaerobic duration predicts biogeochemical consequences of oxygen depletion in lakes

Lake deoxygenation is of growing concern because it threatens ecosystem services delivery. Complete deoxygenation, anoxia, is projected to prolong and expand in lakes, promoting the production or release of nutrients, greenhouse gases and metals from the water column and sediments. Accumulation of these compounds cannot be easily predicted thus hindering our capacity to forecast the ecological consequences of global changes on aquatic ecosystems. Here, we used monitoring data of four lakes to develop a novel tool, anaerobic duration, to study anaerobic processes in lake waters. Anaerobic duration explained, as a single predictor, 21–60% of the variation for ammonium, phosphorus and a dissolved organic matter fluorophore. Anaerobic duration could be modeled using only two oxygen profiles and lake bathymetry, making it an easily applicable tool to interpret and extrapolate biogeochemical data. This novel tool thus has the potential to transform widely available oxygen profiles into an ecologically meaningful variable.

Lakes provide essential ecosystem services (Jane et al. 2021), several of which are threatened by anthropogenic activities. Both eutrophication and global warming critically affect dissolved oxygen (DO) in lakes through higher hypolimnetic oxygen demand (Müller et al. 2012) or reduced hypolimnetic ventilation (Bartosiewicz et al. 2019). These pressures thus threaten more lake hypolimnia to become or stay anoxic for longer episodes (Matzinger et al. 2010;Jenny et al. 2016). DO depletion in hypolimnia have far-reaching ecological consequences, including the accumulation of reduced compounds toxic to organisms, loss of habitats and intensified production of greenhouse gases (Jane et al. 2021, and references therein). Anoxia is also associated with phosphorus release from redox-sensitive sediment components (Hupfer and Lewandowski 2008). Monitoring and forecasting of these threats can be improved through modeling of oxygen dynamics.
Models aim to condense the spatiotemporal oxygen dynamics into indices. Bulk indices reflect the areal or volumetric hypolimnetic oxygen demand (AHOD and VHOD) (Strom 1931;Rosa and Burns 1987). Higher-resolution models use DO depletion from several hypolimnetic water layers to differentiate benthic and pelagic demand (J z , Livingstone and Imboden 1996). These DO consumption models are usually tuned to single lakes and rely on oxygen data and lake bathymetry and are convenient tools to study drivers, onset and extent of anoxia, but not its consequences. Approaches to estimate consequences of anoxia such as phosphorus accumulation consider the lake-scale proportion (anoxic factor [AF]) of either sediment surface (Nürnberg 1984) or water volume (Foley et al. 2012) affected by anoxia, but are insufficient to deconstruct the temporal sequence of anoxiarelated processes (Matzinger et al. 2010). Currently, no modeling approach includes all elements required to fully address lake anoxia, including time, benthic and pelagic prokaryotic activity and their associated metabolite dynamics in the water column.
Lake water column chemistry reflects the material taken-up and released by photosynthetic and heterotrophic organisms. Anoxia is usually restricted to the aphotic zone, dominated by heterotrophs that consume, process and release dissolved organic matter (DOM) including fluorescent (FDOM; Burdige et al. 2004) compounds. Under anoxic conditions, organic matter is not only the precursor of anaerobic metabolic products like ammonium (NH 4 + ) and methane, but also enables the accumulation of oxygen-sensitive toxins such as methylmercury and hydrogen sulfide (Ach a et al. 2018). Thus, the extent and duration of anaerobic conditions determine the quantity of these substances that will be introduced to surface waters in subsequent turnover events. However, easily predicting their accumulation remains a challenge.
In this study, we explored how to predict the accumulation of common anaerobic metabolite solely from lake oxygen data. To this end, we used anaerobic duration (AD), a novel continuous variable containing spatiotemporal information which indicates how long discrete hypolimnion layers were DO depleted. This aggregation is possible because vertical exchange in stratified hypolimnetic waters is considered negligible (Rippey and McSorley 2009), and major solutes originate from the water column or conterminous sediments (Livingstone and Imboden 1996). AD calculation requires daily oxygen profiles, either observed or modeled, and in contrast to AF, AD attributes a value for each depth and day under anoxia and scales with DO dynamics.
To establish AD as a widely applicable tool to study anoxia, we used the Livingstone and Imboden (1996) model between J z and sediment-area-to-water-volume ratio (α(z)) and tested linear and non-linear equations to best model daily oxygen profiles. We calculated AD from DO measurements obtained from either a limited number of profiles or from several loggers, two widely available lake monitoring data formats, and compared their modeled AD (AD mod ) to observed AD (AD obs ) in a specifically instrumented lake, running high-resolution DO profiler in (bi-)daily casts. Finally, we used AD to predict anaerobic metabolites, including NH 4 + , soluble reactive phosphorus (SRP), total phosphorus (TP), and a FDOM component in four lakes. AD proves to be an easily modeled variable describing anaerobic lake biogeochemistry that can accurately predict a wide array of compound accumulation. We argue that not only the compounds used herein, but all compounds produced as by-products of anaerobic metabolism may be modeled with high precision and low effort.

AD calculation
The AD concept transforms oxygen data below a specified threshold (DO threshold ) into a spatiotemporal variable that reflects the time passed since a parcel of water crossed the threshold. Here, we used a DO threshold of 0.5 mg L À1 based on the best R 2 and Akaike information criterion (AIC) of a sensitivity analysis of various thresholds (Table S2), but others can LaBrie et al.
Novel tool describing anoxia in lakes be used depending on research question. Following DO decrease in hypolimnetic waters, AD increases for every consecutive timestep (Day) between the first (i) and last (n) day of stratification a water stratum DO concentration (DO z ) is below DO threshold (Eq. 1). If DO z increases above DO threshold , AD returns to 0. This results in a matrix of values (depths, dates) graphically represented in Fig. S1. Similarly, oxic duration (OD; days) of hypolimnetic water can be calculated for each depth-date pair as the time since stratification onset while DO z > DO threshold .
AD calculation requires DO concentration in high spatiotemporal resolution (daily and every meter). Since measurements in this resolution are rarely available, daily oxygen values may be interpolated from J z , themselves calculated from a few oxygen profiles or loggers (Livingstone and Imboden 1996). To assess if AD may be accurately predicted from all common measurement practices, we subsampled observations from the full Arendsee dataset (see below).

Modeling oxygen
We calculated J z using the slope of linear regressions between DO concentration and time during stagnation ("casts"), either with manually selected dates or a two breakpoints segmented regression (Muggeo 2008); both yielded similar values (Fig. S2). We also calculated J z for all oxic profile casts taken 28 days apart as a reduced sampling "two oxic profiles" scenario, and for all three hypolimnetic DO loggers combinations ("loggers"). J z are reported in Tables S3-S6. We then compared three different types of equations to best describe the relationship between J z and α(z) (see Supporting Information for α(z) calculation). We used a linear fit, and to capture the asymptotic behavior, both log-linear and exponential-plateau fits (Eq. 2), where b (non-negative) and k are fitting parameters, and J z,max is the maximum J z fitted as a random parameter by the nlsLM function (Elzhov et al. 2016). The best model was chosen using R 2 , lowest rootmean-square error (RMSE) and AIC.
Using these equations, we modeled daily DO profile series by assuming that the entire hypolimnion was fully oxygenated at onset of summer stratification. Because onset date is unknown, we needed to calibrate dates. To do so, we assigned dates to the DO profile series by matching a measured profile using lowest RMSE. For Arendsee, we had daily measurements to calculate AD obs , but needed to use AD mod for the other lakes. We calculated ADs (Eq. 1) and compared the modeled and observed first date (as day of year, DOY) on which each stratum became hypoxic (DO < 2mg L À1 ) to assess quality of AD mod .

Chemical data
Concentrations of SRP and NH 4 + were determined photometrically using molybdenum blue (Murphy and Riley 1962) and indophenol (Bolleter et al. 1961), respectively, by segmented flow analysis (Scan++; Skalar Analytical) for Arendsee and Stechlin. SRP and NH 4 + values taken at non-integer depths were rounded (floored at 0.5 m). For Mendota, Exo2 FDOM was measured at excitation wavelength 365 nm, emission wavelength 480 nm and expressed in quinine sulfate units. This peak (F 365/480 , F ex/em ) usually indicates terrestrially derived recalcitrant compounds (Coble 1996). Delavan TP data were taken from Jane et al. (2021). All reported statistical parameters were statistically significant at p value < 0.01 unless specified. All statistics and modeling were performed using R version 4.2.2 (R Core Team 2022), and all data are available at https://fred.igb-berlin.de/data/ package/819 (LaBrie et al. 2023).

Modeling AD
As AD is calculated from daily DO profiles, we first assessed which equation best fitted J z to α(z) to model profile time series. This is critical to interpolate missing depths from irregular profile measurements or depth-discrete loggers. The best fit was generally the exponential-plateau, followed by the loglinear equation (Figs. 1, S3; Table S7). When excluding deeper layers by subsampling J z to simulate partial anaerobic conditions, we found that the exponential-plateau model was slightly less resilient than the log-linear relationship (Fig. S4). Overall, the linear fit was inferior to the other equations (Figs. 1, S3) and loggers had lower fits regardless of equations. Subsampling three loggers or using only two oxic profiles provided reliable J z -α(z) relationships in most cases (Fig. S5).
Using linear, log-linear and exponential-plateau fits, we then compared estimates of the date of hypoxia first appearance. For all years and most data source-model type pairs, the log-linear and exponential-plateau models yielded good results, whereas the linear relationship did not (Figs. 1b, S6-S9).

Nutrients and FDOM
We analyzed patterns in SRP, TP, and NH 4 + from hypolimnetic waters in relation to OD, AD obs (Arendsee) and AD mod (Stechlin and Delavan) (Fig. 2). The influence of OD on nutrients varied between phosphorus and nitrogen and among lakes. In Stechlin and Arendsee, SRP slowly accumulated during the oxic stagnant period at 0.29 AE 0.03 and 0.22 AE 0.03 μg SRP L À1 d À1 (R 2 = 0.25 and R 2 = 0.40, Fig. 2a,e), respectively, but not in Delavan (R 2 = 0.02, p = 0.53; Fig. 2i). In all lakes (Fig. 2b,f,k), the accumulation rate with AD was higher at 0.65 AE 0.02 μg SRP L À1 d À1 for Stechlin (R 2 = 0.21), 0.88 AE 0.01 μg SRP L À1 d À1 for Arendsee (R 2 = 0.53) and at 5.6 AE 0.8 μg TP L À1 d À1 for Delavan (R 2 = 0.60). Ammonium was highly dynamic in Stechlin (Fig. 2c,d). In the deepest layers in 2020 and 2021, there was an increase during the oxic period in Stechlin, which was rapidly removed after 3 weeks under hypoxia (< 2 mgO 2 L À1 ), whereas in other depths and years it remained low. Under anaerobic conditions, there was an increase at 3.4 AE 0.4 μg NH 4 + L À1 d À1 , and a subsequent decrease in late 2019. In Arendsee (Fig. 2g,h), there was an exponential decrease during the oxic period (k = 0.018 AE 0.006 d À1 , R 2 = 0.22) followed by strong increase once anaerobic conditions settled (7.2 AE 1.3 μg NH 4 + L À1 d À1 , R 2 = 0.50). We extracted phosphorus accumulation rates from both OD and AD relationships for each year and lake (N = 9). We observed a decreasing trend in phosphorus-OD slopes with trophic status and an increasing one for phosphorus-AD (Fig. 3). We explored patterns of the fluorescent component F 365/480 in Mendota using OD and AD mod . This fluorophore increased during the oxic period at 0.0316 AE 0.0008 (R 2 = 0.94) and at 0.0275 AE 0.0019 (R 2 = 0.74) QSE d À1 in 2018 and 2020, respectively (Fig. 4a,c). The relationships with AD mod were weaker with R 2 of 0.54 and 0.66 for 2018 and 2020, respectively (Fig. 4b,d). Anaerobic conditions seemed to have a reproducible, stabilizing influence on this fluorophore for about 30 days; a behavior not well reflected by linear models and only visualized using AD.

Discussion
Heterotrophic respiration and oxidation of reduced compounds drive oxygen loss in lake hypolimnia (Steinsberger et al. 2020). To study and predict the specific lake water biogeochemistry that begins with deoxygenation, we developed a novel, ecologically meaningful variable, AD. By regressing limnological variables (e.g., SRP, NH 4 + ) over AD, we obtain the overall metabolism of the anaerobic hypolimnion and a daily accumulation rate for each compound. In turn, this enables prediction of said compound total loads in a warmer climate where stratification and anoxia prolong (Jane et al. 2023), provided that destratification dynamics can be adequately modeled. To this end, we built on the deductive approach  + , c, d, g, h) on the right column, as a function of oxic duration (a, c, e, g, i) and anaerobic duration (b, d, f, h, k) in lake Stechlin, Arendsee, and Delavan, respectively. Colored symbols represent different years for the same lake but note the difference among lakes. SRP: soluble reactive phosphorus; TP: total phosphorus. Points in parenthesis in panels (b) and (d) are thought to represent a slow reaeration process oxidizing these compounds before oxygen accumulation; in other panels, they are outliers. None of these points are included in linear regressions. Dashed line in panel (i) represent a statistically non-significant regression line.
that deconstructs oxygen depletion rates using lake bathymetry (Livingstone and Imboden 1996).
Our approaches are reminiscent of other oxygen metrics AHOD and VHOD (Strom 1931;Rosa and Burns 1987), and the AF (Nürnberg 1984), yet they are not interchangeable. For instance, AHOD and VHOD bulk hypolimnetic oxygen, normalize according to bathymetry, to then study drivers of oxygen consumption including nutrient load (e.g., Cornett and Rigler 1980;Rippey and McSorley 2009). AF extends this approach beyond the oxic/anoxic transition and extends to the lake-scale, where it is used to predict internal phosphorus loading (Nürnberg 1984) or fish richness (Nürnberg 1995). In contrast, AD is designed to study within-lake dynamics of reduced solutes. By extracting material accumulation rates, AD can be used to compare lakes of different trophic status (Fig. 3). Because P release is partially driven by anoxia, AD can also predict phosphorus concentrations in lakes (Fig. 2). Once AD is shown to work with other by-products of microbial anaerobic metabolism, it will provide a unique tool to assess whole-hypolimnion anaerobic metabolism.

Compound accumulation rates during anoxia
We examined the progression of anaerobic processes in lake hypolimnia over time and derived metabolite accumulation rates using AD. AD is in essence a space-for-time transformation that tracks compounds accumulation independent from water column position. Because water layers differ in their sediment-to-water ratio, we expected a stronger coupling of AD with pelagic-driven anaerobic metabolites than with those related to sediment processes. The fluorescent component F 365/480 that increased during lake Mendota hypolimnetic anoxia is a dominant fluorescence signature in anoxic marine waters (Loginova et al. 2016) which suggests the production of recalcitrant DOM from both pelagic processing and sediment release (Dadi et al. 2017;Lau and del Giorgio 2020). AD revealed two stages of anaerobic F 365/480 transformation: a month-long stabilization followed by a strong increase of this component. This behavior could be related to the change of electron acceptors, but it remains unexplained by current conceptualization and would be missed without AD inspection, hence offering a new lens to investigate carbon dynamics.
We also found good, positive relationships with TP, SRP, and NH 4 + , three compounds mostly related to sediment processes (Carey et al. 2022), and different relationships with OD. The close phosphorus concentrations between the end of OD and early AD values are not necessarily expected as they are from shallow and deep hypolimnetic depths, respectively. For Stechlin, the relationship with NH 4 + was quite dynamic.
For instance, NH 4 + accumulated during the oxic period in the deepest layers, possibly related to zooplankton release or inhibition of oxidation. The accumulated NH 4 + was rapidly consumed during the first 3 weeks prior to anaerobic conditions, presumably by nitrifying bacteria. In late 2019 early 2020, NH 4 + decreased under anoxia, as did SRP, suggesting a slow reaeration process where phosphorus and reduced compounds were oxidized and prevented DO accumulation, thus prolonging anoxia until mid-winter.   Fig. 2 (open) or their mean (full). Note that there was no oxic phosphorus accumulation in lake Delavan (slopes not statistically different from 0, panel (a)) and that the accumulation rates were based on TP for Delavan vs. SRP for Arendsee and Stechlin. This can lead to different material responses in anaerobic conditions and should be interpreted with caution as TP may include more material than SRP. Letters above the points represent statistically significant different distributions according to ANOVA and Tukey test.
In Arendsee, there was an oxic release of SRP that continued, stronger, under anaerobic conditions, when organic phosphorus diagenetic release was unobstructed by SRPcapturing iron oxides at the sediment surface (Hupfer and Lewandowski 2008). Classically, this process should be more pronounced in water layers with a higher α(z), yet we found a linear increase with depth-independent AD. Similarly, NH 4 + production is dominated by sediment diagenesis, but also a product of pelagic microbial dissolved organic nitrogen turnover (Berman et al. 1999) and accumulates in the absence of oxygen. Indeed, NH 4 + was exponentially decreasing with OD and abruptly increased once anaerobic conditions started, as indicated with the relationship with AD. The depth-independent behavior of anaerobic metabolites may be caused by limited molecular diffusivity, akin to oxygen consumption in sediments (Rippey and McSorley 2009). We estimated a lake's ADs based on a one-dimensional conceptualization of the water column, assuming that horizontal diffusivity greatly exceeds vertical diffusivity (Quay et al. 1980). Larger vertical diffusivity such as during seiching and horizontal turbulence promote vertical exchange among water layers, decreasing between-layer concentration differences (Rippey and McSorley 2009). Hence, the observed rates calculated herein would be slightly underestimated and therefore represent conservative estimates. For these reasons, this space-for-time analysis can capture the lake-scale substrate behavior under anoxic conditions.
Oxygen: From observations to modeled profile series We presented various equations to calculate AD from widely available oxygen monitoring data. The exponentialplateau provided the best fit between J z and α(z) for most lakes and years, and the best prediction of first day of hypoxia. This comparison assumed a monotonic DO decrease even though physical mixing may oxygenate upper hypolimnetic strata during stratification (Burns 1995) which would reset AD to zero. A reset of AD would be reflected in lake water composition as reduced compounds are oxygen-sensitive and rapidly oxidize (Bauer and Kappler 2009), as observed in Stechlin. This simplification of monotonous decline seemed adequate for deep lakes' hypolimnia but may need critical evaluation in shallower and more wind-exposed lakes. We note that an unconstrained exponential-plateau equation is particularly sensitive to J z at large α(z). However, constraining the equation with a priori knowledge of the system narrowed model results, then similar to the more robust log-linear equation. In addition, less comprehensive sampling scheme, such as using only three discrete loggers or two profiles, are also adequate to model J z . Common lake sampling practices therefore allow daily oxygen profiles modeling and we recommend using the most plausible fit as the objective is to accurately predict ADs.

Future perspectives
Anoxia is pervasive and of growing concern in aquatic ecosystems worldwide (Rabalais et al. 2010), promoted by various anthropogenic activities including eutrophication and browning (Brothers et al. 2014;Jenny et al. 2016). We argue that AD can be used across aquatic ecosystems to predict critical consequences of anoxia, but this remains to be proven. Indeed, AD can be tuned to specific DO threshold of interest (Table S8) to analyze the effects of anoxia on organisms (Elshout et al. 2013), greenhouse gases (Richardson et al. 2009), and toxic substances (Ach a et al. 2018). As lakespecific production rates can be modeled using AD from limited observations, these production rates will provide valuable information to study drivers and trends of anaerobic metabolism and aid in assessing aquatic ecosystems health under global change.