Modeling sediment oxygen demand in a highly productive lake under various trophic scenarios

Hypolimnetic oxygen depletion in lakes is a widespread problem and is mainly controlled by the sediment oxygen uptake (SOU) and flux of reduced substances out of the sediments (Fred). Especially in eutrophic lakes, Fred may constitute a major fraction of the areal hypolimnetic mineralization rate, but its size and source is often poorly understood. Using a diagenetic reaction-transport model supported by a large data set of sediment porewater concentrations, bulk sediment core data and lake monitoring data, the behavior of Fred was simulated in eutrophic Lake Baldegg. Transient boundary conditions for the gross sedimentation of total organic carbon and for hypolimnetic O2 concentrations were applied to simulate the eutrophication and re-oligotrophication history of the lake. According to the model, Fred is dominated by methanogenesis, where up to70% to the total CH4 is produced from sediments older than 20 years deposited during the time of permanent anoxia between 1890 and 1982. An implementation of simplified seasonal variations of the upper boundary conditions showed that their consideration is not necessary for the assessment of annual average fluxes in long-term simulations. Four lake management scenarios were then implemented to investigate the future development of Fred and SOU until 2050 under different boundary conditions. A comparison of three trophic scenarios showed that further reduction of the lake productivity to at least a mesotrophic state is required to significantly decrease Fred and SOU from the present state. Conversely, a termination of artificial aeration at the present trophic state would result in high rates of organic matter deposition and a long-term increase of Fred from the sediments of Lake Baldegg.

interface and its consumption within the top sediment layer) and the flux of reduced substances out of the sediment (F red ) [1][2][3]. The SOU is dominated by the rapid aerobic mineralization of freshly deposited organic matter (OM) and contributions from the oxidation of reduced substances such as CH 4 and NH 4 + . In lakes with a high nutrient load, O 2 levels usually decrease to zero within the first millimeters of the sediment [4], while in oligotrophic lakes the O 2 penetration depth can reach up to several centimeters [5]. The ensuing anaerobic degradation processes generate reduced substances such as CH 4 , NH 4 + , Mn(II), Fe(II) and S(-II), which once in contact with O 2 , will be oxidized. Hence, the combined fluxes of reduced substances (F red ) from the sediment can be regarded as a negative O 2 flux. Especially in the hypolimnion of shallow lakes, F red can constitute the main sink for ambient O 2 with up to 80% of the total hypolimnetic O 2 consumption [2,3].
Artificial aeration systems were installed in many anoxic lakes to increase hypolimnetic O 2 levels. Surprisingly, SOD increased in the pursuit of artificial aeration, presumably because the additional O 2 was consumed to mineralize large parts of recently deposited OM and to oxidize F red generated by "legacy" OM [3,6], which was largely buried during hypolimnetic anoxia. This phenomenon where sediments continue to be a major sink for DO, in spite of decades of artificial aeration, is called "sediment memory effect" [7]. However, it is still unclear under which circumstances hypolimnetic O 2 depletion rates will start to decrease as the relative importance of SOU and F red and their possible resilience to re-oligotrophication measures are poorly understood. To better understand these feedbacks and to navigate re-oligotrophication campaigns and aeration systems, a multi component early diagenesis model has been developed on the basis of the Van Cappellen and Wang [8], Couture, Shafei [9] and Dittrich, Wehrli [10] models. The model uses kinetic rate constants for biogeochemical reactions and processes coupled with a one-dimensional transport model using parameters mostly universal for the aquatic environment. This type of reaction network reflects the current state of the art of modeling [9][10][11] and enables the investigation of both O 2 sinks, F red and SOU. The model also allows simulations of the transient behavior of the O 2 fluxes under varying boundary conditions due to seasonal variations or long-term trends such as eutrophication, re-oligotrophication or climate change.
This paper investigates the controls of long-term O 2 depletion in the hypolimnion of the Swiss eutrophic Lake Baldegg. To this end, we model SOD of four scenarios for variations in total organic carbon (TOC) gross sedimentation rate and O 2 concentrations, and depict the resulting fluxes of reduced substances out of the sediment and the flux of O 2 into the sediments. The different modeling scenarios are based on sediment and sediment porewater data from Lake Baldegg. Lake monitoring data from Lake Baldegg was used to define the upper boundary conditions for the "Status Quo" (SQ) and "No artificial aeration" (NoAa) scenarios. Data from neighboring mesotrophic Lake Hallwil and from oligotrophic Lake Aegeri were utilized for the "Mesotrophic" and the "Oligotrophic" scenario, respectively.
In a first step, a model with constant boundary conditions was applied to reproduce average sediment porewater concentration. As SOD is governed by the quality and quantity of TOC deposited at the sediment surface and the overlying O 2 concentrations [1,2,12], the different transient input parameters of the modeling scenarios show changes in SOD and help identifying lake restoration goals. Hence, in a second step, as lake primary production and other parameters such as hypolimnetic O 2 , NO 3 and SO 4 2-concentrations and concentrations of reduced substances at the sediment-water interface (NH 4 + , Mn(II), Fe(II)) can show seasonal variations [13], their respective upper boundary concentrations were modeled with simplified cosine functions to assess the possible impact of seasonal variations on long term O 2 depletion.

