OMEN-SED(-RCM) (v1.1): A pseudo reactive continuum representation of organic matter degradation dynamics for OMEN-SED

OMEN-SED is a one-dimensional, analytical reaction-transport model for early diagenesis in marine sediments. It explicitly resolves organic matter (OM) degradation and associated biogeochemical terminal electron acceptor, reduced species and nutrient dynamics in porous media under steady state conditions. OMEN-SED has been specifically designed for the coupling to global Earth System Models and the analytical solution of the coupled set of mass conservation equations ensures the computational efficiency required for such a coupling. To find an analytic solution, OMEN-SED expresses all 5 explicitly resolved biogeochemical processes as a function of OM degradation. The original version of OMEN-SED contains a relatively simple description of OM degradation based on two reactive OM classes, a so-called 2-G model. However, such a simplified approach does not fully account for the widely observed continuous decrease in organic matter reactivity with burial depth/time. The reactive continuum model that accounts for the continuous distribution of organic compounds over the reactive spectrum represents an alternative and more realistic description, but cannot be easily incorporated within the 10 general OMEN-SED framework. Here, we extend the diagenetic framework of OMEN-SED with a multi-G approximation of the reactive continuum model (RCM) of organic matter degradation by using a finite, but large number of OM fractions, each characterized by a distinct reactivity. The RCM and its multi-G approximation are fully constrained by only two free parameters, a and ν, that control the initial distribution of OM compounds over the reactivity spectrum. The new model is not only able to reproduce observed pore water profiles, sediment-water interface fluxes and redox zonation across a wide range 15 of depositional environments, but also provides a more realistic description of anaerobic degradation pathways. The added functionality extends the applicability of OMEN-SED to a broader range of environments and timescales, while requiring fewer parameters to simulate a wider spectrum of OM reactivities. 1 https://doi.org/10.5194/gmd-2021-4 Preprint. Discussion started: 8 March 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
The degradation of organic matter (OM) in marine sediments is a key component of the global carbon cycle and climate (Berner, 1991;Archer and Maier-Reimer, 1994;Ridgwell and Zeebe, 2005). It is the main engine behind the complex and dynamic network of biogeochemical reactions in marine sediments and, thus, controls benthic carbon and nutrient recycling, as well as carbon burial (e.g. Arndt et al., 2013;LaRowe et al., 2020b). As a consequence, it exerts an important influence 5 on marine primary production, the strength of the primary long-term sink for atmospheric CO 2 and, ultimately, the size of the largest carbon reservoir within the Earth system (e.g. Berner, 1980;Soetaert et al., 2000;Aller, 2014;LaRowe et al., 2020a). Understanding, quantifying and predicting organic matter degradation and associated biogeochemical dynamics in marine sediments is thus critical for our ability to understand past, present and future biogeochemical cycles and climate dynamics. 10 Benthic biogeochemical dynamics are driven by a complex and dynamic interplay of transport and reaction processes that operate over different temporal and spatial scales. This underlying complexity compromises our ability to understand, quantify and predict diagenetic dynamics. In this respect, reaction-transport models (RTMs) that account for transport (advection, molecular diffusion, bioturbation, and bioirrigation) and reaction (production and consumption) processes, are, in combination with observational data, powerful analytic, diagnostic and predictive tools. Because of its prominent role, the description 15 and parametrization of OM degradation in diagenetic models is critical for their ability to capture and predict sediment-water exchange and burial fluxes. The rate of OM degradation and, thus, the intensity of the associated biogeochemical cycling is first and foremost driven by the supply of OM at the sediment-water interface. However, it is not only the quantity of that flux that determines benthic biogeochemical rates, but also its quality (Middelburg et al., 1993;Arndt et al., 2013;LaRowe et al., 2020b). In particular, the benthic redox zonation and benthic recycling fluxes are largely driven by the susceptibility of OM to 20 microbial degradation (Freitas et al., 2020).
Existing OM degradation models can be broadly divided into two different classes according to their description of apparent reactivity of the bulk OM, k(t), and its evolution over time which depends on a complex and dynamic interplay of many different factors (Arndt et al., 2013). Discrete models, also called multi-G models, divide the bulk OM into several discrete fractions, OM i , each characterized by a different reactivity, k i (Berner, 1980;Jørgensen, 1978). The apparent reactivity of the 25 bulk OM, k(t), and its evolution over time is then given by the weighted sum of the individual reactivities k i : OM is degraded according to first-order degradation kinetics and the change in OM concentration with burial time/depth is thus given by: (RCM) apply a gamma distribution because of its mathematical flexibility and its ability to capture observed OM degradation dynamics (Boudreau and Ruddick, 1991). Organic matter degradation can then also be formulated as a first-order degradation rate law with an apparent OM reactivity that decreases with burial time/depth, k(t), (Boudreau and Ruddick, 1991): with: where a and ν are free, positive parameters that control the shape of the initial distribution of OM compounds over the reactivity spectrum and, thus the evolution of reactivity with burial time/depth.
The application of continuum models in the framework of local, one-dimensional diagenetic models is straight forward 15 as these models generally solve the coupled reaction-transport equations numerically. However, computational power rapidly becomes a limiting factor if a coupled set of reaction-transport equations were to be numerically solved on a global scale and over longer timescales. Therefore, global diagenetic models designed for large scale applications and the coupling to Earth System Models are generally either highly simplistic, limited in scope or rely on the analytical solution of the reaction-transport equation (Hülse et al., , 2018. However, finding an analytical solution for complex reaction networks and depth-varying 20 parameters is difficult and requires several simplifying assumptions. For instance, the recently developed one-dimensional and numerically efficient diagenetic model, Organic Matter ENabled SEDiment model (OMEN-SED) assumes steady-state conditions and accounts for depths changes in parameters and reaction processes by dividing the sediment into a set of discrete, yet dynamic biogeochemical zones (Hülse et al., 2018). The considered diagenetic processes are, either expressed as a function of OM degradation or as zero-or first-order reactions. In the original version of OMEN-SED, bulk OM degradation is described 25 by a multi-G model (i.e. as a sum of first-order degradation rate laws). Consequently, related biogeochemical processes, such as different metabolic pathways, can be expressed as a series of exponential terms (i.e. the analytical solution of the OM conservation equation), for which a general analytical solution can be easily found (Hülse et al., 2018). However, the key limitation of such a multi-G type approach is that it converges to a constant apparent OM reactivity at depth and thus fails to capture the widely observed continuous decrease in OM reactivity with burial time/depth (Middelburg and Meysman, 2007).

30
Although this limitation can be somewhat mitigated by introducing a larger number of compound classes, each new compound class i requires constraining two new parameters (k i and OM i ). Yet, because of the difficulty associated with identifying discrete classes based on available observation data, the total number of classes that can be reasonably well constrained is typically restricted to a maximum of three (Jørgensen, 1978;Middelburg, 1989). In addition, due to the lack of a theoretical framework that would allow to constrain OM reactivity on a global scale, the choice of OM degradation model parameters k i and OM i is particularly challenging for global scale applications (Hülse et al., 2018). Reported degradation rate constants in marine sediment have been shown to vary by about 10 orders of magnitude (Middelburg et al., 1993;Arndt et al., 2013). Their spatial variability on the global scale, as well as their response to changing environmental conditions over, for instance, past 5 extreme climate events or to projected climate change is largely unknown. Furthermore, the timescales of OM degradation in the water column (days to weeks) are vastly different than those found in marine sediments (10 -10 6 years). Consequentially, model parameters of OM degradation cannot be directly inferred from the pelagic parametrization. The selected, discrete OM pools of the pelagic model are optimized for the representation of OM degradation dynamics on short-time scales and can thus not be applied to describe the slower, long-term OM degradation dynamics in the underlying sediment (Hülse et al., 2018). 10 Continuum models, on the other hand, merely require constraining the two free parameters that define the shape of the initial distribution and the continuous decrease of OM reactivity with time/depth. Therefore, the RCM approach captures a wider range of OM reactivity scenarios over both short and long timescales, while using fewer parameters than multi-G models (Middelburg, 1989;Jørgensen, 1978). In addition, the RCM is in better agreement with our theoretical and qualitative understanding of the OM degradation process (Aller and Blair, 2006). The application of a RCM approach in a global diagenetic 15 model designed for large scale applications and the coupling to Earth System Models would ease model parametrization requirements. Yet, the RCM cannot be directly applied in OMEN-SED since OMEN-SED's analytical solution relies on the exponential decrease of OM concentration (see Hülse et al., 2018, for details).
Therefore, we here developed a multi-G approximation of the RCM and directly integrated it into the mathematical framework of OMEN-SED. We first provide a short summary of OMEN-SED, including a description of the general model approach 20 and the generic algorithm used to match internal boundary conditions and to determine the integration constants for the analytical solutions. We provide a detailed description of the newly integrated RCM approximation for organic matter degradation, followed by an evaluation of the performance of OMEN-SED-RCM by 1) comparing simulated depth profiles of OM, terminal electron acceptors (TEAs) and metabolic by-products with observations from selected sites, as well as 2) benthic exchange   Berner, 1980;Boudreau, 1997): where C i is the concentration of biogeochemical species i, ξ equals the porosity φ for solute species (which is defined as φ = volume porewater volume porewater + solid sediment ) and (1 − φ) for solid species. The term z denotes the sediment depth, t is the time, w is the sedimentation rate, D i is the apparent diffusion coefficient of dissolved species i, and j R j i represents the sum of all biogeochemical rates j affecting species i.
OMEN-SED accounts for both the advective and diffusive transport of dissolved and solid species. Sedimentation is sim-10 ulated as a constant rate, w, and neglects the effect of sediment compaction (i.e. dφ dz = 0) due to mathematical constraints (see Hülse et al., 2018). Fick's law of molecular diffusion describes the transport of dissolved species with a species-specific apparent diffusion coefficient D mol,i . The effect of bioturbation is parametrized by a diffusive term (Boudreau, 1996) with a constant bioturbation coefficient D bio within the bioturbated zone (z ≤ z bio ). The reaction network implemented in OMEN-SED accounts for the most important primary and secondary redox reactions, equilibrium reactions, mineral dissolution and 15 precipitation. Adsorption and desorption processes associated with OM degradation dynamics that affect dissolved and solid species are also explicitly resolved in the model. The resulting set of non-linear, depth-varying coupled equations can per se only be solved numerically. To find an analytic solution, OMEN-SED assumes steady-state conditions and expresses all explicitly resolved biogeochemical processes as a function of OM degradation, as zero-or first-order reactions. In addition, it divides the sediment into a bioturbated (z < z bio ) and a non-bioturbated zone (z bio ≤ z), as well as an superimposed set of 20 discrete, yet dynamic redox zones (Hülse et al., 2018): 1) an oxic zone (0 < z ≤ z ox ); 2) a denitrification (or nitrogenous) zone (z ox < z ≤ z NO3 ); 3) a sulfate reduction zone (z NO3 < z ≤ z SO4 ); and finally, 4) a methanogenic zone (z > z SO4 ) (Fig. 1).
OMEN-SED allows bioturbation to occur in the anoxic zones of the sediment (here all zones z > z ox combined), and so does not limit the maximum bioturbation depth to be restricted to a specific redox zone. The extent of each redox zone is dynamic, and the depth of each redox boundary varies in response to changes in ocean boundary conditions and forcing. The zones are 25 interlinked by imposing continuity in concentrations and fluxes at the dynamic internal boundaries (z b ∈ z bio , z ox , z NO3 , z SO4 , (compare e.g. Ruardij and Van Raaphorst, 1995;Billen, 1982)). The general equation (eq. 5) is applied to each redox zone with depth invariant parameters, but parameters and the formulation of the reaction term ( j R j i ) in Eq. 5 can vary between the redox zones. OMEN-SED only resolves the most pertinent reaction process within each redox zone and, thus, simplifies the mathematical description of the reaction network, while retaining its biogeochemical complexity (Hülse et al., 2018). accounts for a continuous, yet dynamic distribution of OM compounds over a range of reactivities and captures the widely observed, continuous decrease in apparent OM reactivity with depth/burial age (Boudreau and Ruddick, 1991, and see Fig. 1).
Within the RCM the rate of organic matter degradation, R OM , is given by: where om(k, t) denotes a probability density function that determines the amount of bulk OM characterized by a reactivity 5 between k and k + dk at time t. The initial distribution of organic compounds, om(k,0), may take different mathematical forms, but cannot be inferred by observations, but a gamma function is often used due to its mathematical properties Boudreau and Ruddick (1991), following Aris (1968) and Ho and Aris (1987). Assuming first order degradation kinetics, the initial (t = 0) distribution of om over k is given by: where OM (0) is the initial OM content (at the sediment water interface, SWI), Γ is the gamma function, a (yr) is the average lifetime of the more reactive components of the mixture and ν is a dimensionless parameter determining the shape of the distribution near k = 0. The free parameters a and ν completely determine the shape of the initial OM compound distribution over the range of k. They thus control the overall reactivity of the bulk OM and its evolution with depth/burial 15 age. High ν and low a values indicate the presence of highly reactive compounds that get rapidly degraded, leading to a rapid decrease in OM reactivity with depth/time. Low ν and high a values, on the other hand, indicate a dominance of less reactive compounds within the initial distribution, resulting in a lower bulk OM reactivity and a slower decrease in reactivity with depth. Although the choice of the gamma function that describes the initial distribution of OM compounds is partly guided by mathematical experience, it has been shown to capture the widely observed dynamics of OM degradation in marine sediments 20 well (Arndt et al., 2013). Assuming steady-state conditions (i.e. ∂OM ∂t = 0), the depth profile of bulk OM concentration is given by Boudreau and Ruddick (1991): where OM (0) is the known concentration at the SWI and t(z) refers to the age of the sediment layer at depth z. 25 Because the OMEN-SED framework relies on an exponential decrease of OM concentration with depth and requires invariant first-order degradation rate constants, we here approximate the reactive continuum model by a multi-G description of OM degradation (see Fig. 1b). In other words, bulk organic matter is represented by large number of i distinct compound classes, OM i , each degraded according to a first-order degradation rate constant, k i . The initial fraction of organic matter in compound class i, OM i ,and its respective degradation rate constant, k i , can be calculated based on the initial probability density function 30 by determining the amount of OM within a given reactivity range k and k + dk at t = 0 (Eq. 6). Because the initial probability 7 https://doi.org/10.5194/gmd-2021-4 Preprint. Discussion started: 8 March 2021 c Author(s) 2021. CC BY 4.0 License. density function is fully described by the two free parameters a and ν, OMEN-SED-RCM only requires defining two parameters instead of the 2 cot i parameters of the multi-G description.
Based on Eq. 7, the initial amount of OM within a discrete reactivity class k can be calculated as: The initial fraction of OM within the reactivity range between 0 and k, i.e. having a reactivity ≤ k at t = 0, is then given by integrating Eq. 9, assuming a, ν, k > 0: where Γ(ν, a · k) denotes the inverse gamma function and Γ(ν, 0) = Γ(ν). To ensure a full coverage of the relevant reactivity space, we recommend choosing a reactivity range from k min = 10 −15 to k max = 10 −log(a)+2 , although the lower k limit can be adapted to the timescales resolved. In addition, the total reactivity range should be divided into a least i=100 equal reactivity bins to ensure an appropriate approximation of the initial OM distribution in case parameter a « 0.01 (Fig. 2). However, 15 because computational cost is not a limiting factor for OMEN-SED-RCM, a larger number of classes can also be applied.
Here, the RCM is approximated by dividing the reactivity range k = [10 −15 , 10 −log(a)+2 ] into n=500 equal reactivity bins, each characterized by a different reactivity constant, k i , thus ensuring a comprehensive approximation of the gamma function defined by the respective a and ν values.
The reactivity bins are defined by k i and k i +dk i at t(0): The least and most reactive fraction, F min with k min = 10 −15 yr −1 and F max with k max = 10 −log(a)+2 yr −1 , are calculated based on the lower and upper incomplete gamma function, respectively (LaRowe et al., 2020a): The derived degradation rate constants for OM i , k i , are then used in the first-order reaction term within the conservation equation for organic matter dynamics (see Hülse et al. (2018) for details about the analytical solution): with D OMi = D bio for z ≤ z bio and D OMi = 0 for z > z bio and OM i = F i · OM (0). cores that were used to benchmark OMEN-SED (Hülse et al., 2018): Santa Barbara Basin (Reimers et al., 1996), Iberian margin and the Nazaré Canyon (Epping et al., 2002). These simulations cover a wide range of different benthic environments reaching from the shallow shelf to the deep-sea, including the highly productive eastern upwelling area on the Iberian Shelf (108 m), the upper slope, organic-rich and non-bioturbated sediments underlying anoxic bottom waters in the Santa Barbara Basin (585 m) and lower slope on the Iberian margin (2213 m), as well as the deep-sea sediments in the Nazaré Canyon (4298 m). For a 10 detailed description of the study sites, see Reimers et al. (1996) and Epping et al. (2002). The respective boundary conditions used in the simulation are provided in Table 1. Additional sediment characteristics (e.g. sedimentation rate, porosity, density), stoichiometric factors and secondary reaction parameters are chosen as in Hülse et al. (2018) (see Tables B1 and C1). The free RCM parameters a and ν are the optimized to find the best fit between observations and model results, while all other parameters are kept at their standard value (Hülse et al., 2018).  More specifically, for the two Iberian margin sites (108 m and 2213 m) OMEN-SED-RCM fits observations generally well. A slight overestimation ofNH + 4 production is simulated for the shallower site. A good agreement between simulation results and observations is also found for the Nazaré Canyon station (4298 m), although simulation results underestimate NH + 4 concentrations. The highly variable observed OM depth profile that results from episodic pulses of OM deposition through the canyon could be equally well fit by several different a-ν combinations, none of these combinations produces a better fit 25 of the NH + 4 profile. Hülse et al. (2018) invoked a number of factors that compromise the model-data fit of the NH + 4 profile at this site. First, Epping et al. (2002), who also did not obtain a better fit with the diagenetic model OMEXDIA, suggest intense bioirrigation activity as a potential driver of higher NH + 4 concentrations. Due to mathematical constrains, OMEN-SED merely includes a simplified description of bioirrigation that would require site specific fitting to account for locally enhanced irrigation rates. In addition, Hülse et al. (2018) showed that the model data fit can be improved by choosing a different NC-30 ratio and, thus, reflecting observed variations in Redfield stoichiometries at the Iberian margin and at the Nazaré Canyon sites (Epping et al., 2002).
Overall, a comparison between OMEN-SED and OMEN-SED-RCM results shows that both models reproduce observations equally well. They mainly differ in simulating the deeper parts of the NH + 4 profiles. These differences can be directly linked to the different descriptions of the OM degradation dynamics. While the 2G-model formulation (i.e. OMEN-SED) results in an abrupt decrease in OM degradation rates once the more reactive fraction is consumed, however the RCM (i.e. OMEN-SED-RCM) accounts for the continuous decrease in OM degradation, thus resulting in a continuous consumption of SO 2− 4 and 5 production of NH − 4 in the first layers under the SWI. Finally, inversely fitted OM reactivity parameters a and ν reveal a dominance of highly reactive to reactive OM compounds in the bulk OM across all sites. The consumption of these reactive compounds in the upper sediment layers results in a rapid decrease in apparent OM reactivity with depth for all sites (see Tables C1). These results are not only fully consistent with the previously determined 2G-Model parameters (Hülse et al., 2018), but also explain the good performance of the 2G Model and 10 the small differences between the RCM and 2G Model for these sites. The observed initially high, but rapidly decreasing OM reactivity can be, in contrast to a more slowly decreasing OM reactivity, adequately reproduced by a 2G-Model description.   following an empirical relationship proposed by Boudreau (1997). They assume that this rate constant is broadly representative  (Boudreau and Ruddick, 1991;Arndt et al., 2013;Freitas et al., 2020). Based on these previous global results, we 20 thus assume that parameter a controls the global variability in OM reactivity and that parameter ν remains roughly constant at ν = 0.15 across different depositional environments. Because parameter a is conceptually related to the average lifetime of the more reactive OM components and, thus, the degree of degradation or freshness of OM, the residence time of OM in the water column should exert, among other factors, a control on its magnitude (Boudreau and Ruddick, 1991;Middelburg, 1989). Boudreau and Ruddick (1991) and (Arndt et al., 2013) found, indeed, a weak trend between parameter a and sedimentation 25 rate: Log 10 (a) = 3.35 − 14.81 · ω, r 2 = 0.46.
Here, we use the global relationship to estimate parameter a as a function of ω for sites deeper than 200m water depths along the global hypsometry by using an empirical relationship between ω and SFD derived from Middelburg et al. (1997) (Table 2).
For high sedimentation rates and thus shallow water depths, the limited available dataset does not fully cover the variability 30 in OM reactivities that results from the dynamic mix of different OM sources in combination with the complex interplay of different controlling factors, such as physical protection mechanisms, microbial community composition or macrobenthic activity (Arndt et al., 2013). As a result eq. 15 tends to underestimate OM reactivity in the coastal ocean. Therefore, we here assume that the shallowest site (SFD<200m) has an apparent OM reactivity that is characteristic for fresh phytoplankton (Boudreau and Ruddick, 1991). As a result, apparent OM reactivities along the global ocean transect decreases over six orders of magnitude from the shallowest to the deepest SFD: k = 53.09 yr −1 (a = 2.8·10 −3 yr) at 100m to k = 7.45 10e −5 yr −1 (a = In addition to OM reactivity, the reoxidation of reduced substances (i.e. γ N H4 , γ H2S ) substantially affects the simulated benthic fluxes across the SWI. Here we assume that a fixed fraction of N H 4 and H 2 S is oxidised (i.e. γ N H4 = γ H2S = 0.95). We This difference can be directly attributed to the different representations of OM degradation dynamics in the three models (i.e. numerical 1-G, analytical 1-G and analytical 500-G model RCM approximation). The 1-G model applied by Thullner et al. 20 (2009) and Hülse et al. (2018) assumes a constant OM reactivity for the entire sediment column and, therefore, overestimates bulk OM reactivity in the deeper, often anoxic sediment layers, resulting in higher TEA uptake. The RCM approximation on the other hand, accounts for the entire reactivity continuum, and thus captures the decrease in apparent OM reactivity with sediment depth. As a result, OM degradation rates decrease faster with sediment depth, resulting in a reduced TEA uptake.  that are more representative of the Eastern Pacific Ocean (see Bohlen et al., 2012). Finally, OMEN-SED-RCM also captures the observed trend in SO 2− 4 uptake fluxes well. Simulation results also agree well with the OMEN-SED results that assume an almost complete re-oxidation of the deep H 2 S flux (i.e. γ H2S = 95% (black line with black circle, hidden by the OMEN-SED-RCM results). In contrast, BRNS model results, as well as OMEN-SED results that assume a weak re-oxidation of the deep ) and temperature at the SWI are constrained based on regionally averaged observational data extracted from the World Ocean Atlas (Boyer et al., 2009). The variables are reported with a 1°×1°horizontal resolution and we assume the last vertical cell at every grid-coordinate represents bottom water solute concentrations. The initial distribution of organic matter compounds across the reactivity spectrum at the SWI (i.e. parameters a and ν) and, therefore, the apparent reactivity of the depositing OM is constrained based on a compilation of inversely determined RCM parameters across a wide range of different depositional environments (Boudreau and Ruddick, 1991;Arndt et al., 2013;Freitas et al., 2020). Previous studies of inverse models showed that the global variability in apparent OM reactivity is mostly controlled by variations in parameter a, while parameter ν remains comparably constant across environments (Boudreau and Ruddick, 1991;Freitas et al., 2020). Building on these findings, we thus assume a globally constant parameter ν value of 0.15. In contrast, 15 parameter a changes globally. Conscious of the risk to over-parametrize the model, we still use the weak trend in parameter a and sedimentation rate found in Arndt et al. (2013), which follows roughly the observed broad global decrease in apparent OM reactivity from coastal to slope and abyssal depositional environments (Arndt et al., 2013). Additionally, the lack of estimates of parameter a in shallow environments makes this a-ω trend especially weak in such areas. 20 However, our global boundary condition on OM concentration is limited to depth > 1000m, hence remaining within the limits of the trend. We assume the a values implicitly accounts for broad patterns in benthic environmental characteristics, such as, among others: OM sources, transport pathways and transit times, oxygen exposure, mineral protection, microbial community structure etc. This approach thus captures the broad, generally observed reactivity trend across broadly defined depositional environments without over-parametrizing the problem. Finally, model transport and all additional reaction rate parameters are 25 chosen according to the original OMEN-SED set-up and summarized in Tables B1 and C1. global DOU observations and their global extrapolation (see Seiter et al., 2005). In addition, the simulated OPDs generally agree well with observations and capture also the widely observed increase in OPD with water depths (see Fig. 6). The highest 30 O 2 uptake rates ( 38 mmol m −2 yr −1 ) and the shallowest OPD are simulated for sediments underlying the western continental margins, the Arabian Sea, parts of the Greenland Sea and Laptev Sea and the equatorial Pacific (Fig. 5 a).

Results
These regions are characterized by the highest apparent OM reactivity, as well as high OM deposition fluxes (Seiter et al., 2005). In these regions, high apparent OM reactivity has been linked to either the efficient transport of large amounts of fresh, marine OM in the form of seasonal pulses or the presence of extended oxygen minimum zones. Both environmental drivers reduce the degree of pelagic OM degradation. Ecosystems that are characterized by a strong seasonality often deliver OM in strong post-bloom seasonal pulses to the benthic environment (Cavan et al., 2017;Cowie, 2005;Vandewiele et al., 2009;5 Cowie et al., 1999;Suthhof et al., 2001;Keil et al., 2016). These pulses favour the formation of fast sinking phytodetrital aggregates that reduce the residence time of sinking OM in the water column and, thus, the degree pelagic degradation of OM.
Pronounced oxygen minimum zones on the other hand reduce the exposure of OM to oxygen (Fischer et al., 2009;Zonneveld et al., 2010;Vandewiele et al., 2009). In both types of oceanic region, the resulting strong coupling between the euphotic zone and the sediment supports high oxygen uptake rates and shallow OPD. In agreement with global observations, OMEN-SED- 10 RCM results reveal that, in particular, the eastern part of the equatorial Pacific reveals high oxygen fluxes into the sediment (>4 mmol m −2 yr −1 ). However, oxygen uptake decreases westwards, as well as north and southwards due to the rapid decrease in OM flux away from the equatorial upwelling area (Seiter et al., 2005;Smith et al., 1997). In the western equatorial Pacific, higher DOU rates are supported by the increased concentrations of OM close to the active margins in southeast Asia, where high riverine loads deliver OM to marine sediments. High oxygen uptake and shallow OPD are also observed in sediments 15 underlying eastern boundary currents. However, here, OM reactivities are slightly lower due to the increased input of terrestrial material and strong lateral OM fluxes that redistribute pre-aged OM across the wide shelf of the passive Atlantic margin (Alonso-González et al., 2009;Mollenhauer et al., 2007Mollenhauer et al., , 2003Arthur et al., 1998;Inthorn et al., 2006;Lovecchio et al., 2018) as well as on the narrow shelf of the active margins of the Pacific Ocean (Muñoz et al., 2004;Venkatesan and Kaplan, 1992).
In agreement with global observations, generally low O 2 fluxes and deeper OPD are simulated for deep open ocean sed-20 iments, like the north (DOU rates < 1 mmol m −2 yr −1 ) and south Pacific (DOU rates < 2 mmol m −2 yr −1 ) and the Indian Ocean (DOU rates = [0.5 -3] mmol m −2 yr −1 ). Here, both the low magnitude and low quality of the depositing OM limits organic matter degradation and oxygen penetrates deep into the sediment (Glud, 2008;Smith et al., 2001). The lowest O 2 uptake rates are generally observed and also simulated for the central, oligotrophic gyre regions of the ocean basins that receive little OM input. Higher O 2 uptake can be observed close to the continental slopes, where lateral OM fluxes support a higher OM 25 deposition (Seiter et al., 2005;Archer and Devol, 1992). Atypically high O 2 uptake rates are observed as well as simulated for areas in the deep south Indian Ocean and south Atlantic Ocean. In these areas, seasonal pulses of primary production in the Antarctic polar front support an enhanced deposition of OM to the sediment and, thus support higher benthic oxygen uptake rates (Rabouille et al., 1998).
Furthermore, the reoxidation of reduced products can make up a substantial fraction of the oxygen consumption in marine 30 sediments. These secondary reactions consume products originating from anaerobic mineralization [Aller, 1990;Boudreau and Canfield, 1993] and can make up to 56% (Jørgensen and Kasten, 2006) of the total oxygen consumption in certain environments; thus cannot be ignored. We take these processes into account by running the full version of OMEN-SED-RCM.
Additionally, most of the global simulation here occurs in deeper parts of the ocean, where oxygen is the main oxidant of OM. Finally, OMEN-SED-RCM simulation results also reflect the widely observed inverse relationship between OPD and oxygen uptake rates (Rullkötter, 2006;Glud, 2008;Wenzhöfer and Glud, 2002). Shallow OPDs of less than 10 cm are observed for O 2 uptake rates of ≥ 4 mmol m −2 yr −1 . In contrast, the deep ocean and the central ocean gyres reveal OPDs of several centimetres to metres. Such deep OPD agree with observations from the central South Pacific. Here, D'Hondt et al. (2015) reported OPD that reach down to the basement at 120 m sediment depth. These deep OPD are directly linked to the extremely low OM 5 deposition, as well as apparent OM reactivities (Røy et al., 2012) that are reflected in low, observed DOU rates (D'Hondt et al., 2009(D'Hondt et al., , 2015Røy et al., 2012;Fischer et al., 2009).   et al. (1989); Canfield et al. (1993); Glud et al. (1994Glud et al. ( , 1998Glud et al. ( , 1999Glud et al. ( , 2003Glud et al. ( , 2009

Methodology
Sulphate is the most important terminal electron acceptor in marine sediments on a global scale (Jørgensen, 1982;Thullner et al., 2009). Its depth profile is a macroscopic manifestation of long-term OM degradation in the deeper anoxic sediment and, thus, crucial for understanding benthic biogeochemical dynamics. Sulfate depth profiles typically show a down core 5 decrease where the gradient is primarily controlled by OM degradation rates and thus both the amount and quality of OM in the deeper sediment. The sediment depth at which downward diffusing sulphate is completely consumed is referred to the sulphate-methane transition zone (SMTZ). It separates the sulphate reduction from the methanogenic zone. The SMTZ is a very dynamic redox boundary and its position is very sensitive to changes in environmental conditions (Regnier et al., 2011;Meister et al., 2013). Here, we use OMEN-SED-RCM to simulate changes in the position of the SMTZ in response to a range 10 of sedimentation rates and OM reactivity parameters and compare the simulated patterns with previously observed/simulated patterns (Regnier et al., 2011;Meister et al., 2013).
To this end we set-up a one-dimensional OMEN-SED-RCM simulation that is representative for a shallow depositional environment with a high OM input and a high sedimentation rate. OM concentration is set to 1.1 wt%. We assume bottom water concentrations of SO 2− 4 = 28 mM and CH 4 = 0 mM at the SWI. The model domain is set to 45 m and thus accounts 15 for the active part of the anoxic sediment layer. We use this OMEN-SED-RCM set-up to explore the response of the SMTZ to variations in sedimentation rates and apparent OM reactivity (i.e. parameter a). An overview of the applied boundary conditions and parameters is given in Table 3. Because, we do not consider active methane venting or upward fluid flow, the location of the SMTZ is mainly controlled by OM degradation.

Results
In agreement with results of Regnier et al. (2011) and Meister et al. (2013), OMEN-SED-RCM simulation results show that the sedimentation rate, ω, and parameter a exert a dominant control on the location of the SMTZ (Fig. 7). OMEN-SED-RCM captures the observed/simulated response of the SMTZ to changing environmental conditions well and results are similar to those found in Regnier et al. (2011) and Meister et al. (2013)  Simulation results also show that OM reactivity exerts the dominant control on the depth of the SMTZ. Changes in sedimentation rate and OM deposition shift (not shown here), but do not change the general pattern. Increasing the deposition flux of OM at the SWI would lead merely to the shallowing of the SMTZ on both ends of the parameter a spectrum (not shown here). 20 Increasing sedimentation rates generally result in shallower SMTZ in the area of high OM reactivities because they increase the flux of OM to the deeper sediment layers. For the same reason, high sedimentation rates also decrease the minimum depth of the SMTZ for intermediate reactivities. The shallowest SMTZ (-5 m) is simulated for the highest ω value and a = 10 2.7 years, while the deepest at -20m occurs when ω is lowest and a = 10 3.6 years. In the area of low reactivities (i.e. high parameter a), sedimentation rates, on the other hand, only exert a limited effect on the location of the SMTZ because higher sedimentation 25 rates merely result in a slight increase in OM reactivity at depth. require an exponential form of the OM degradation term, the RCM is approximated by a multi-G approximation, applied here 5 with 500 fractions (Boudreau and Ruddick, 1991;Dale et al., 2015). The RCM and its multi-G approximation (i.e. reaction rate constants and initial fraction of each OM class) is fully constrained by the two free RCM parameters a and ν that control the initial distribution of OM compounds over the reactivity spectrum.
We show that the new version of OMEN-SED is not only able to reproduce observed pore water profiles across a wide range of depositional environments and captures observed global patterns of TEA-fluxes, oxygen penetration depths, biogeochemical 10 reaction rates, but also accounts for the widely observed continuous decrease in OM reactivity with sediment depths/burial time and, thus provides a more realistic description of anaerobic degradation pathways. This added functionality offers an alternative to the common, but simpler 2G-model description implemented in the original model, extending the model's applicability to a wider range of environments and timescales, while requiring fewer parameters to describe a wider spectrum of OM reactivity.
These improvements were implemented while maintaining the computational advantages of the original version.