Study site
The study was conducted with data from eutrophic Lake Baldegg located on the Swiss plateau. The lake with an area of 5.2 km 2 has a maximum depth of 66 m, a mean hypolimnion depth of 27.6 m and water residence time of 4.3 years [14]. The catchment of 73 km 2 is dominated by intensive agriculture and pig farms with a high surplus of manure. In 1982, the total phosphorus (TP) concentration reached 520 μgP L -1 . The lake is still eutrophic today but has recovered to *25 μgP L -1 . In spite of artificial aeration with O 2 during the stratified season and forced circulation in winter, the bottom waters become sub-to anoxic towards the end of summer stagnation. Although P loads and TP concentrations have dramatically decreased over the last 35 years, primary production and TOC gross sedimentation rate did not show a downward trend [3]. At the deepest point, anoxic conditions as evidenced by the formation of varves occurred from 1885 until the start of the artificial aeration system in 1982 [15].

Sediment data
We used a sediment porewater data set collected during two years of extensive field campaigns in 2014 and 2015. Sixteen sediment cores were retrieved for sediment porewater analysis (CH 4  , dissolved inorganic phosphorus (DIP)) and several additional cores were taken for bulk sediment analysis such as dating ( 210 Pb, 137 Cs, varve counting), TOC content, total nitrogen (TN) total phosphorus (TP), and total sulfur and physical parameters such as water content and porosity [2]. The sediment cores were taken at the deepest part of the lake at 64 m water depth (47˚11´53´´N, 8˚15´36´´E), four times during a course of a year (March, June, September, November). Details of the sediment core retrieval and analytical procedures are given in Steinsberger, Schmid [2].
Total Fe and Mn contents of the sediment were taken from Schaller, Moor [16]. Sediment trap material was collected bi-weekly and was analyzed for TN, TP and TOC gross sedimentation rates.

Model formulation
The general model approach is outlined in Van Cappellen and Wang [8] and Couture, Shafei [9]. A system of partial differential equations corresponding to early diagenesis equations is automatically generated in MATLAB and solved by MATLAB's built-in solver pdepe. The general, one-dimensional partial differential equation for the approximation of the temporal variations of the concentrations of soluble and solid substances in the sediments was used based on deposition, molecular diffusion, bioturbation and transformation processes [17]: Here, C i is the concentration of solid (μmol g -1 ) or solute (μmol cm -3 ) i, x is the positive downward position along the 1-D vertical domain (x = 0 corresponding to the sediment water interface (SWI)), and t is time. ∑εr i (z,t,C i , . . .) the sum of the sources and sinks of all species i, including the rates of all biogeochemical reactions producing or consuming species, as well as the non-local transport processes that remove or add dissolved species, most notably bioturbation. ε = φ for solutes and ε = 1φ for solids with φ defined as the sediment porosity. For solids, D i is the bioturbation coefficient D b (cm 2 yr -1 ). For solutes, D i is calculated as D b + D sed with D sed = D mo /(1-log(φ 2 ) the tortuosity corrected molecular diffusion coefficient (D mol ) at in situ temperature and salinity [18]. For solids and solutes, @equals to the burial velocity of the sediments (cm yr -1 ).
The upper boundary condition for solutes at the SWI (x = 0) is given by the flux through the diffusive boundary layer (DBL) above the SWI: Here, B i (t) is the concentration of the solute in the bottom boundary layer of the water column, and can be time-dependent in non-steady state simulations to reflect long-term or seasonal changes in bottom water concentrations with d DBL as the thickness of the DBL. Potential geochemical reactions in the DBL are neglected.
For solid-bound species (e.g. total organic carbon or Fe/Mn solid species), the flux continuity condition at x = 0 is: Where J i (t) is the (transient) depositional flux of any given solid-bound species. ρ b (1-φ) is the denominator ensuring consistency among the units of J i (t) (μmol cm -2 yr -1 ) and C i (μmol g -1 for solid species) with ρ b dry density (g cm -3 ).The depth dependence of the bioturbation coefficient D b was approximated from Katsev and Dittrich [19]: is the value at the SWI, τ b shows the characteristic depth half interval in which most of the reduction of D b occurs, while H is the depth of the steepest gradient of D b in the sediment [20]. D bmin is added to account for very small disturbances in the sediment and to increase numerical stability, however with a very low value (0.01 cm -2 yr -1 ).
At the lower boundary x L , zero gradients (equivalent to zero diffusive fluxes) are imposed for all solute and solid species: The grid points between the sediment surface and the maximum simulated sediment depth (here 50 cm) were calculated using a curvature function. This function increases the number of grid points in the upper sediment where production/consumption rates are highest. In the top centimeter the resolution varied between 0.04 mm and 0.13 mm. In deeper less reactive sediments, fewer grid points were used to increase simulation speed with a resolution reaching 2.9 mm at 30 cm sediment depth and 4.8 mm at 50 cm sediment depth. Modeled time steps are automatically chosen by MATLABs PDE solver according to the chosen process rates, species concentrations and tolerance levels. The system of PDEs (partial differential equations) defined in section 2.3 is solved in Matlab using the function pdepe within the temporal t o < t < t f and spatial (x o < x < x f ) domains. The pdepe function is designed to solve initial-boundary value problems consisting of systems of parabolic and elliptic PDEs in space and time. The numerical method is based on a piecewise nonlinear Petrov-Galerkin method with secondorder accuracy. This method solves the ordinary differential equations (ODE) from the spatial discretization of the PDEs, using the built-in Matlab ODE solver to obtain approximated solutions at specific times within a defined time interval. Furthermore, the Matlab code evaluates the boundary values at each time step separately and therefore enables transient boundary conditions. This feature becomes particularly useful when simulating the fate of compounds whose inputs are or can be changed by anthropogenic activity like artificial aeration or re-oligotrophication efforts e.g. a change in TOC gross sedimentation rate.

Modeled species
Twelve solute and ten solid species are modeled as listed in Table 1. SOU is calculated as the flux of O 2 into the sediment and F red as the sum of the fluxes of reduced substances out of the sediment (NH 4 + , CH 4 , Mn(II), Fe(II), S(-II)). Based on the respective redox stoichiometry, the fluxes of reduced compounds (J x ) were converted to fluxes of O 2 equivalents and denoted in gO 2 m -2 d -1 (Eq 6).
Both measured and modeled porewater fluxes were calculated from the vertical porewater concentration gradients by Fick's first law. Measured porewater fluxes were further evaluated by a one-dimensional reaction-transport model [21] that was adapted from Epping and Helder [22]. Although many studies use three OM pools (reactive, slow reactive and non-reactive) [10,19], this model uses only two OM pools (OM 1 , reactive; and OM 2 , non-reactive). While the total amount of settled organic matter or TOC can be accurately determined, the ratio of reactive to non-reactive TOC is not known. After setting the main reaction rates, the ratio of reactive to non-reactive TOC was determined by comparing simulated and measured CH 4 , NH 4 , TOC and TN concentrations. Best results were found for 2/3 reactive and 1/3 non-reactive TOC. Similarly, the boundary conditions of total Fe and Mn were estimated from sediment trap analyses assuming that the predominant fraction of these metals is reactive (F oh1 , 82% of total Fe; and M o1 63% of total Mn), while a smaller part is less reactive (F oh2 , 13% of total Fe; and M o2 8% of total Mn) with the remaining fraction being non-reactive and consequently not modeled. The relative fractions of reactive and non-reactive Fe and Mn were calibrated to match the measured Fe(II) and Mn(II) porewater profiles and total solid Fe and Mn sediment profiles measured by Schaller, Moor [23] The input of both organic matter pools is characterized by the measured elemental C (cx1): N (cy1): P (cz1): S (cs1) ratio of 106: 8: 0.5: 0.5. Other elements contained in the OM are discarded. Hence, by the turnover of one OM 1 106 carbon (cx1), 8 nitrogen (cy1), 0.5 phosphorous (cz1) and 0.5 sulfur (cs1) are liberated (see S1 Table).

Seasonal modeling
Seasonal changes of TOC gross sedimentation rate and concentrations in O 2 , NO 3 -, SO 4 2-, NH 4 + , CH 4 , Fe(II), Mn(II) were observed during the two year-long field campaign [2]. These naturally occurring seasonal variations could influence the hypolimnetic O 2 depletion rate. As the model allows for temporally varying upper boundary conditions, seasonal variations can be readily appended. However, direct data inputs of e.g. highly fluctuating TOC gross sedimentation rates captured in a daily and/or biweekly basis, result in numerical instability. Instead, the boundary conditions were multiplied with cosine functions with different amplitudes and phase shifts as simple approximations for the measured seasonal fluctuations (see Fig 1 and Table 2) [10]. The same cosines functions were applied to all scenarios leading to slight overestimations of hypolimnetic O 2 levels in the "oligotrophic scenario".

Biogeochemical reactions
All biogeochemical reactions are listed in S2 Table. The mineralization processes are calculated with the stoichiometric coefficients of the OM cx 1 (C), cy 1 (N), cz 1 (P) and cs 1   are produced. All reduced substances are prone to aerobic oxidation by secondary reactions e.g. nitrification or anaerobic oxidations such as anaerobic methane oxidation (AMO) or oxidation by other oxides such as ferric iron. Further, various more stable forms of Fe-and Mncontaining mineral phases are formed such as iron sulfides, pyrite, vivianite, and manganese carbonates.

Simulated scenarios
The model and the respective modeling scenarios are based on historic observations, lake monitoring data, bulk sediment parameters and sediment porewater data of Lake Baldegg [2,23,25,26]. To calculate the different scenarios, two other lakes of comparable size and morphometry were chosen where complete limnological data sets were available. Similar to Lake Baldegg, neighboring Lake Hallwil was aerated since the 1980s but already shows decreased primary production and was therefore chosen as representative for the mesotrophic scenario. Lake Aegeri, where total phosphors concentrations are below 10 μg L -1 due to the low fraction of agriculture in its catchment, was chosen as representative for the oligotrophic scenario. The initial conditions were set to match the presumed natural oligotrophic state of the lake with average DO concentrations of 5.76 mg O 2 L -1 available at the sediment water interface and low TOC gross sedimentation rate and were modeled for 350 years. All boundary conditions can be found in Table 1 and all reaction rates in Table 3. The transient boundary conditions of the model were then used to simulate the advancing eutrophication of Lake Baldegg starting in 1860 with decreasing O 2 concentrations until permanent anoxic conditions were reached in 1890. These anoxic conditions are preserved in the sediments by the onset of sediment varves [15,25]. The anoxic phase in the hypolimnion of Lake Baldegg lasted at least from 1890 until 1982 and was only reversed by the initiation of artificial aeration. During the anoxic phase, the boundary conditions of all reduced substances, especially for NH 4 + , rose as they started accumulating in the hypolimnetic waters and became an increasing environmental problem. Furthermore, this condition led to enhanced geochemical focusing as more reduced manganese and iron were exported to the sediments of the hypolimnion [27]. Based on the available lake O 2 monitoring data, following the onset of aeration boundary O 2 concentrations were set to 3.2 mg O 2 L -1 between 1982 and 2018. The simulations were then continued from 2018 to 2050 with four scenarios with different upper boundary conditions to represent different trophic and oxic conditions in Lake Baldegg until the year 2050: 1. Status Quo (SQ) scenario, in which O 2 concentrations (average 3.2 mg O 2 L -1 ) and TOC gross sedimentation rates (107 g C m -2 yr -1 ) remain constant at current levels.
2. Mesotrophic production (M) scenario based on lake monitoring data from neighboring mesotrophic Lake Hallwil [2], which has been undergoing a re-oligotrophication phase since the early 1980s. Similar to Lake Baldegg, lake TP concentrations in Lake Hallwil decreased from 235 μg P L -1 in 1976 to 90 μg P L -1 in the 1990s and further to 12 μg P L -1 in 2015. However in Lake Hallwil, this change was accompanied by a drastic decrease of primary production and hence TOC gross sedimentation from over 100 g C m -2 yr -1 [28] to around 60 g C m -2 yr -1 [2]. Hence, for the M scenario the TOC gross sedimentation rate was reduced by 44% of the present day value and O 2 concentrations were increased to an annual average of 4.8 mg O 2 L -1. as observed in Lake Hallwil.

No artificial aeration scenario (NoAa):
As artificial aeration systems are expensive due to maintenance costs, the reduction of maintenance costs or termination of the aeration systems are regularly debated. However, it is highly likely that without artificial aeration Lake Baldegg would rapidly turn anoxic during summer stagnation, and depending on the intensity of winter mixing, most probably become permanently anoxic within a few years, similar to its status before the onset of artificial aeration in 1982. Therefore, the O 2 concentration is set to 0 mg L -1 . Similarly, NO 3 and SO 4 2concentrations are decreased. The TOC gross sedimentation is assumed to be not affected by this and is kept constant at the present-day level. Furthermore, the boundary conditions of the reduced substances NH 4 + , CH 4 , Fe(II), Mn(II) and S(-II) are increased to 300 μmol L -1 , 100 μmol L -1 , 80 μmol L -1 , 80 μmol L -1 and 30 μmol L -1 , respectively, to account for their likely accumulation in the hypolimnion (see Table 1). The underlying premise of the SQ, M and O scenarios is that artificial aeration remains active to supply sufficient hypolimnetic DO. The detailed boundary conditions for each scenario are given in Table 1

Model calibration
The values of the model parameters (Table 3) were initially set either based on literature values or field observations. Boundary conditions were subsequently modified during the calibration process, using a manual trial-and-error approach (Table 1). Further, reaction rates taken from similar model studies by Van Cappellen and Wang [8], Dittrich, Wehrli [10] and Katsev and Dittrich [19] were used as initial rates. First, by using these initial rates, relative fractions of reactive to non-reactive pools of OM, Fe, and Mn were determined by visually comparing the modeled and measured porewater and sediment profiles. Then, process rates were calibrated to further improve, reproduce and match the measured porewater concentration profiles and sediment concentrations without using numerical fitting procedures (see asterisks Table 3).
The sensitivity of the model results to individual process rates was assessed by step-wise increasing/decreasing (±1%) the respective rate and comparing the resulting AHM to that from the baseline run. Sensitivities in the following are expressed as relative change in the model output divided by relative change of the model parameter. The sensitivity of O 2 consumption and F red production to most of the process rates was low (below <0.01), but adjusting these rates was required for reproducing the observed concentrations of individual species Based on the fluxes across the sediment-water interface, the total mass within the sediment and the amounts leaving the modeled sediment depth range at the bottom, the total mass balance of each element species was calculated. Mass balances of all species were preserved. Carbon and nitrogen mass balances are shown in S1 Fig

Porewater profiles
In agreement with sediment porewater measurements [2], modeled porewater concentrations of NH 4 + and CH 4 increased with increasing sediment depth. This was mainly a result of the aerobic mineralization at the top of the sediment and the dominating process of methanogenesis in the deeper parts of the sediment. While modeled and measured CH 4 concentrations were similar in the top few of centimeters, the deviation increased with increasing sediment depth. Reproduction of the observed shape of the CH 4 profiles required an additional CH 4 sink below the oxic zone or changes of the C/N ratio in the OM pool. Several anaerobic oxidation processes of CH 4 such as oxidation by ferric iron oxides [30], Mn(IV)-oxides [31] or humic compounds [32] might play a role here but are not implemented in the reaction network. Modeled and measured porewater profiles of Fe(II) and Mn(II) were also in reasonable agreement, despite their complex shapes due to various precipitation and dissolution reactions. The Fe(II) peak is caused by the initial reduction of iron (hydr-)oxides. The resulting ferrous iron diffuses towards the sediment surface where it is rapidly re-oxidized to various forms of ferric (hydr-)oxides creating an iron redox cycle. Although the major reactions for Fe(III) and Mn(IV) precipitation and reductive dissolution were included in the model, we acknowledge that without simultaneous pH calculations their significance remains limited. However, as Fe(II) and Mn(II) concentrations and hence fluxes only contributed a very small fraction to the total F red [3,7], their values had a negligible effect on the overall picture. concentrations of around 245 μmol L -1 in the sediment overlying water and decreasing to zero within 2 mm below the diffusive boundary layer. Similar rapid drops in porewater concentrations were observed for NO 3 - and SO 4 2- , where porewater concentrations decrease to zero within 5 mm and 15 mm, respectively.

Fluxes of reduced substances (F red )
The porewater fluxes were estimated from measured porewater concentration profiles using Fick's first law of diffusion. The fluxes where then expressed as O 2 equivalents according to their oxidation stoichiometry and summed up to F red [3]. The F red flux resulting from the sediment model for 2015 was 0.36 gO 2 m -2 d -1 . The modeled value is close to the value determined with the steady-state 1-D model by Epping and Helder [22] of 0.49±0.09 gO 2 m -2 d -1 in 2015 [2]. Likewise modeled and measured CH 4 fluxes were in perfect agreement with 0. 22

Sediment TOC
To reproduce the observed vertical profiles of TOC in the sediment, we calibrated the ratio between reactive and non-reactive OM. A good agreement between simulated and observed TOC was reached under the assumption that two thirds of the total TOC gross sedimentation are reactive with an absolute error of~7.5% between the measured and modeled TOC values (Fig 4). This assumption seems reasonable for a highly productive lake such as Lake Baldegg. A large fraction of the reactive TOC pool (OM 1 ) was mineralized directly at the sediment surface, and decreased by about 50% within the top 10 centimeters of the sediment. Deeper in the sediment, the largest TOC pool consisted of non-reactive TOC, and no reactive TOC was left below 33 cm.

Simulated scenarios
The calibrated model is an excellent tool for assessing the future development of the O 2 consumption rate in the lake applying different management strategies. The transient boundary conditions of the model were adjusted to simulate four different scenarios up to the year 2050. Each setup started with the same initial conditions in 1850, and boundary conditions were only changed after 2018 to account for the four abovementioned scenarios: "Status Quo" (SQ), Mesotrophic production (M), "Oligotrophic production" (O) and "No artificial aeration" (NoAa). The underlying premise for all settings (except NoAa) is that artificial aeration remains active to supply sufficient hypolimnetic O 2 . Furthermore, all scenarios were modeled twice, once with no seasonal variations and once with simplified seasonal variations of the changing input parameters described in Table 1.

SQ scenario (Status quo continued)
The modeled concentration profiles of NH 4 + and CH 4 for the year 2050 in the various scenarios are shown in Fig 5. In the SQ scenario, NH 4 + fluxes are projected to decrease by 8.6% from 2015 to 2050. In contrast, CH 4 fluxes (37%) decline more rapidly until 2050. NH 4 + is produced during each mineralization step of organic matter, while CH 4 is solely the result of methanogenesis, which occurs deeper in the sediment and is fueled by the buried reactive TOC (see Fig 6).
In 2015, the methanogenesis in sediments older than 20 years contributed~20% of NH 4 + production but~70% of the CH 4 production. In contrast, 68% of the NH 4 + released from the sediments originate from OM deposited during the last two years and, therefore, its contribution indirectly reflects the current primary production. By 2050, the contribution to NH 4 + and CH 4 production by sediments older than 20 years will decrease by~55%, whereas almost no change in sediments younger than 20 years occurs (see Fig 6). This suggests that the upper sediment strata will remain in quasi steady state with the current TOC gross sedimentation and O 2 concentration, while the stock of buried reactive TOC from the era of hypertrophy and anoxia decreases over time (see Fig 7). The SQ assumption shows that the present share of NH 4 + and CH 4 production from sediments older than 20 years is 53% and will decrease to 39% by the year 2050 (Fig 6) These results are in agreement with previous studies. Matzinger, Müller [7] calculated that in Türlersee and Pfäffikersee, deep buried organic rich sediments (>20 years) contributed up to 29% of the CH 4 and NH 4 + fluxes from the sediment. Total F red will decrease by 23% from 0.35 gO 2 m -2 d -1 (2015) to 0.27 gO 2 m -2 d -1 (2050) until 2050 (Fig 8A), while SOU remains unchanged at 0.36 gO 2 m -2 d -1 (Fig 8B). Hence, SOD will decrease only by 13% until 2050 (see Fig 8C), largely as a consequence of the diminishing pool of reactive TOC.

Reoligotrophication scenarios
The mesotrophic (M) and oligotrophic (O) scenarios consider the potential decrease of DO consumption if further measures for reoligotrophication are implemented in the lake and its catchment. We assumed the TOC flux to be reduced to 60 gC m -2 yr -1 for the M scenario and to 36 gC m -2 yr -1 for the O scenario. Simultaneously, average O 2 concentrations are increased from 3.2 mgO 2 L -1 (in SQ) to 4.8 mgO 2 L -1 (in M) and to 6.4 mgO 2 L -1 (in O) over a 12-year period from 2018 to 2030. These changes lead to decreased NH 4 + and CH 4 production rates in the sediments (Fig 6). Sediment CH 4 production decreases by~58% for both the M and O scenarios compared to the SQ scenario. Interestingly, the eventual CH 4 production is not smaller for the O than for the M scenario, probably because already in the mesotrophic state, TEA concentrations are sufficient to mineralize enough TOC and thus minimize the supply for methanogenesis. In contrast to the present state and the SQ scenario, mineralization of reactive TOC is largely completed within the top few centimeters and only a small fraction of reactive TOC is available for methanogenesis deeper in the sediment (see Figs 6 and 7). In both reoligotrophication scenarios, the main fraction of CH 4 (88%) originates from sediments older than 20 years. NH 4 + production decreases by 38% in the M scenario and by 59% in the O scenario.
In contrast to CH 4 production, the largest decrease in NH 4 + production occurs in the top two centimeters of the sediment and reflects the decreased deposition and increased initial mineralization of OM (see Fig 6).  While SOU is predicted to increase from 0.36 gO 2 m -2 d -1 in 2015 to 0.38 gO 2 m -2 d -1 by 2050, F red decreases from 0.36 gO 2 m -2 d -1 to only 0.08 gO 2 m -2 d -1 in the M scenario, to 0.05 gO 2 m -2 d -1 in the O scenario and only contributes 12% and 18% to SOD, respectively (see Fig  8A and 8B). Hence, in the reoligotrophication scenarios SOD decreases by over 36% and is mainly controlled by SOU and, therefore, by the lake's current primary production.

No artificial aeration scenario (NoAa)
Without measures to reduce primary production and without artificial hypolimnetic aeration promoting aerobic mineralization of TOC, large amounts of reactive TOC would be buried in the sediments of Lake Baldegg (see Fig 7). Under the resulting anoxic conditions bioturbation rapidly decreases (see Table 1). The high deposition of reactive TOC (see Fig 6) increases the production of NH 4 + and CH 4 in all sediment depths (see Fig 5). As a result, methanogenesis dominates TOC mineralization in 2050 with up to 93% in contrast to only about 30% in 2015. Furthermore, fluxes of NH 4 + and CH 4 increase by 23% and 191%. In total, F red doubles from 0.36 gO 2 m -2 d -1 to 0.72 gO 2 m -2 d -1 by 2050 (see Fig 8A). The projected increase of the NH 4 + flux from the sediments is rather modest due to the increasing concentration of NH 4 + in the bottom water (300 μmol L -1 ). The boundary condition for CH 4 was set to 100 μmol L -1 , however, we acknowledge that much higher concentrations might be possible. In Lake Rotsee, a highly productive and eutrophic lake, CH 4 concentrations reach the gas saturation concentration directly below the sediment surface during summer stagnation, and CH 4 is therefore additionally lost to the hypolimnion by bubble formation [33]. SOU obviously is zero during permanent anoxia (see Fig 8B), and SOD consists only of F red (0.72 gO 2 m -2 d -1 ). However, once aeration is restarted, a large spike in SOD, similar to the spike following the begin of aeration in 1982, is expected to occur. The large TOC pool buried during the anoxic phase (2020-2050) would cause high F red values for decades (similar to the reactive TOC pool buried during the anoxic phase between 1890 to 1982). This would impede any efforts of ensuing re-oligotrophication measures and artificial aeration.

Seasonal modeling
In the seasonal models, concentrations and fluxes of TOC, NH 4 + , CH 4 , Mn(II) and Fe(II) increased during summer stratification, while concentrations of O 2 , NO 3 and SO 4 2decreased (see Fig 1). Consequently, during the time of the highest TOC deposition rate, the lowest concentrations of TEAs occurred. This resulted in the highest fluxes of F red at the end of summer stratification (Fig 9) as increasingly more material is degraded anaerobically at that time. Lowest F red values occurred during winter overturn, where highest O 2 concentrations and lowest TOC gross sedimentation prevailed. Hence, during winter overturn the deposited TOC is increasingly exposed to and mineralized by O 2 thereby diminishing the storage of OM in the sediment. Since most of the CH 4 was produced in the deeper parts of the sediment, its rate was independent of seasonal variations (Fig 10). It decreased only when old deposits of reactive TOC were exhausted. However, seasonal variations of CH 4 fluxes are induced by varying O 2 penetration depths (see Fig 8). In contrast to CH 4 , only~30% of the NH 4 + is produced during methanogenesis (Fig 10). Aerobic mineralization of OM (� 37%) and dissimilatory reduction of nitrate to ammonium (� 27%) together produce up to 64% of NH 4 + . As O 2 concentrations in the hypolimnion decrease during stratification, the production of NH 4 + via aerobic mineralization diminishes likewise. The seasonal impact of other TEAs such as NO 3 -, Fe(III)-(hydr) oxides, Mn(IV)-oxides and SO 4 2on the overall NH 4 + production is negligible. Conversely, the fraction of aerobically degraded TOC increases at the end of each year. In contrast to CH 4 , both the production and fluxes of NH 4 + show distinct seasonal variations (Fig 10).
Seasonal variations of fluxes of reduced substances were observed in several lakes [2,7,34]. The coincidence of highest TOC gross accumulation and lowest TEA concentrations in the middle of any given year leads to higher accumulation rates of reactive TOC and thus to increased burial of reactive TOC in the sediments (see Fig 11) [12]. However, this effect is only apparent in the seasonal models.
As more reactive TOC is buried, slightly higher fluxes of F red are encountered in the seasonal models compared to non-seasonal models (see Fig 9).  In conclusion, simplified seasonal variations of the upper boundary condition cause seasonal variations in the NH 4 + production, slightly higher accumulation rates of reactive TOC, variations in F red , increased SOU and consequently, slightly higher O 2 depletion rates. Although the overall impact of seasonal boundary variation is negligible, they are important to explain observed porewater profiles of species such as NH 4 + , Fe(II) and Mn(II), which have been shown to vary seasonally, especially in the top five centimeters of the sediments (see Fig 2)[2].

Conclusion
A multicomponent early diagenesis model with transient boundary conditions was developed and calibrated using historic observations, monitoring data and recent state of the art porewater data from the eutrophic and artificially aerated Lake Baldegg. The model was applied to simulate four different management scenarios, including a continuation of the status quo, two re-oligotrophication scenarios, and a scenario where lake aeration is stopped. Simulations with seasonal boundary conditions yielded similar O 2 depletion rates as the non-seasonal model, implying that it is not necessary to consider seasonality and thus increase computation time for long-term projections of sediment oxygen demand.
All modeled scenarios show the direct impact of the TOC gross sedimentation rate and the bottom water O 2 concentration on the hypolimnetic O 2 depletion rates. While the upper boundary O 2 concentration is doubled between the scenarios (SQ to O), the sediment oxygen uptake (SOU) differs only slightly. However, less reactive TOC is buried which significantly decreases F red . Although older sediments (deposited during the hypertrophic and anoxic period) dominate F red due to methanogenesis, this "sediment memory effect" [7] or "legacy carbon effect" [6] is overcome once primary production decreases. This is shown in both reoligotrophication scenarios, where SOD is mostly controlled by SOU. Interestingly, already a change to mesotrophic production yields a SOD rate comparable to a low productive state. Hence, with an adequate O 2 concentration, a return to mesotrophic productivity could be sufficient for a long-term decline of SOD. In contrast, a termination of the aeration program without measures to decrease TOC gross sedimentation would yield long-lasting damaging effects on the O 2 consumption in the hypolimnion of Lake Baldegg as large amounts of buried reactive TOC would lead to a high production of CH 4