OMEN-SED 1.0: a novel, numerically efficient organic matter sediment diagenesis module for coupling to Earth system models

We present the first version of OMEN-SED (Organic Matter ENabled SEDiment model), a new, one-dimensional analytical early diagenetic model resolving organic matter cycling and the associated biogeochemical dynamics in marine sediments designed to be coupled to Earth system models. OMEN-SED explicitly describes organic matter (OM) cycling and the associated dynamics of the most important terminal electron acceptors (i.e. O₂ , NO₃, SO₄) and methane (CH₄), related reduced substances (NH₄, H₂S), macronutrients (PO₄) and associated pore water quantities (ALK, DIC). Its reaction network accounts for the most important primary and secondary redox reactions, equilibrium reactions, mineral dissolution and precipitation, as well as adsorption and desorption processes associated with OM dynamics that affect the dissolved and solid species explicitly resolved in the model. To represent a redox-dependent sedimentary P cycle we also include a representation of the formation and burial of Fe-bound P and authigenic Ca–P minerals. Thus, OMEN-SED is able to capture the main features of diagenetic dynamics in marine sediments and therefore offers similar predictive abilities as a complex, numerical diagenetic model. Yet, its computational efficiency allows for its coupling to global Earth system models and therefore the investigation of coupled global biogeochemical dynamics over a wide range of climate-relevant timescales. This paper provides a detailed description of the new sediment model, an extensive sensitivity analysis and an evaluation of OMEN-SED's performance through comprehensive comparisons with observations and results from a more complex numerical model. We find that solid-phase and dissolved pore water profiles for different ocean depths are reproduced with good accuracy and simulated terminal electron acceptor fluxes fall well within the range of globally observed fluxes. Finally, we illustrate its application in an Earth system model framework by coupling OMEN-SED to the Earth system model cGENIE and tune the OM degradation rate constants to optimise the fit of simulated benthic OM contents to global observations. We find that the simulated sediment characteristics of the coupled model framework, such as OM degradation rates, oxygen penetration depths and sediment–water interface fluxes, are generally in good agreement with observations and in line with what one would expect on a global scale. Coupled to an Earth system model, OMEN-SED is thus a powerful tool that will not only help elucidate the role of benthic–pelagic exchange processes in the evolution and the termination of a wide range of climate events, but will also allow for a direct comparison of model output with the sedimentary record – the most important climate archive on Earth.

1 Introduction

Marine surface sediments are key components in the Earth system. They host the largest carbon reservoir within the surficial Earth system, provide the primary long-term sink for atmospheric CO₂, recycle nutrients and represent the most important geochemical archive used for deciphering past changes in biogeochemical cycles and climate (e.g. Berner, 1991; Archer and Maier-Reimer, 1994; Ridgwell and Zeebe, 2005; Arndt et al., 2013). Physical and chemical processes in sediments (i.e. diagenetic processes) depend on the water column and vice versa: diagenesis is controlled by the external supply of solid material (e.g. organic matter, calcium carbonate, opal) from the water column and is affected by overlying bottom water concentrations of solutes. At the same time, sediments impact the water column directly either through the short- and long-term storage of deposited material or the diagenetic processing of deposited material and the transport of terminal electron acceptors (e.g. O₂, SO₄) into the sediments, as well as metabolic products (e.g. nutrients, DIC) to the overlying bottom waters. This so-called benthic–pelagic coupling is essential for understanding global biogeochemical cycles and climate (e.g. Archer and Maier-Reimer, 1994; Archer et al., 2000; Soetaert et al., 2000; Mackenzie, 2005).

The biological primary production of organic matter (OM, generally represented in its simple form CH₂O in Reaction R1) and the reverse process of degradation can be written in a greatly simplified reaction as

$$\begin{array}{}\text{(R1)}& {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}.\end{array}$$

On geological timescales, production of OM is generally greater than degradation, which results in some organic matter being buried in marine sediments and oxygen accumulating in the atmosphere. Thus, the burial of OM deep into the sediment leads to net oxygen input to and CO₂ removal from the atmosphere (Berner, 2004). On shorter timescales, the upper few metres of the sediments where early diagenesis occurs are specifically important, as this zone controls whether a substance is recycled to the water column or buried for a longer period of time in the deeper sediments (Hensen et al., 2006). Most biogeochemical cycles and reactions in this part of marine sediments can be related either directly or indirectly to the degradation of organic matter (Middelburg et al., 1993; Arndt et al., 2013). Oxygen and nitrate, for instance, the highest energy-yielding electron acceptors, are preferentially consumed in the course of the degradation of organic matter, resulting in the release of ammonium and phosphorus to the pore water. As such, the degradation of OM in the sediments can profoundly affect the oxygen and nutrient inventory of the ocean and thus primary productivity (Van Cappellen and Ingall, 1994; Lenton and Watson, 2000). Furthermore, organic matter degradation releases metabolic CO₂ to the pore water, causing it to have a lower pH and carbonate ion concentration, thus provoking the dissolution of calcium carbonate CaCO₃ (Emerson and Bender, 1981).

Benthic nutrient recycling from marine sediments has been suggested to play a key role for climate and ocean biogeochemistry throughout Earth history. For example, feedbacks between phosphorus storage and erosion from shelf sediments and marine productivity have been hypothesised to play an important role for glacial–interglacial atmospheric CO₂ changes (Broecker, 1982; Ruttenberg, 1993). Furthermore, benthic nutrient recycling from anoxic sediments has been invoked to explain the occurrence of more extreme events in Earth history, for instance oceanic anoxic events (OAEs; e.g. Van Cappellen and Ingall, 1994; Mort et al., 2007; Tsandev and Slomp, 2009). OAEs represent severe disturbances of the global carbon, oxygen and nutrient cycles of the ocean and are usually characterised by widespread bottom water anoxia and photic zone euxinia (Jenkyns, 2010). One way to explain the genesis and persistence of OAEs is increased oxygen demand due to enhanced primary productivity. Increased nutrient inputs to fuel primary productivity may in turn have come from marine sediments as the burial efficiency of phosphorus declines when bottom waters become anoxic (Ingall and Jahnke, 1994; Van Cappellen and Ingall, 1994). The recovery from OAE-like conditions is thought to involve the permanent removal of excess CO₂ from the atmosphere and ocean by burying carbon in the form of organic matter in marine sediments (e.g. Arthur et al., 1988; Jarvis et al., 2011), which is consistent with the geological record of widespread black shale formation (Stein et al., 1986). Models capable of simulating not only the expansion and intensification of oxygen minimum zones, but also of predicting how the underlying sediments interact are hence needed.

Quantifications of diagenetic processes in the sediments are possible through the application of idealised mathematical representations, or so-called diagenetic models (see e.g. Berner, 1980; Boudreau, 1997). A plethora of different approaches have been developed, mainly following two distinct directions (see Arndt et al., 2013, for an overview). The first involves state-of-the-art vertically resolved numerical models simulating the entire suite of essential coupled redox and equilibrium reactions within marine sediments (e.g. BRNS, Aguilera et al., 2005; CANDI, Boudreau, 1996; MEDIA, Meysman et al., 2003; MUDS, Archer et al., 2002; STEADYSED, Van Cappellen and Wang, 1996). These “complete”, multi-component steady-state or non-steady-state models thus resolve the resulting characteristic redox zonation of marine sediments by explicitly accounting for oxic OM degradation, denitrification, oxidation by manganese and iron (hydr)oxides, sulfate reduction and methanogenesis as well as the reoxidation of reduced byproducts (i.e. NH₄, Mn²⁺, Fe²⁺, H₂S, CH₄; see e.g. Regnier et al., 2011). Furthermore, they incorporate various mineral dissolution and precipitation reactions, as well as fast equilibrium sorption processes, for example of NH₄, PO₄ and metal ions (i.e. Mn²⁺, Fe²⁺ and Mg²⁺; compare Van Cappellen and Wang, 1996; Meysman et al., 2003). Modelled, depth-dependent transport processes usually comprise advection, diffusion, bioturbation and bio-irrigation. This group of diagenetic models generally describes OM degradation via a so-called multi-G approach (Jørgensen, 1978; Berner, 1980), thus dividing the bulk organic matter pool into a number of compound classes that are characterised by different degradabilities k_i. Alternative approaches, so-called continuum models (Middelburg, 1989; Boudreau and Ruddick, 1991), assume a continuous distribution of reactive types but, although conceptually superior, are much less popular (Arndt et al., 2013). These complex, multi-component models have great potential for quantifying diagenetic dynamics at sites where comprehensive observational datasets are available to constrain its model parameters (see e.g. Boudreau et al., 1998; Wang and Van Cappellen, 1996; Thullner et al., 2009, for applications). However, due to the high degree of coupled processes and depth-varying parameters, the diagenetic equation needs to be solved numerically, thus resulting in a very high computational demand and consequently rendering their application in an Earth system model (ESM) framework with a large number of grid points prohibitive. Additionally, their global applicability is seriously compromised by the restricted transferability of model parameters from one site to the global scale (Arndt et al., 2013).

The second group of diagenetic models emerged during the early days of diagenetic modelling when computing power was severely restricted (e.g. Berner, 1964). These models solve the diagenetic equation analytically, thus providing an alternative and computationally more efficient approach. Finding an analytical solution, especially when complex reaction networks are to be considered, is not straightforward and analytical models are thus usually less sophisticated and comprehensive than numerical models and generally require the assumption of steady-state conditions. It has been shown that the complexity of the reaction network can be reduced by dividing the sediment column into distinct zones and accounting for the most pertinent biogeochemical processes within each zone, thus increasing the likelihood of finding an analytical solution without oversimplifying the problem. Analytical approaches with distinct biogeochemical zones were implemented and used in the 1970s and 1980s to describe observed pore water profiles (e.g. Vanderborght and Billen, 1975; Vanderborght et al., 1977; Billen, 1982; Goloway and Bender, 1982; Boudreau and Westrich, 1984) and later for inclusion into multi-box ecosystem models (e.g. Ruardij and Van Raaphorst, 1995; Gypens et al., 2008) and global Earth system models (Tromp et al., 1995). However, in addition to the oxic zone these models generally only describe one anoxic zone explicitly, either a denitrification (Vanderborght and Billen, 1975; Billen, 1982; Goloway and Bender, 1982; Ruardij and Van Raaphorst, 1995; Gypens et al., 2008) or a sulfate reduction zone (Boudreau and Westrich, 1984; Tromp et al., 1995). Furthermore, the approaches of Vanderborght and Billen (1975), Goloway and Bender (1982) and Tromp et al. (1995) do not explicitly account for reduced species (i.e. NH₄ and H₂S).

In most current ESMs sediment–water dynamics are either neglected or treated in a very simplistic way (Soetaert et al., 2000; Hülse et al., 2017). Most Earth system models of intermediate complexity (EMICs) and also some of the higher-resolution Earth system–climate models represent the sediment–water interface either as a reflective or a conservative–semi-reflective boundary (Hülse et al., 2017). Thus, all particulate material deposited on the sea floor is either instantaneously consumed (reflective boundary), or a fixed fraction is buried in the sediments (conservative–semi-reflective boundary). Both highly simplified approaches furthermore completely neglect the exchange of solute species through the sediment–water interface and therefore cannot resolve the complex benthic–pelagic coupling. However, due to their computational efficiency, both representations are often used in global biogeochemical models (e.g. Najjar et al., 2007; Ridgwell et al., 2007; Goosse et al., 2010). Analytical diagenetic models represent the most complex description of diagenetic dynamics in Earth system models. Examples of global ESMs employing a vertically resolved diagenetic model are NorESM (Tjiputra et al., 2013) and HAMOCC (Palastanga et al., 2011; Ilyina et al., 2013), both using a version of Heinze et al. (1999). None of the EMICs reviewed by Hülse et al. (2017) use such a sediment representation. DCESS (Shaffer et al., 2008) and MBM (Munhoven, 2007) are box models employing a vertically resolved diagenetic model. These analytic models account for the most important transport processes (i.e. advection, bioturbation and molecular diffusion) through basic parameterisations and include fewer biogeochemical reactions, which are generally restricted to the upper, bioturbated 10 cm of the sediments. Pore water species explicitly represented in DCESS (Shaffer et al., 2008) and the HAMOCC model of Heinze et al. (1999) and Palastanga et al. (2011) are restricted to DIC, TA, PO₄ and O₂. The MEDUSA model (Munhoven, 2007) considers CO₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$, ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and O₂. Other species produced or consumed during OM degradation are neglected. Thus, with oxygen being the only TEA explicitly modelled, the influence of reduced species is only implicitly included in the boundary conditions for O₂. A newer version of the HAMOCC model is a notable exception, as Ilyina et al. (2013) include NO₃ and denitrification explicitly. Furthermore, the version of Palastanga et al. (2011) represents a redox-dependent explicit sedimentary phosphorus cycle. Yet, the reoxidation of reduced byproducts, so-called secondary redox reactions (e.g. oxidation of NH₄, H₂S or CH₄), or sorption processes are not included in any of the discussed models. Furthermore, these global models assume that the sedimentary organic matter pool is composed of just a single compound class which is either degraded with a globally invariant degradation rate constant (Munhoven, 2007) or a fixed rate constant depending on local oxygen concentrations (Shaffer et al., 2008; Palastanga et al., 2011).

Obviously, such a simplification of the OM pool can neither account for the observed vast structural complexity in natural organic matter and its resulting different degradation rates nor for the rapid decrease in OM degradability in the uppermost centimetres of the sediments (Arndt et al., 2013). It has been suggested that at least a 3G approach is necessary to accurately represent organic matter dynamics in this part of the sediments where most OM is degraded (e.g. Soetaert et al., 1996). Even more restrictive is the use of O₂ as the only TEA and the complete absence of reduced substances and related secondary redox reactions. For the majority of the modern sediments (i.e. in the deep ocean) O₂ is the primary electron acceptor; however, recent model and data studies have reported that sulfate reduction is the dominant degradation pathway on a global average (with contributions of 55–76 % Canfield et al., 2005; Jørgensen and Kasten, 2006; Thullner et al., 2009). Oxygen becomes progressively less important as TEA with decreasing sea-floor depth and sulfate reduction has been shown to account for 83 % of OM degradation in coastal sediments (Krumins et al., 2013). In these environments most O₂ is used to reoxidise reduced substances produced during anaerobic degradation (Canfield et al., 2005; Thullner et al., 2009). Thus, the in situ production of e.g. NO₃ and SO₄ through the oxidation of NH₄ and H₂S forms an important sink for O₂ which is entirely neglected in current sediment representations in global models. In addition, the lack of anoxic degradation pathways in these models limits their application to oxic oceans. Currently no analytical sediment model exists that can be used under anoxic conditions. Due to the lack of an appropriate sedimentary P cycle (with the exception of the HAMOCC version of Palastanga et al., 2011), no current global ESM is able to model the redox-dependent P release from marine sediments and its implications for primary productivity, global biogeochemical cycles and climate. A sediment model suitable for coupling to an ESM and enabling a wide range of paleo-questions to be addressed has to provide a robust quantification of organic (and inorganic) carbon burial fluxes, benthic uptake and return fluxes of oxygen, growth-limiting nutrients and reduced species, as well as anoxic degradation pathways. As a consequence, the reaction network must account for the most important primary and secondary redox reactions, equilibrium reactions, mineral precipitation and dissolution, and adsorption and desorption, resulting in a complex set of coupled reaction–transport equations.

Therefore, we developed the Organic Matter ENabled SEDiment model (OMEN-SED), a new, one-dimensional, numerically efficient diagenetic model. OMEN-SED builds upon and stands in the tradition of earlier stand-alone, analytical diagenetic models (Vanderborght et al., 1977; Billen, 1982; Goloway and Bender, 1982; Boudreau, 1991) and analytical diagenetic models developed for coupling to regional-scale ecosystem or global Earth system models (Ruardij and Van Raaphorst, 1995; Tromp et al., 1995; Heinze et al., 1999; Gypens et al., 2008).

OMEN-SED is the first analytical model to explicitly describe OM cycling and the associated dynamics of the most important TEAs (i.e. O₂, NO₃, SO₄), related reduced substances (NH₄, H₂S), the full suite of secondary redox reactions, macronutrients (PO₄) and the associated pore water quantities (ALK, DIC). To represent a redox-dependent sedimentary P cycle we consider the formation and burial of Fe-bound P and authigenic Ca–P minerals. Thus, while OMEN-SED captures most of the features of a complex, numerical diagenetic model, its computational efficiency allows for coupling to global Earth system models and therefore the investigation of coupled global biogeochemical dynamics over different timescales. Here, the model is presented as a 2G approach; however, OMEN-SED can be easily extended to a multi-G approach. The first part of the paper provides a detailed description of OMEN-SED (Sect. 2). This includes descriptions of the general model approach (Sect. 2.1), the conservation equations for all explicitly represented biogeochemical tracers (Sect. 2.2) and a summary of global relationships used to constrain reaction and transport parameters in OMEN-SED (Sect. 2.4). In addition, a generic algorithm is described which is used to match internal boundary conditions and to determine the integration constants for the analytical solutions (Sect. 2.3). In order to validate the stand-alone version of OMEN-SED, the second part of the paper performs an extensive sensitivity analysis for the most important model parameters, and resulting sediment–water interface fluxes are compared with a global database (Sect. 3.1). In addition, the results of the stand-alone model are compared with observed pore water profiles from different ocean depths (Sect. 3.2), and OMEN-SED simulations of TEA fluxes along a typical ocean transect are compared with observations and results from a complete, numerical diagenetic model (Sect. 3.3). Thereafter, OMEN-SED is coupled to the carbon-centric version of the GENIE Earth system model (cGENIE; Ridgwell et al., 2007, Sect. 4.1). Sensitivity studies are carried out using this coupled model and modelled organic matter concentrations in the surface sediments are compared to a global database (Seiter et al., 2004, Sect. 4.2). We finally discuss potential applicabilities of OMEN-SED and critically analyse model limitations (Sect. 5).

2 Model description

OMEN-SED is implemented as a FORTRAN version that can be easily coupled to any pelagic, biogeochemical model via the coupling routine OMEN_SED_main. In addition, OMEN-SED exists as a stand-alone version implemented in MATLAB and the entire model can be executed on a standard personal computer in less than 0.1 s. The source code of both the FORTRAN and the MATLAB stand-alone version and instructions for executing OMEN-SED and for plotting model results are available as a Supplement to this paper.

The following section provides a detailed description of OMEN-SED and the fundamental equations underlying the model are highlighted. Tables 1 and A1 summarise the biogeochemical reaction network and Tables 9 and 10 provide a glossary of model parameters along with their respective units.

2.1 General model approach

In OMEN-SED, the calculation of benthic uptake, recycling and burial fluxes is based on the vertically resolved conservation equation for solid and dissolved species in porous media (e.g. Berner, 1980; Boudreau, 1997):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\xi }{\displaystyle \frac{\partial C_i}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left(-\mathit{\xi }{D}_i{\displaystyle \frac{\partial C_i}{\partial z}}+\mathit{\xi }w{C}_i\right)+\mathit{\xi }\sum _j{R}_i^j,\end{array}$$

where C_i is the concentration of biogeochemical species i, and ξ equals the porosity ϕ for solute species and (1−ϕ) for solid species. The term z is the sediment depth, t denotes the time, D_i is the apparent diffusion coefficient of species i, w is the advection rate and ${\sum }_j{R}_i^j$ represents the sum of all biogeochemical rates j affecting species i.

Figure 1 Schematic of the different modelled species and zones in OMEN-SED. Here the case $z_{\mathrm{ox}}<z_{\mathrm{bio}}<z_{{\mathrm{NO}}_{\mathrm{3}}}<z_{{\mathrm{SO}}_{\mathrm{4}}}$ is shown.

OMEN-SED accounts for both the advective and the diffusive transport of solid and dissolved species. They are buried in the sediment according to a constant advection rate w, thus neglecting the effect of sediment compaction (i.e. $\frac{\partial \mathit{\varphi }}{\partial z}=\mathrm{0}$) due to mathematical constraints. The molecular diffusion of dissolved species is described by Fick's law applying a species-specific apparent diffusion coefficient, D_mol,i. In addition, the activity of infaunal organisms in the bioturbated zone is simulated using a diffusive term (e.g. Boudreau, 1986), with a constant bioturbation coefficient D_bio in the bioturbated zone, while D_bio is set to zero below the maximum bioturbation depth, z_bio. The pumping activity by burrow-dwelling animals and the resulting ventilation of tubes, the so-called bio-irrigation, is encapsulated in a factor f_ir that enhances the molecular diffusion coefficient (hence, $D_{i,\mathrm{0}}=D_{\mathrm{mol},i}\cdot f_{\mathrm{ir}}$; Soetaert et al., 1996). The reaction network of OMEN-SED accounts for the most important primary and secondary redox reactions, equilibrium reactions, mineral dissolution and precipitation, and the adsorption and desorption processes associated with OM dynamics that affect the dissolved and solid species explicitly resolved in the model. Tables 1 and A1 provide a summary of the reactions and biogeochemical tracers considered in OMEN-SED together with their respective reaction stoichiometries.

Table 1 Reactions and biogeochemical tracers implemented in the reaction network of OMEN-SED. The primary and secondary redox reactions are listed in the sequence they occur with increasing sediment depth.

All parameters in Eq. (1), apart from porosity and burial rate, may vary with sediment depth and many reaction rate expressions depend on the concentration of other species. Expressing Eq. (1) for a set of chemical species thus results in a non-linear, coupled set of equations that can only be solved numerically. However, OMEN-SED is designed for coupling to Earth system models and therefore cannot afford a computationally expensive numerical solution. Instead, similar to early analytical diagenetic models, a computationally efficient analytical solution to Eq. (1) can be derived by (1) assuming steady-state conditions (i.e.$\frac{\partial C_i}{\partial t}=\mathrm{0}$) and (2) reducing the vertical variability in parameters and reaction rate expressions by dividing the sediment column into a number of functional biogeochemical zones (Fig. 1; compare e.g.  Billen, 1982; Goloway and Bender, 1982; Ruardij and Van Raaphorst, 1995; Tromp et al., 1995; Gypens et al., 2008, for similar solutions). More specifically, OMEN-SED follows Berner (1980) by dividing the sediment column into (I) a bioturbated and (II) a non-bioturbated zone defined by an imposed, constant bioturbation depth z_bio (Fig. 1). Furthermore, it resolves the dynamic redox stratification of marine sediments by dividing the sediment into (1) an oxic zone delineated by the oxygen penetration depth z_ox; (2) a denitrification (or nitrogenous) zone situated between z_ox and the nitrate penetration depth $z_{{\mathrm{NO}}_{\mathrm{3}}}$; (3) a sulfate reduction zone situated between $z_{{\mathrm{NO}}_{\mathrm{3}}}$ and the sulfate penetration depth $z_{{\mathrm{SO}}_{\mathrm{4}}}$; and (4) a methanogenic zone situated below $z_{{\mathrm{SO}}_{\mathrm{4}}}$ (Fig. 1). Although in each of these zones Eq. (1) is applied with depth invariant parameters, parameter values may differ across zones. The biogeochemical zones are linked by stating continuity in both concentrations and fluxes at the dynamic, internal boundaries ($z_{\mathrm{b}}\in \mathit{\{}{z}_{\mathrm{bio}},z_{\mathrm{ox}},z_{{\mathrm{NO}}_{\mathrm{3}}},z_{{\mathrm{SO}}_{\mathrm{4}}}\mathit{\}}$; compare e.g. Billen, 1982; Ruardij and Van Raaphorst, 1995). Note that these boundaries are dynamic because their depth varies in response to changing ocean boundary conditions and forcings (see Sect. 2.3.1 for details). Furthermore, the maximum bioturbation depth is not restricted to a specific biogeochemical zone, and hence OMEN-SED allows bioturbation to occur in the anoxic zones of the sediment (here all zones z>z_ox combined).

The formulation of the reaction term in Eq. (1) varies between zones and encapsulates the most pertinent reaction processes within the respective zone (see Sect. 2.2), thus simplifying the mathematical description of the reaction network while retaining most of its biogeochemical complexity. One such simplification is that solid-phase iron and manganese oxidants and their reductants are not considered in the reaction network. All consumption and production processes of dissolved species related to the degradation of organic matter are a function of the organic matter concentration. Because organic matter degradation is described as a first-order degradation, these processes can be expressed as a series of exponential terms (${\sum }_j{\mathit{\alpha }}_j\exp (-{\mathit{\beta }}_{j}z)$; see Eq. 2). In addition, slow adsorption and desorption and mineral precipitation processes can be expressed as zero- or first-order (reversible) reactions (Q_m or k_l⋅C_i in Eq. 2). Fast adsorption is described as an instantaneous equilibrium reaction using a constant adsorption coefficient K_i. The reoxidation of reduced substances is accounted for implicitly by adding a (consumption and production) flux to the internal boundary conditions (see Sect. 2.2.2, 2.2.3 and 2.2.4). This simplification has been used previously by Gypens et al. (2008) for nitrate and ammonium and can be justified as it has been shown that reoxidation mainly occurs within a thin layer at the oxic–anoxic interface (Soetaert et al., 1996). The general reaction–transport equation underlying OMEN-SED is thus given by

$$\begin{array}{ll}{\displaystyle \frac{\partial C_i}{\partial t}}& {\displaystyle }=\mathrm{0}={\displaystyle \frac{D_i}{\mathrm{1}+K_i}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{C}_i}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial C_i}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+K_i}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }\left(\sum _j{\mathit{\alpha }}_j\exp (-{\mathit{\beta }}_{j}z)+\sum _l{k}_l\cdot C_i-\sum _m{Q}_m\right),\end{array}$$

where 1∕β_j can be interpreted as the length scale and α_j as the relative importance (or the magnitude at z=0) of reaction j (Boudreau, 1997); k_l represents generic first-order reaction rate constants and Q_m represents zero-order (or constant) reaction rates.

The analytical solution to Eq. (2) is of the general form

$$\begin{array}{ll}{\displaystyle }{C}_i\left(z\right)& {\displaystyle }=A\cdot \exp \left(az\right)+B\cdot \exp \left(bz\right)+\sum _j{\displaystyle \frac{{\mathit{\alpha }}_j}{D{\mathit{\beta }}_j^{\mathrm{2}}-w{\mathit{\beta }}_j-{\sum }_l{k}_l}}\\ \text{(3)}& {\displaystyle }& {\displaystyle }\cdot \exp (-{\mathit{\beta }}_{j}z)+{\displaystyle \frac{{\sum }_m{Q}_m}{{\sum }_l{k}_l}},\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }a={\displaystyle \frac{w-\sqrt{w^{\mathrm{2}}+\mathrm{4}\cdot D\cdot {\sum }_l{k}_l}}{\mathrm{2}\cdot D}},\\ \text{(4)}& {\displaystyle }& {\displaystyle }b={\displaystyle \frac{w+\sqrt{w^{\mathrm{2}}+\mathrm{4}\cdot D\cdot {\sum }_l{k}_l}}{\mathrm{2}\cdot D}}\end{array}$$

where A and B are integration constants that can be determined by applying a set of internal boundary conditions (see Sect. 2.3) and $D=\frac{D_i}{\mathrm{1}+K_i}$.

Based on Eq. (2) and its analytical solution Eq. (3), OMEN-SED returns the fraction of particulate organic carbon (POC) buried in the sediment, f_POC, and the benthic uptake and return fluxes $F_{C_i}$ of dissolved species C_i (in mol cm⁻² yr⁻¹) in response to the dynamic interplay of transport and reaction processes under changing boundary conditions and forcings:

$$\begin{array}{}\text{(5)}& {\displaystyle }{f}_{\mathrm{POC}}& {\displaystyle }={\displaystyle \frac{\mathrm{POC}\left(z_{\max }\right)}{\mathrm{POC}\left(\mathrm{0}\right)}},\text{(6)}& {\displaystyle }{F}_{C_i}& {\displaystyle }=\mathit{\varphi }\left(\mathrm{0}\right)\left(D_i{\displaystyle \frac{\partial C_i\left(z\right)}{\partial z}}{|}_{z=\mathrm{0}}-w\cdot C_i\left(\mathrm{0}\right)\right),\end{array}$$

where w is the deposition rate, D_i is the diffusion coefficient and POC(0), POC(z_max) and C_i(0) denote the concentration of POC and dissolved species i at the sediment–water interface (SWI) and at the lower sediment boundary, respectively.

2.2 Conservation equations and analytical solution

The following sections provide a detailed description of the conservation equations and analytical solutions for each chemical species that is resolved in this version of OMEN-SED.

2.2.1 Organic matter or particulate organic carbon (POC)

In marine sediments, organic matter (or in the following called particulate organic carbon, POC) is degraded by heterotrophic activity coupled to the sequential utilisation of terminal electron acceptors (TEAs) according to the free energy gain of the half-reaction (${\mathrm{O}}_{\mathrm{2}}>{\mathrm{NO}}_{\mathrm{3}}^{-}>{\mathrm{MnO}}_{\mathrm{2}}>\mathrm{Fe}(\mathrm{OH}{)}_{\mathrm{3}}>{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$; e.g. Stumm and Morgan, 2012). Once all TEAs are depleted, organic matter is degraded via methanogenesis. Here, organic matter degradation is described via a multi-G model approach (Jørgensen, 1978), dividing the bulk OM into a number i of discrete compound classes POC_i characterised by class-specific first-order degradation rate constants k_i. The conservation equation for organic matter dynamics is thus given by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial {\mathrm{POC}}_i}{\partial t}}=\mathrm{0}=D_{{\mathrm{POC}}_i}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{POC}}_i}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{POC}}_i}{\partial z}}-k_i\cdot {\mathrm{POC}}_i,\end{array}$$

with $D_{{\mathrm{POC}}_i}=D_{\mathrm{bio}}$ for z≤z_bio and $D_{{\mathrm{POC}}_i}=\mathrm{0}$ for z>z_bio. Integration of Eq. (7) yields the following general solutions for the bioturbated and non-bioturbated layers.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(I) Bioturbated zone}(z\le z_{\mathrm{bio}})\\ \text{(8)}& {\displaystyle }& {\displaystyle }{\mathrm{POC}}_i^{\mathrm{I}}\left(z\right)=A_{\mathrm{1}i}\cdot \exp \left(a_{\mathrm{1}i}z\right)+B_{\mathrm{1}i}\cdot \exp \left(b_{\mathrm{1}i}z\right){\displaystyle }& {\displaystyle }\text{(II) Non-bioturbated zone}(z_{\mathrm{bio}}<z)\\ \text{(9)}& {\displaystyle }& {\displaystyle }{\mathrm{POC}}_i^{\mathrm{II}}\left(z\right)=A_{\mathrm{2}i}\cdot \exp \left(a_{\mathrm{2}i}z\right)\end{array}$$

In the above equations,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_{\mathrm{1}i}={\displaystyle \frac{w-\sqrt{w^{\mathrm{2}}+\mathrm{4}\cdot D_{{\mathrm{POC}}_i}\cdot k_i}}{\mathrm{2}\cdot D_{{\mathrm{POC}}_i}}},\\ \text{(10)}& {\displaystyle }& {\displaystyle }{b}_{\mathrm{1}i}={\displaystyle \frac{w+\sqrt{w^{\mathrm{2}}+\mathrm{4}\cdot D_{{\mathrm{POC}}_i}\cdot k_i}}{\mathrm{2}\cdot D_{{\mathrm{POC}}_i}}},a_{\mathrm{2}i}=-{\displaystyle \frac{k_i}{w}}.\end{array}$$

Determining the integration constants ($A_{\mathrm{1},i},B_{\mathrm{1},i},A_{\mathrm{2},i}$) requires the definition of a set of boundary conditions (Table 2). For organic matter, OMEN-SED applies a known concentration at the sediment–water interface and assumes continuity across the bottom of the bioturbated zone, z_bio. When OMEN-SED is coupled to an ESM, the POC depositional flux from the coupled ocean model is converted to a concentration by solving the flux divergence Eq. (51). The integration constants ($A_{\mathrm{1},i},B_{\mathrm{1},i},A_{\mathrm{2},i}$) are thus given by the following.

$$\begin{array}{ll}\text{(11)}& {\displaystyle }{B}_{\mathrm{1}i}& {\displaystyle }\stackrel{\mathrm{BC}\mathrm{1})}{=}{\mathrm{POC}}_{\mathrm{0}i}-A_{\mathrm{1}i}{\displaystyle }{A}_{\mathrm{2}i}& {\displaystyle }\stackrel{\mathrm{BC}\mathrm{2})}{=}{\displaystyle \frac{A_{\mathrm{1}i}\cdot \exp \left(a_{\mathrm{1}i}{z}_{\mathrm{bio}}^{-}\right)+B_{\mathrm{1}i}\cdot \exp \left(b_{\mathrm{1}i}{z}_{\mathrm{bio}}^{-}\right)}{\exp \left(a_{\mathrm{2}i}{z}_{\mathrm{bio}}^{+}\right)}}\\ {\displaystyle }{A}_{\mathrm{1}i}& {\displaystyle }\stackrel{\mathrm{BC}\mathrm{3})}{=}-{\displaystyle \frac{B_{\mathrm{1}i}{b}_{\mathrm{1}i}\cdot \exp \left(b_{\mathrm{1}i}{z}_{\mathrm{bio}}^{-}\right)}{a_{\mathrm{1}i}\cdot \exp \left(a_{\mathrm{1}i}{z}_{\mathrm{bio}}^{-}\right)}}\end{array}$$

See Sect. 2.3.1 for further details on how to find the analytical solution.

Table 2 OM Boundary conditions applied in OMEN-SED. For the boundaries we define $z_{\mathrm{bio}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathrm{bio}}-h)$ and $z_{\mathrm{bio}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathrm{bio}}+h)$.

2.2.2 Oxygen

OMEN-SED explicitly accounts for oxygen consumption by the aerobic degradation of organic matter within the oxic zone and the oxidation of reduced species (i.e. NH₄, H₂S) produced in the anoxic zones of the sediment. In the oxic zone (z<z_ox), aerobic degradation consumes oxygen with a fixed O₂:C ratio (O₂C, Table 10). A predefined fraction, ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, of the ammonium produced during the aerobic degradation of OM is nitrified to nitrate, consuming 2 moles of oxygen per mole of ammonium produced. In addition, OMEN-SED implicitly accounts for oxygen consumption due to the oxidation of reduced species (NH₄, H₂S) produced below the oxic zone through the flux boundary condition at the dynamically calculated oxygen penetration depth z_ox (see Sect. 2.4.2 for details). All oxygen consumption processes can thus be formulated as a function of organic matter degradation. The conservation equation for oxygen is given by

$$\begin{array}{ll}{\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial t}}& {\displaystyle }=\mathrm{0}=D_{{\mathrm{O}}_{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial z}}-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\sum _i{k}_i\\ \text{(12)}& {\displaystyle }& {\displaystyle }\cdot [{\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i]\cdot {\mathrm{POC}}_i\left(z\right).\end{array}$$

For illustrative purposes, we here substitute the analytical solution for the POC depth profile and provide the analytical solution. The remaining paragraphs only outline the general equation, whose analytical solution can be derived in an identical manner. Substituting Eqs. (8) and (9) for POC_i(z) and Eq. (11) for B_1i gives the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(I) Bioturbated zone}(z\le z_{\mathrm{bio}})\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}\stackrel{\mathrm{8}\mathit{\&}\mathrm{11}}{=}{D}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial z}}-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\sum _i{k}_i\\ {\displaystyle }& {\displaystyle }\cdot [{\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i]\cdot (A_{\mathrm{1}i}\cdot [\exp \left(a_{\mathrm{1}i}z\right)-\exp \left(b_{\mathrm{1}i}z\right)]\\ {\displaystyle }& {\displaystyle }+{\mathrm{POC}}_{\mathrm{0}i}\cdot \exp \left(b_{\mathrm{1}i}z\right))\\ {\displaystyle }& {\displaystyle }\text{(II) Non-bioturbated zone}(z_{\mathrm{bio}}<z<z_{\mathrm{ox}})\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}\stackrel{9}{=}{D}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial z}}-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\sum _i{k}_i\\ {\displaystyle }& {\displaystyle }\cdot [{\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i]\cdot (A_{\mathrm{2}i}\cdot \exp (a_{\mathrm{2}i}z\left)\right)\end{array}$$

$D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}$ and $D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}$ denote the O₂ diffusion coefficient for the bioturbated and non-bioturbated zone, respectively. The term $\frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}$ accounts for the volume conversion from solid to dissolved phase and NC_i is the nitrogen to carbon ratio in POC. Integration yields the following analytical solution for each zone.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(I) Bioturbated zone}(z\le z_{\mathrm{bio}}):\\ {\displaystyle }& {\displaystyle }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}\left(z\right)=A_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}}+B_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}}\cdot \exp \left(b_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}}z\right)+\sum _i{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{I}}\cdot \exp \left(a_{\mathrm{1}i}z\right)\\ \text{(13)}& {\displaystyle }& {\displaystyle }+\sum _i{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{II}}\cdot \exp \left(b_{\mathrm{1}i}z\right)+\sum _i{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{III}}\cdot \exp \left(b_{\mathrm{1}i}z\right){\displaystyle }& {\displaystyle }\text{(II) Non-bioturbated zone}(z_{\mathrm{bio}}<z<z_{\mathrm{ox}})\\ {\displaystyle }& {\displaystyle }{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}\left(z\right)=A_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}+B_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}\cdot \exp \left(b_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}z\right)\\ \text{(14)}& {\displaystyle }& {\displaystyle }+\sum _i{\mathrm{\Phi }}_{i,\mathrm{2}}^{\mathrm{I}}\cdot \exp \left(a_{\mathrm{2}i}z\right)\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{b}_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}}={\displaystyle \frac{w}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}}},b_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}={\displaystyle \frac{w}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}}}\\ {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{I}}={\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot {\displaystyle \frac{k_i\cdot ({\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i)\cdot A_{\mathrm{1}i}}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}(-a_{\mathrm{1}i}{)}^{\mathrm{2}}-w\cdot (-a_{\mathrm{1}i})}},\\ {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{II}}=-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot {\displaystyle \frac{k_i\cdot ({\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_{\mathrm{i}})\cdot A_{\mathrm{1}i}}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}(-b_{\mathrm{1}i}{)}^{\mathrm{2}}-w\cdot (-b_{\mathrm{1}i})}}\\ {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{1},i}^{\mathrm{III}}={\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot {\displaystyle \frac{k_i\cdot ({\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i)\cdot {\mathrm{POC}}_{\mathrm{0}i}}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{I}}(-b_{\mathrm{1}i}{)}^{\mathrm{2}}-w\cdot (-b_{\mathrm{1}i})}}\\ {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{i,\mathrm{2}}^{\mathrm{I}}:={\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot {\displaystyle \frac{k_i\cdot ({\mathrm{O}}_{\mathrm{2}}\mathrm{C}+\mathrm{2}{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NC}}_i)\cdot A_{\mathrm{2}i}}{D_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{II}}(-a_{\mathrm{2}i}{)}^{\mathrm{2}}-w\cdot (-a_{\mathrm{2}i})}}\end{array}$$

Determining the four integration constants ($A_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}},B_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{1}},A_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}},B_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}$; see Sect. 2.3 for details) and the a priori unknown oxygen penetration depth requires the definition of five boundary conditions (see Table 3). At the sediment–water interface, OMEN-SED applies a Dirichlet condition (i.e. known concentration) and assumes concentration and flux continuity across the bottom of the bioturbated zone, z_bio. The oxygen penetration depth z_ox marks the lower boundary and is dynamically calculated as the depth at which O₂(z)=0. Therefore, OMEN-SED applies a Dirichlet boundary condition O₂(z_ox)=0. In addition, a flux boundary is applied that implicitly accounts for the oxygen consumption by the partial oxidation of NH₄ and H₂S diffusing into the oxic zone from below (BC 4.2; Table 3). It is assumed that respective fractions (${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$ and ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$) are directly reoxidised at the oxic–anoxic interface and the remaining fraction escapes reoxidation (or is precipitated as pyrite, γ_FeS). OMEN-SED iteratively solves for z_ox by first testing if there is oxygen left at z_max (i.e. O₂(z_max)>0). If that is not the case, it determines the root for the flux boundary condition 4.2 (Table 3). If z_ox=z_max, a zero diffusive flux boundary condition is applied as a lower boundary condition.

Table 3 Boundary conditions for oxygen. For the boundaries we define $z_{\mathrm{bio}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathrm{bio}}-h)$ and $z_{\mathrm{bio}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathrm{bio}}+h)$.

Note: ${\stackrel{\mathrm{̃}}{z}}_{{\mathrm{NO}}_{\mathrm{3}}}=z_{\mathrm{ox}}$ as the upper boundary here, as $z_{{\mathrm{NO}}_{\mathrm{3}}}$ is not known at this point.

2.2.3 Nitrate and ammonium

Nitrogen dynamics in OMEN-SED are controlled by the metabolic production of ammonium, nitrification, denitrification and ammonium adsorption. Ammonium is produced by organic matter degradation in both the oxic and anoxic zones, while denitrification consumes nitrate in the denitrification zone with a fixed NO₃:C ratio (NO₃C; Table 10). Anaerobic ammonium oxidation (anammox) is implicitly included in the model. The organic nitrogen released during denitrification is assumed to be directly oxidised with nitrite to N₂ through a coupling between denitrification and anammox.

The adsorption of ammonium to sediment particles is formulated as an equilibrium process with a constant equilibrium adsorption coefficient $K_{{\mathrm{NH}}_{\mathrm{4}}}$, thus assuming that the adsorption is fast compared to the characteristic timescales of transport processes (Wang and Van Cappellen, 1996). In addition, a defined fraction, ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, of metabolically produced ammonium is directly nitrified to nitrate in the oxic zone, while the nitrification of upward-diffusing ammonium produced in the sulfidic and methanic zones is implicitly accounted for in the boundary conditions. The conservation equations for ammonium and nitrate are thus given by the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(1) Oxic zone}(z\le z_{\mathrm{ox}})\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{NO}}_{\mathrm{3}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{I}}}{\partial z}}\\ \text{(15)}& {\displaystyle }& {\displaystyle }+{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{NC}}_i\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}={\displaystyle \frac{D_{{\mathrm{NH}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z}}\\ \text{(16)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}\cdot {\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{NC}}_i\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(2) Denitrification (or nitrogenous) zone}(z_{\mathrm{ox}}<z\le z_{{\mathrm{NO}}_{\mathrm{3}}})\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{NO}}_{\mathrm{3}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{II}}}{\partial z}}\\ \text{(17)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}{\mathrm{NO}}_{\mathrm{3}}\mathrm{C}\cdot \sum _i{k}_i\cdot {\mathrm{POC}}_i\left(z\right)\text{(18)}& {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}={\displaystyle \frac{D_{{\mathrm{NH}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z}}{\displaystyle }& {\displaystyle }\text{(3) Sulfidic and methanic zone}(z_{{\mathrm{NO}}_{\mathrm{3}}}<z\le z_{\max })\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{III}}}{\partial t}}=\mathrm{0}={\displaystyle \frac{D_{{\mathrm{NH}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{III}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{NH}}_{\mathrm{4}}}^{\mathrm{III}}}{\partial z}}\\ \text{(19)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}\cdot {\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{NC}}_i\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right)\end{array}$$

$D_{{\mathrm{NO}}_{\mathrm{3}}}$ and $D_{{\mathrm{NH}}_{\mathrm{4}}}$ denote the diffusion coefficients for NO₃ and NH₄, which depend on the bioturbation status of the respective geochemical zone (compare to Sect. 2.3.1). Integration of Eqs. (15)–(19) yields the analytical solutions, which are not further developed here but follow the procedure outlined in Sect. 2.2.2 for oxygen (also see Sect. 2.3.1 for more details on how to find the analytical solution). Table 4 summarises the boundary conditions applied in OMEN-SED to solve Eqs. (15)–(19) and to find the a priori unknown nitrate penetration depth, $z_{{\mathrm{NO}}_{\mathrm{3}}}$. The model assumes known bottom water concentrations for both NO₃ and NH₄, the complete consumption of nitrate at the nitrate penetration depth (in the case $z_{{\mathrm{NO}}_{\mathrm{3}}}<z_{\max }$) and no change in ammonium flux at z_max. In addition, concentration and diffusive flux continuity across z_bio and z_ox is considered for NO₃ and NH₄. Furthermore, the reoxidation of upward-diffusing reduced ammonium is accounted for in the oxic–anoxic boundary condition for nitrate and ammonium. OMEN-SED iteratively solves for $z_{{\mathrm{NO}}_{\mathrm{3}}}$ by first testing if there is nitrate left at z_max (i.e. NO₃(z_max)>0) and otherwise by finding the root for the flux boundary condition 6.2 (Table 4).

Table 4 Boundary conditions for nitrate and ammonium. For the boundaries we define $z_{\mathit{\_}\mathit{\_}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}-h)$ and $z_{\mathit{\_}\mathit{\_}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}+h)$.

2.2.4 Sulfate and sulfide

Below the denitrification zone ($z>z_{{\mathrm{NO}}_{\mathrm{3}}}$), organic matter degradation is coupled to sulfate reduction, consuming sulfate and producing hydrogen sulfide with a fixed SO₄:C ratio (SO₄C; Table 10). In addition, the anaerobic oxidation of upward-diffusing methane (AOM) produced below the sulfate penetration and the associated consumption of sulfate and production of sulfide, as well as the production of sulfate and consumption of sulfide through sulfide oxidation, are implicitly accounted for through the boundary conditions (Table 5). In the sulfidic zone a defined fraction of sulfide, γ_FeS, can be precipitated as pyrite (in the presented simulations γ_FeS=0.0 as we do not want to make any assumptions about pyrite precipitation). It can be safely assumed that almost all CH₄ is oxidised anaerobically in the sediments (e.g. Reeburgh, 2007, suggests up to 90 %), except for active (very localised) sites and slope failure, which can, in theory, be accounted for through the ${\mathit{\gamma }}_{{\mathrm{CH}}_{\mathrm{4}}}$ term. The conservation equations for sulfate and sulfide are thus given by the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(1) Oxic and nitrogenous zone}(z\le z_{{\mathrm{NO}}_{\mathrm{3}}})\\ \text{(20)}& {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{SO}}_{\mathrm{4}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z}}\text{(21)}& {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{I}}}{\partial z}}{\displaystyle }& {\displaystyle }\text{(2) Sulfidic zone}(z_{{\mathrm{NO}}_{\mathrm{3}}}<z\le z_{{\mathrm{SO}}_{\mathrm{4}}})\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{SO}}_{\mathrm{4}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z}}\\ \text{(22)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{SO}}_{\mathrm{4}}\mathrm{C}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{II}}}{\partial z}}\\ \text{(23)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot (\mathrm{1}-{\mathit{\gamma }}_{\mathrm{FeS}})\cdot \sum _i{\mathrm{SO}}_{\mathrm{4}}\mathrm{C}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(3) Methanic zone}(z_{{\mathrm{SO}}_{\mathrm{4}}}<z\le z_{\max })\\ \text{(24)}& {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{III}}}{\partial t}}=\mathrm{0}=D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{III}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{\mathrm{III}}}{\partial z}}\end{array}$$

$D_{{\mathrm{SO}}_{\mathrm{4}}}$ and $D_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$ denote the diffusion coefficients for SO₄ and H₂S, which depend on the bioturbation status of the respective geochemical zone (compare to Sect. 2.3.1). Integration of Eqs. (20)–(24) yields the analytical solution and Table 5 summarises the boundary conditions applied. OMEN-SED assumes known concentrations at the sediment–water interface and continuity across the bioturbation depth and the nitrate penetration depth. The reoxidation of reduced H₂S to SO₄ is accounted for implicitly via the oxic–anoxic boundary condition for both species, while the reduction of SO₄ and the associated production of H₂S via AOM is accounted for through the respective boundary conditions at $z_{{\mathrm{SO}}_{\mathrm{4}}}$. In the case $z_{{\mathrm{SO}}_{\mathrm{4}}}<z_{\max }$, OMEN-SED assumes zero sulfate concentration at $z_{{\mathrm{SO}}_{\mathrm{4}}}$ and its diffusive flux must equal the amount of methane produced below (with a methane to carbon ratio of MC); in the case $z_{{\mathrm{SO}}_{\mathrm{4}}}=z_{\max }$, a zero diffusive flux condition for sulfate is considered. OMEN-SED iteratively solves for $z_{{\mathrm{SO}}_{\mathrm{4}}}$ by first testing if there is sulfate left at z_max (i.e. SO₄(z_max)>0) and otherwise by finding the root for the flux boundary condition 8.2 (Table 5). At the lower boundary z_max zero diffusive flux of H₂S is considered.

Table 5 Boundary conditions for sulfate and sulfide. For the boundaries we define $z_{\mathit{\_}\mathit{\_}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}-h)$ and $z_{\mathit{\_}\mathit{\_}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}+h)$.

2.2.5 Phosphate

Figure 2 A schematic of the sedimentary P cycle in OMEN-SED. Red numbers represent kinetic rate constants for phosphorus dynamics (compare to Table 10; p_i represents the uptake rate of PO₄ via primary production in shallow environments). Adapted from Slomp et al. (1996).

The biogeochemical description of phosphorus (P) dynamics builds on earlier models developed by Slomp et al. (1996) and accounts for phosphorus recycling through organic matter degradation, adsorption onto sediments and iron(III) hydroxides (Fe-bound P), as well as carbonate fluorapatite (CFA or authigenic P) formation (see Fig. 2 for a schematic overview of the sedimentary P cycle). In the oxic zone of the sediment, PO₄ liberated through organic matter degradation can adsorb to iron(III) hydroxides forming Fe-bound P (or FeP; Slomp et al., 1998). Below the oxic zone, PO₄ is not only produced via organic matter degradation but can also be released from the Fe-bound P pool due to the reduction of iron(III) hydroxides under anoxic conditions. Furthermore, in these zones phosphate concentrations build up and pore waters can thus become supersaturated with respect to carbonate fluorapatite, thus triggering the authigenic formation of CFA (Van Cappellen and Berner, 1988). Phosphorus bound in these authigenic minerals represents a permanent sink for reactive phosphorus (Slomp et al., 1996). As for ammonium, the adsorption of P to the sediment matrix is treated as an equilibrium process parameterised with dimensionless adsorption coefficients for the oxic and anoxic zone, respectively ($K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}$, $K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}$; Slomp et al., 1998). The adsorption and desorption of P to iron(III) hydroxides and authigenic fluorapatite formation are described as first-order reactions with rate constants k_s, k_m and k_a, respectively (Table 10). The rate of the respective process is calculated as the product of the rate constant and the difference between the current concentration (of PO₄ and FeP) and an equilibrium or asymptotic concentration (Slomp et al., 1996). The asymptotic Fe-bound P concentration is FeP^∞ and the equilibrium concentrations for P sorption and authigenic fluorapatite formation are ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{s}}$ and ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{a}}$, respectively (Table 10). The last term in Eqs. (25) and (26) represents the sorption of PO₄ to FeP in the oxic zone, the last term in Eqs. (27) and (28) is the release of PO₄ from the FeP pool and the fourth term in Eq. (28) represents the permanent loss of PO₄ to authigenic fluorapatite formation. The conservation equations for phosphate and Fe-bound P are thus given by the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(1) Oxic zone}(z\le z_{\mathrm{ox}})\\ {\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial t}}={\displaystyle \frac{D_{{\mathrm{PO}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{I}}}{\partial z}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}}}\sum _i\left({\mathrm{PC}}_i\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right)\right)\\ \text{(25)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{k_s}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}}}({{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{I}}-{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{s}}){\displaystyle }& {\displaystyle \frac{\partial {\mathrm{FeP}}^{\mathrm{I}}}{\partial t}}=D_{\mathrm{FeP}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{FeP}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{FeP}}^{\mathrm{I}}}{\partial z}}\\ \text{(26)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathit{\varphi }}{\mathrm{1}-\mathit{\varphi }}}{k}_s({{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{I}}-{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{s}}){\displaystyle }& {\displaystyle }\text{(2) Anoxic zones}(z_{\mathrm{ox}}<z\le z_{\max })\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{FeP}}^{\mathrm{II}}}{\partial t}}=D_{\mathrm{FeP}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{FeP}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{FeP}}^{\mathrm{II}}}{\partial z}}\\ \text{(27)}& {\displaystyle }& {\displaystyle }-k_m({\mathrm{FeP}}^{\mathrm{II}}-{\mathrm{FeP}}^{\mathrm{\infty }}){\displaystyle }& {\displaystyle \frac{\partial {{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial t}}={\displaystyle \frac{D_{{\mathrm{PO}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{II}}}{\partial z}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}}}\sum _i\left({\mathrm{PC}}_i\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right)\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{k_a}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}}}({{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{II}}-{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{a}})\\ \text{(28)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{(\mathrm{1}-\mathit{\varphi })}{\mathit{\varphi }}}{\displaystyle \frac{k_m}{\mathrm{1}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}}}({\mathrm{FeP}}^{\mathrm{II}}-{\mathrm{FeP}}^{\mathrm{\infty }})\end{array}$$

$D_{{\mathrm{PO}}_{\mathrm{4}}}$ denotes the diffusion coefficient for PO₄, which depends on the bioturbation status of the respective geochemical zone, and D_FeP=D_bio for z≤z_bio and D_FeP=0 for z>z_bio (compare to Sect. 2.3.1). Integration of Eqs. (25)–(28) yields the analytical solution and Table 6 summarises the boundary conditions applied in OMEN-SED. The model assumes known bottom water concentrations and equal concentrations and diffusive fluxes at z_bio and z_ox for both species. Additionally, OMEN-SED considers no change in phosphate flux and an asymptotic Fe-bound P concentration at z_max.

Table 6 Boundary conditions for phosphate and Fe-bound P (FeP). For the boundaries we define $z_{\mathit{\_}\mathit{\_}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}-h)$ and $z_{\mathit{\_}\mathit{\_}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}+h)$.

2.2.6 Dissolved inorganic carbon (DIC)

OMEN-SED accounts for the production of dissolved inorganic carbon (DIC) through organic matter degradation and methane oxidation. Organic matter degradation produces dissolved inorganic carbon with a stoichiometric DIC:C ratio of 1:2 in the methanic zone and 1:1 in the rest of the sediment column (DICC^II and DICC^I, respectively). DIC production through methane oxidation is implicitly taken into account through the boundary condition at $z_{{\mathrm{SO}}_{\mathrm{4}}}$. A mechanistic description of DIC production from CaCO₃ dissolution would lead to significant mathematical problems and is therefore not included in the current version of OMEN-SED. The conservation equations for DIC are thus given by the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(1) Oxic, nitrogenous and sulfidic zone}(z\le z_{{\mathrm{SO}}_{\mathrm{4}}})\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{DIC}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}=D_{\mathrm{DIC}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{DIC}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{DIC}}^{\mathrm{I}}}{\partial z}}+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\\ \text{(29)}& {\displaystyle }& {\displaystyle }\cdot \sum _i{\mathrm{DICC}}^{\mathrm{I}}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(2) Methanic zone}(z_{{\mathrm{SO}}_{\mathrm{4}}}<z\le z_{\max })\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{DIC}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}=D_{\mathrm{DIC}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{DIC}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{DIC}}^{\mathrm{II}}}{\partial z}}+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\\ \text{(30)}& {\displaystyle }& {\displaystyle }\cdot \sum _i{\mathrm{DICC}}^{\mathrm{II}}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right)\end{array}$$

D_DIC denotes the diffusion coefficient for DIC (taking the values for ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ from Schulz, 2006), which depends on the bioturbation status of the respective geochemical zone. Integration of Eqs. (29) and (30) yields the analytical solution and Table 7 summarises the boundary conditions applied in OMEN-SED. A Dirichlet condition is applied at the sediment–water interface. In addition, the model assumes a zero diffusive flux through the lower boundary z_max and continuity across the bottom of the bioturbated zone and the sulfate penetration depth. An additional flux boundary condition at $z_{{\mathrm{SO}}_{\mathrm{4}}}$ implicitly accounts for DIC production through the anaerobic oxidation of methane (Table 7 Eq. 5).

Table 7 Boundary conditions for DIC. For the boundaries we define $z_{\mathit{\_}\mathit{\_}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}-h)$ and $z_{\mathit{\_}\mathit{\_}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}+h)$.

2.2.7 Alkalinity

Organic matter degradation and secondary redox reactions exert a complex influence on alkalinity (e.g. Jourabchi et al., 2005; Wolf-Gladrow et al., 2007; Krumins et al., 2013). To model alkalinity, OMEN-SED divides the sediment column into four geochemical zones in which different equations describe the biogeochemical processes using variable stoichiometric coefficients (compare values in Table 10). Above z_ox, the combined effects of NH₄ and P release due to aerobic OM degradation increases alkalinity according to ALK^OX, whereas nitrification decreases alkalinity with stoichiometry ALK^NIT. In the remaining three zones anaerobic OM degradation generally results in an increase in alkalinity, with the exact magnitude depending on the nature of the terminal electron acceptor used (i.e. ALK^DEN, ALK^SUL, ALK^MET). In addition, the effect of secondary redox reactions, such as nitrification, sulfide and methane oxidation, and pyrite precipitation are implicitly accounted for in the boundary conditions. Note that the alkalinity description in the current version of OMEN-SED does not account for CaCO₃ dissolution and precipitation due to the mathematical complexity of the problem. In OMEN-SED, the conservation equations for alkalinity are thus given by the following.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{(1) Oxic zone}(z\le z_{\mathrm{ox}})\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{I}}}{\partial t}}=\mathrm{0}=D_{\mathrm{ALK}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{ALK}}^{\mathrm{I}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{I}}}{\partial z}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i\left({\mathrm{ALK}}^{\mathrm{NIT}}\cdot {\displaystyle \frac{{\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{NH}}_{\mathrm{4}}}}}{\mathrm{NC}}_i+{\mathrm{ALK}}^{\mathrm{OX}}\right)\\ \text{(31)}& {\displaystyle }& {\displaystyle }\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(2) Denitrification or nitrogenous zone}(z_{\mathrm{ox}}<z\le z_{{\mathrm{NO}}_{\mathrm{3}}})\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{II}}}{\partial t}}=\mathrm{0}=D_{\mathrm{ALK}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{ALK}}^{\mathrm{II}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{II}}}{\partial z}}\\ \text{(32)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{ALK}}^{\mathrm{DEN}}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(3) Sulfidic zone}(z_{{\mathrm{NO}}_{\mathrm{3}}}<z\le z_{{\mathrm{SO}}_{\mathrm{4}}})\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{III}}}{\partial t}}=\mathrm{0}=D_{\mathrm{ALK}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{ALK}}^{\mathrm{III}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{III}}}{\partial z}}\\ \text{(33)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{ALK}}^{\mathrm{SUL}}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right){\displaystyle }& {\displaystyle }\text{(4) Methanic zone}(z_{{\mathrm{SO}}_{\mathrm{4}}}<z\le z_{\max })\\ {\displaystyle }& {\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{IV}}}{\partial t}}=\mathrm{0}=D_{\mathrm{ALK}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathrm{ALK}}^{\mathrm{IV}}}{\partial z^{\mathrm{2}}}}-w{\displaystyle \frac{\partial {\mathrm{ALK}}^{\mathrm{IV}}}{\partial z}}\\ \text{(34)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}-\mathit{\varphi }}{\mathit{\varphi }}}\cdot \sum _i{\mathrm{ALK}}^{\mathrm{MET}}\cdot k_i\cdot {\mathrm{POC}}_i\left(z\right)\end{array}$$

D_ALK denotes the diffusion coefficient for alkalinity (taking the values for ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ from Schulz, 2006), which depends on the bioturbation status of the respective geochemical zone. Integration of Eqs. (31)–(34) yields the analytical solution and Table 8 summarises the boundary conditions applied in OMEN-SED. A Dirichlet boundary condition is applied at the sediment–water interface. The decrease in alkalinity due to the oxidation of reduced species produced in the anoxic zones and due to the precipitation of pyrite (with stoichiometry ALK^NIT, ${\mathrm{ALK}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$ and ALK^FeS) is implicitly taken into account through the flux boundary condition at z_ox (Table 8 Eq. 5). Furthermore, the oxidation of methane by sulfate reduction increases alkalinity with stoichiometry ALK^AOM, which is accounted for through the flux boundary condition at $z_{{\mathrm{SO}}_{\mathrm{4}}}$ (Table 8 Eq. 9). At the lower boundary z_max a zero diffusive flux condition is applied.

Table 8 Boundary conditions for alkalinity. For the boundaries we define $z_{\mathit{\_}\mathit{\_}}^{-}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}-h)$ and $z_{\mathit{\_}\mathit{\_}}^{+}={lim}_{h\to \mathrm{0}}(z_{\mathit{\_}\mathit{\_}}+h)$.

2.3 Determination of integration constants

The integration constants of all general analytical solutions derived above change in response to changing boundary conditions. Thus, OMEN-SED has to redetermine integration constants for each dynamic zone (i.e. $z_{\mathrm{ox}},z_{\mathrm{bio}},z_{{\mathrm{NO}}_{\mathrm{3}}}$ and $z_{{\mathrm{SO}}_{\mathrm{4}}}$) at every time step for all biogeochemical species. The bioturbation boundary poses a particular challenge as it can theoretically occur in any of the dynamic geochemical zones (Fig. 3). Therefore, in order to generalise and simplify this recurring boundary matching problem, an independent, generic algorithm (generic boundary condition matching) is implemented (rather than using multiple fully-worked-out algebraic solutions for each possible case and every biogeochemical species). As a consequence, the algorithm only has to solve a two-simultaneous-equation problem.

2.3.1 Generic boundary condition matching (GBCM)

As discussed in Sect. 2.1, the solution to the general steady-state transport-reaction equation (Eq. 2) for a generic species C is of the general form

$$\begin{array}{ll}{\displaystyle }C\left(z\right)& {\displaystyle }=A\cdot \exp \left(az\right)+B\cdot \exp \left(bz\right)+\sum _j{\displaystyle \frac{{\mathit{\alpha }}_j}{D{\mathit{\beta }}_j^{\mathrm{2}}-w{\mathit{\beta }}_j-k}}\\ \text{(35)}& {\displaystyle }& {\displaystyle }\cdot \exp (-{\mathit{\beta }}_{j}z)+{\displaystyle \frac{Q}{k}}\end{array}$$

and can therefore be expressed as

$$\begin{array}{}\text{(36)}& {\displaystyle }C\left(z\right)=A\cdot E\left(z\right)+B\cdot F\left(z\right)+G\left(z\right),\end{array}$$

where E(z) and F(z) are the homogeneous solutions to the ODE, G(z) the particular integral (collectively called the basis functions), and A and B are the integration constants that must be determined using the boundary conditions (shown in Fig. 3 for the whole sediment column).

Each internal boundary matching problem (i.e. excluding z=0 and z=z_max) involves matching continuity and flux for the two solutions to the respective reaction–transport equation above (C_U(z) “upper”) and below (C_L(z) “lower”) the dynamic boundary at z=z_b.

$$\begin{array}{}\text{(37)}& {\displaystyle }{C}_{\mathrm{U}}\left(z\right)& {\displaystyle }=A_{\mathrm{U}}\cdot E_{\mathrm{U}}\left(z\right)+B_{\mathrm{U}}\cdot F_{\mathrm{U}}\left(z\right)+G_{\mathrm{U}}\left(z\right)\text{(38)}& {\displaystyle }{C}_{\mathrm{L}}\left(z\right)& {\displaystyle }=A_{\mathrm{L}}\cdot E_{\mathrm{L}}\left(z\right)+B_{\mathrm{L}}\cdot F_{\mathrm{L}}\left(z\right)+G_{\mathrm{L}}\left(z\right)\end{array}$$

OMEN-SED generally applies concentration continuity and flux boundary conditions at its internal dynamic boundaries.

Continuity (where for generality we allow a discontinuity V_b):

$$\begin{array}{ll}\text{(39)}& {\displaystyle }& {\displaystyle }{C}_{\mathrm{U}}\left(z_{\mathrm{b}}\right)=C_{\mathrm{L}}\left(z_{\mathrm{b}}\right)+V_b.{\displaystyle }& {\displaystyle }\text{Flux}:\\ \text{(40)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{U}}{C}_{\mathrm{U}}^{\prime }\left(z_{\mathrm{b}}\right)+w{C}_{\mathrm{U}}\left(z_{\mathrm{b}}\right)=D_{\mathrm{L}}{C}_{\mathrm{L}}^{\prime }\left(z_{\mathrm{b}}\right)+w{C}_{\mathrm{L}}\left(z_{\mathrm{b}}\right)+F_b,\end{array}$$

where w is advection, D represents the diffusion coefficients and F_b is any flux discontinuity (e.g. resulting from secondary redox reactions). Considering that the advective flux above and below the boundary is equal (i.e. wC_U(z_b)=wC_L(z_b)) and substituting the general ODE solutions (37) and (38), the boundary conditions can be represented as two equations connecting the four integration constants:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(\begin{array}{cc}{E}_{\mathrm{U}}& F_{\mathrm{U}}\\ D_{\mathrm{U}}{E}_{\mathrm{U}}^{\prime }& D_{\mathrm{U}}{F}_{\mathrm{U}}^{\prime }\end{array}\right)\left(\begin{array}{c}{A}_{\mathrm{U}}\\ B_{\mathrm{U}}\end{array}\right)=\left(\begin{array}{cc}{E}_{\mathrm{L}}& F_{\mathrm{L}}\\ D_{\mathrm{L}}{E}_{\mathrm{L}}^{\prime }& D_{\mathrm{L}}{F}_{\mathrm{L}}^{\prime }\end{array}\right)\left(\begin{array}{c}{A}_{\mathrm{L}}\\ B_{\mathrm{L}}\end{array}\right)\\ \text{(41)}& {\displaystyle }& {\displaystyle }+\left(\begin{array}{c}{G}_{\mathrm{L}}-G_{\mathrm{U}}+V_b\\ D_{\mathrm{L}}{G}_{\mathrm{L}}^{\prime }-D_{\mathrm{U}}{G}_{\mathrm{U}}^{\prime }+F_b-w{V}_b\end{array}\right),\end{array}$$

where the ODE solutions E, F and G are all evaluated at z_b. Equation (41) can now be solved to give A_U and B_U as a function of the integration constants from the layer below (A_L and B_L), thereby constructing a piecewise solution for both layers with just two integration constants (this is implemented in the function benthic_utils.matchsoln of OMEN-SED):

$$\begin{array}{}\text{(42)}& {\displaystyle }\left(\begin{array}{c}{A}_{\mathrm{U}}\\ B_{\mathrm{U}}\end{array}\right)=\left(\begin{array}{cc}{c}_{\mathrm{1}}& c_{\mathrm{2}}\\ c_{\mathrm{3}}& c_{\mathrm{4}}\end{array}\right)\left(\begin{array}{c}{A}_{\mathrm{L}}\\ B_{\mathrm{L}}\end{array}\right)+\left(\begin{array}{c}{d}_{\mathrm{1}}\\ d_{\mathrm{2}}\end{array}\right).\end{array}$$

Using Eq. (42), C_U(z) in Eq. (37) can now be rewritten as a function of A_L and B_L (implemented in benthic_utils.xformsoln):

$$\begin{array}{ll}{\displaystyle }{C}_{\mathrm{U}}\left(z\right)=& {\displaystyle }(c_{\mathrm{1}}{A}_{\mathrm{L}}+c_{\mathrm{2}}{B}_{\mathrm{L}}+d_{\mathrm{1}})\cdot E_{\mathrm{U}}\left(z\right)+(c_{\mathrm{3}}{A}_{\mathrm{L}}+c_{\mathrm{4}}{B}_{\mathrm{L}}+d_{\mathrm{2}})\\ \text{(43)}& {\displaystyle }& {\displaystyle }\cdot F_{\mathrm{U}}\left(z\right)+G_{\mathrm{U}}\left(z\right),\end{array}$$

and hence the “transformed” basis functions $E_{\mathrm{U}}^{*}\left(z\right),F_{\mathrm{U}}^{*}\left(z\right)$ and $G_{\mathrm{U}}^{*}\left(z\right)$ can be defined such that

$$\begin{array}{}\text{(44)}& {\displaystyle }{C}_{\mathrm{U}}\left(z\right)=A_{\mathrm{L}}\cdot E_{\mathrm{U}}^{*}\left(z\right)+B_{\mathrm{L}}\cdot F_{\mathrm{U}}^{*}\left(z\right)+G_{\mathrm{U}}^{*}\left(z\right),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }{E}_{\mathrm{U}}^{*}\left(z\right)& {\displaystyle }=c_{\mathrm{1}}{E}_{\mathrm{U}}\left(z\right)+c_{\mathrm{3}}{F}_{\mathrm{U}}\left(z\right)\\ \text{(45)}& {\displaystyle }{F}_{\mathrm{U}}^{*}\left(z\right)& {\displaystyle }=c_{\mathrm{2}}{E}_{\mathrm{U}}\left(z\right)+c_{\mathrm{4}}{F}_{\mathrm{U}}\left(z\right){\displaystyle }{G}_{\mathrm{U}}^{*}\left(z\right)& {\displaystyle }=G_{\mathrm{U}}\left(z\right)+d_{\mathrm{1}}{E}_{\mathrm{U}}\left(z\right)+d_{\mathrm{2}}{F}_{\mathrm{U}}\left(z\right).\end{array}$$

Figure 3 Schematic of the generic boundary condition matching (GBCM) problem. Shown are the resulting integration constants (A_i, B_i) and ODE solutions ($E_i,F_i,G_i$) for the different sediment layers and the bioturbation boundary (possible locations are indicated by the green vertical arrow).

Equations (42), (44) and (45) can now be consecutively applied for each of the dynamic biogeochemical zone boundaries (Fig. 3) starting at the bottom of the sediment column. The net result is a piecewise solution for the whole sediment column with just two integration constants (coming from the lowest layer), which can then be solved for by applying the boundary conditions at the sediment–water interface and the bottom of the sediments.

2.3.2 Abstracting out the bioturbation boundary

The bioturbation boundary affects the diffusion coefficient of the modelled solutes and the conservation equation of organic matter (and thereby the exact form of each reaction–transport equation). This boundary is particularly inconvenient as it can in principle occur in the middle of any of the dynamically shifting biogeochemical zones and therefore generate multiple cases (Fig. 3). The GBCM algorithm described above is thus not only used to construct a piecewise solution for the whole sediment column, but also to abstract out the bioturbation boundary. For each biogeochemical zone the “bioturbation status” is initially tested (i.e. fully bioturbated, fully non-bioturbated or crossing the bioturbation boundary). Therefore, the upper and lower boundaries for the different zones (e.g. for the nitrogenous zone z_U=zox, $z_{\mathrm{L}}=z_{{\mathrm{NO}}_{\mathrm{3}}}$) and the respective reactive terms and diffusion coefficients (bioturbated and non-bioturbated) are passed over to the routine zTOC.prepfg_l12 in which the bioturbation status is determined. In the case that the bioturbation depth is located within this zone (i.e. $z_{\mathrm{U}}<z_{\mathrm{bio}}<z_{\mathrm{L}}$) a piecewise solution for this layer is constructed. Therefore, the reactive terms and diffusion coefficients are handed over to the routines zTOC.calcfg_l1 and zTOC.calcfg_l2, which calculate the basis functions ($E_{\mathrm{U}},F_{\mathrm{U}},G_{\mathrm{U}}$ and $E_{\mathrm{L}},F_{\mathrm{L}},G_{\mathrm{L}}$) and their derivatives for the bioturbated and the non-bioturbated part of this specific geochemical zone. The concentration and flux for both solutions at z_bio are matched and the coefficients $c_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},c_{\mathrm{4}},d_{\mathrm{1}}$ and d₂ (as in Eq. 42) are calculated by the routine benthic_utils.matchsoln. These coefficients and the “bioturbation status” of the layer are passed back to the main GBCM algorithm in which they can be used by the routine benthic_utils.xformsoln to calculate the “transformed” basis functions ($E_{\mathrm{U}}^{*}\left(z\right),F_{\mathrm{U}}^{*}\left(z\right),G_{\mathrm{U}}^{*}\left(z\right)$) such that both layers are expressed in the same basis (compare Eqs. 43–45).

For instance, in the case of sulfate, zTOC.prepfg_l12 is called three times before the actual profile is calculated (once per zone: oxic, nitrogenous, sulfidic) and hands back the information about the “bioturbation status” of the three layers and the coefficients $c_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},c_{\mathrm{4}},d_{\mathrm{1}}$ and d₂ for the biogeochemical zone including the bioturbation depth. When calculating the complete piecewise solution for the sediment column, this information is passed to the function zTOC.calcfg_l12, which sorts out the correct solution type to use. The main GBCM algorithm therefore never needs to know whether it is dealing with a piecewise solution (i.e. matched across the bioturbation boundary) or a “simple” solution (i.e. the layer is fully bioturbated or fully non-bioturbated).

2.4 Model parameters

The following section provides a summary of global relationships used to constrain reaction and transport parameters in OMEN-SED. Table 9 synthesises sediment and transport parameters, while Table 10 provides an overview of all biogeochemical parameters used in OMEN-SED.

2.4.1 Transport parameters

The burial of sediments and pore water is directly related to the accumulation of new material on the sea floor (i.e. sedimentation, Burdige, 2006). This results in a downward advective flux of older sediment material and pore water in relation to the sediment–water interface. When coupled to an ocean model, its sedimentation flux can be readily used in OMEN-SED. The stand-alone version of OMEN-SED uses the empirical global relationship between the sediment accumulation rate (w in cm yr⁻¹) and sea-floor depth (z in m) of Burwicz et al. (2011):

$$\begin{array}{}\text{(46)}& {\displaystyle }w={\displaystyle \frac{w_{\mathrm{1}}}{\mathrm{1}+(\frac{z}{z\mathrm{1}}{)}^{c\mathrm{1}}}}+{\displaystyle \frac{w_{\mathrm{2}}}{\mathrm{1}+(\frac{z}{z\mathrm{2}}{)}^{c\mathrm{2}}}},\end{array}$$

with parameter values as found in the original study (i.e. w₁=0.117 cm yr⁻¹, w₂=0.006 cm yr⁻¹, z1=200 m, z2=4000 m, c1=3, c2=10). As an option we include the parameterisation of Middelburg et al. (1997):

$$\begin{array}{}\text{(47)}& {\displaystyle }w=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{0.87478367}-\mathrm{0.00043512}\cdot z}.\end{array}$$

As mentioned before (Sect. 2.1), the diffusion coefficient of species i is calculated as $D_i=D_{i,\mathrm{0}}+D_{\mathrm{bio}}=D_{\mathrm{mol},i}\cdot f_{\mathrm{ir}}+D_{\mathrm{bio}}$ for dissolved species and D_i=D_bio for solid species. The bioturbation coefficient D_bio (cm² yr⁻¹) is constant in the bioturbated zone and also follows the empirical relationship by Middelburg et al. (1997):

$$\begin{array}{}\text{(48)}& {\displaystyle }{D}_{\mathrm{bio}}=\mathrm{5.2}\times {\mathrm{10}}^{\mathrm{0.76241122}-\mathrm{0.00039724}\cdot z}.\end{array}$$

Observations indicate that bioturbation is largely restricted to the upper 10 cm of the sediments and is only marginally related to sea-floor depth (e.g. Boudreau, 1998; Teal et al., 2010). Therefore, OMEN-SED imposes a globally invariant bioturbation depth z_bio of 10 cm. In the case that the bottom water oxygen concentration is low (here < 4.5 nmol cm⁻³, which is often used to define suboxic waters; e.g. Morrison et al., 1999; Karstensen et al., 2008) infaunal activity is assumed to cease and z_bio=0.01 cm. We choose a low value unequal to zero in order to simplify the implementation of the model. This approach ensures that the sediment column always consists of a bioturbated (even though very small for the low oxygen condition) and a non-bioturbated zone, and thus the same GBCM algorithm can be used to solve the conservation equations. Furthermore, when OMEN-SED is coupled to an Earth system model the same method can be used to convert the POC depositional flux into an SWI concentration (i.e. the flux needs to be converted assuming bioturbation; see Sect. 4.1).

Bio-irrigation (i.e. the pumping activity by burrow-dwelling animals) exchanges burrow water with overlying water and may enhance the SWI flux of solutes (Aller, 1984, 1988). Several approaches exist to incorporate this into a 1-D diagenetic model, for instance as a non-local transport–exchange process (Boudreau, 1984; Emerson et al., 1984) or as an enhancement factor of the molecular diffusion coefficient (Devol and Christensen, 1993; Soetaert et al., 1996). In OMEN-SED the latter approach is applied and the apparent “bio-diffusion” coefficient is calculated as $D_{i,\mathrm{0}}=D_{\mathrm{mol},i}\cdot f_{\mathrm{ir}}$. Soetaert et al. (1996) derived an empirical relationship between f_ir and sea-floor depth ($f_{\mathrm{ir}}=\mathrm{Min}\mathit{\{}\mathrm{1};\mathrm{15.9}\cdot z^{-\mathrm{0.43}}\mathit{\}}$) based on observations from Archer and Devol (1992) and Devol and Christensen (1993), which is used in OMEN-SED. The specific molecular diffusion coefficients D_mol,i are corrected for sediment porosity ϕ and tortuosity F and are linearly interpolated for an ambient temperature T (in ^∘C) using zero-degree coefficients $D_i^{\mathrm{0}}$ and temperature-dependent diffusion coefficients $D_i^T$ (Soetaert et al., 1996):

$$D_{\mathrm{mol},i}=(D_i^{\mathrm{0}}+D_i^T\cdot T)\cdot {\displaystyle \frac{\mathrm{1}}{\mathit{\varphi }\cdot F}}.$$

Tortuosity can be expressed in terms of porosity as $F=\frac{\mathrm{1}}{{\mathit{\varphi }}^m}$ (Ullman and Aller, 1982) with the exponent m varying according to the type of sediment (here m=3 is used representing muddy sediments with high porosity). Values for $D_i^T$ and $D_i^{\mathrm{0}}$ are summarised in Table 9 and are adapted from Li and Gregory (1974), Schulz (2006) and Gypens et al. (2008).

Middelburg et al. (1997)Boudreau (1998)Teal et al. (2010)Middelburg et al. (1997)Soetaert et al. (1996)(Li and Gregory, 1974; Schulz, 2006; Gypens et al., 2008)

Table 9 Sediment characteristics and transport parameters.

Note: DIC and ALK coefficients are the values of ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ from Schulz (2006).

2.4.2 Stoichiometries and reaction parameters

The first-order organic matter degradation constants of compound class i, k_i (yr⁻¹), are assumed invariant along the sediment column and therefore independent of the nature of the terminal electron acceptor. The rate constants can be altered manually to fit observed sediment profiles (compare modelled profiles in Sect. 3.2) or related to a master variable provided by a coupled Earth system model (e.g. sedimentation rate; see Sect. 4.2). The partitioning of the bulk OM pool into reactivity classes (f_i) needs to be specified manually in the stand-alone version or can be provided by the ESM. Organic matter degradation releases N, P and DIC to the pore water using Redfield molar ratios (Redfield, 1963) and consumes TEA with specific stoichiometries (O₂C, NO₃C, SO₄C) as summarised in Table 10. Table A1 in the Appendix provides a list of reactions and their stoichiometries as implemented in OMEN-SED. The effect of OM degradation, secondary redox reactions and pyrite precipitation on total alkalinity is also accounted for via reaction-specific stoichiometries representing the release of NH₄, H₂S and P and is based on Jourabchi et al. (2005). In reality, the reoxidation of reduced substances produced during OM degradation may be incomplete. Yet, in OMEN-SED we have to assume their complete, instantaneous reoxidation at z_ox to allow for an analytical solution. In order to relax this assumption for cases in which it can be justified, we include a “switch” to allow part of the NH₄, H₂S and CH₄ flux to escape reoxidation. The secondary redox parameters (i.e. ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$, ${\mathit{\gamma }}_{{\mathrm{CH}}_{\mathrm{4}}}$) therefore account for the fraction of reduced substances that are reoxidised and would ideally be parameterised, for instance, in relation to bottom water oxygen concentration or oxygen penetration depth (z_ox). Gypens et al. (2008), for example, expressed ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$ as a function of oxygen penetration depth (${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}=\mathrm{0.243}\cdot \ln \left(z_{\mathrm{ox}}\right)+\mathrm{1.8479}$) based on a fitting exercise to a numerical model and showed that the fraction varies between 0.2 for z_ox=0.1 cm and 1.0 for z_ox>3 cm. Due to mathematical constraints in OMEN-SED for finding an analytical solution to the model equations these fractions take constant values generally representing oxygenated deep sea conditions. However, when coupled to an ESM, ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$ becomes dependent on the bottom water oxygenation state. That is, ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\mathrm{1.0}$ for oxic bottom waters and a user-defined value ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}<\mathrm{1.0}$ for anoxic bottom waters. The parameter γ_FeS represents the fraction of sulfide that is precipitated as pyrite in the sulfidic zone. The majority of H₂S produced by sulfate reduction is reoxidised, but it is estimated that ∼10–25 % is eventually buried as pyrite (Bottrell and Newton, 2006). However, this fraction can vary significantly over geological timescales (Berner, 1984). If a user does not want to make any assumptions about pyrite precipitation, it can be set to zero (as in the results presented here). The instantaneous equilibrium adsorption coefficients of NH₄ and PO₄ ($K_{{\mathrm{NH}}_{\mathrm{4}}}$, $K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}$, $K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}$) are based on Wang and Van Cappellen (1996) and Slomp et al. (1998), respectively. The first-order rate constants for the sorption of PO₄ to Fe oxides (k_s), the release of PO₄ from Fe-bound P due to Fe-oxide reduction (k_m) and authigenic CFA precipitation (k_a), as well as the pore water equilibrium concentrations for P sorption and CFA precipitation (${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{s}}$, ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{a}}$) and the asymptotic concentration for Fe-bound P (FeP^∞) are taken from Slomp et al. (1996). See Table 10 for a complete summary of the parameters and their values.

(Wang and Van Cappellen, 1996; Slomp et al., 1998)(Slomp et al., 1996)

Table 10 Values for biogeochemical parameters used in OMEN-SED. The variables x, y and z denote the elemental ratio of carbon, nitrogen and phosphorus of the degrading organic matter (here set to C : N : P $=\mathrm{106}:\mathrm{16}:\mathrm{1}$).

3 Stand-alone sensitivity analysis and case studies

3.1 Sensitivity analysis

3.1.1 Methodology

Model parameters implicitly account for processes that are not explicitly resolved. They are notoriously difficult to constrain and thus a primary source of uncertainty for numerical and analytical models, in particular on the global scale and/or in data-poor areas. A comprehensive sensitivity analysis can help quantify this uncertainty and identify the most sensitive parameters. More specifically, sensitivity analysis is used to investigate how the variations in the outputs ($y_{\mathrm{1}},\mathrm{\dots },y_N$) of a model can be attributed to variations in the different input parameters ($x_{\mathrm{1}},\mathrm{\dots },x_M$; Pianosi et al., 2016). Different types of sensitivity indices, which quantify the relative influence of parameter x_i on output y_j with a scalar S_i,j (for $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },M\mathit{\}}$ and $j\in \mathit{\{}\mathrm{1},\mathrm{\dots },N\mathit{\}}$), can be calculated, ranging from simple one-at-a-time methods to statistical evaluations of the output distribution (e.g. variance-based or density-based approaches; Pianosi et al., 2016). The latter indices take values between zero and one ($S_{i,j}\in [\mathrm{0},\mathrm{1}]$), where zero indicates a non-influential parameter and a higher value a more influential parameter. Here, sensitivity analysis is used mainly to identify which parameters have the largest impact on the different model outputs and therefore require more careful calibration. As the probability density functions of our model outputs (i.e. the resulting SWI fluxes) are generally highly skewed towards extreme organic matter degradation rates (not shown), variance-based sensitivity indices may not be a suitable proxy for output uncertainty (Pianosi et al., 2016). Hence, instead the density-based PAWN method by Pianosi and Wagener (2015) is employed which considers the entire conditional and unconditional cumulative distribution function (CDF) of the model output rather than its variance only. The unconditional CDF, F_y(y), of output y is obtained when all uncertain parameters ($x_{\mathrm{1}},\mathrm{\dots },x_M$) are varied simultaneously, and the conditional CDFs, $F_{y|x_i}\left(y\right)$, are obtained when all inputs but the ith parameter are varied (i.e. x_i is fixed to a so-called conditioning value). The sensitivity index of parameter i is measured by the distance between the two CDFs using the Kolmogorov–Smirnov statistic (Kolmogorov, 1933; Smirnov, 1939), i.e.

$$\begin{array}{}\text{(49)}& {\displaystyle }{S}_i=\underset{x_i}{max}\underset{y}{max}\left|F_y\right(y)-F_{y|x_i}(y\left)\right|.\end{array}$$

Since $F_{y|x_i}\left(y\right)$ accounts for what happens when the variability due to x_i is removed, the distance between the two CDFs provides a measure of the effects of x_i on the output y. Due to the model complexity it is impossible to compute the sensitivity indices analytically. Therefore, they are approximated from a Latin hypercube sampling of parameter inputs and calculated outputs. For a brief description of the methodology, see Fig. 4. For more details, we refer the interested reader to Pianosi and Wagener (2015).

Figure 4 (a) Schematic of the PAWN method, plotting an uncertain parameter (x_i) against a generic model output (y). Red dots represent points for calculating the unconditional CDF (NU, here 15), and grey dots are points for calculating each conditional CDF (NC, here 10) with n=2 conditioning points as an example. The user can change the values of NU, NC and n. The number of model evaluations equals N_eval= NU + n⋅ NC ⋅ M, where M is the number of uncertain input parameters. (b, c) Two examples of CDFs for the model-calculated SWI flux of NO₃ using NU = 200, NC = 100 and n=10. The red lines are the unconditional distribution functions F_y(NO₃) and the grey lines are the conditional distribution functions $F_{y|x_i}\left({\mathrm{NO}}_{\mathrm{3}}\right)$ at different fixed values of the input parameters k₁ (b) and $K_{{\mathrm{NH}}_{\mathrm{4}}}$ (c). As the maximal distance between conditional CDFs and unconditional CDF is greater for k₁, this parameter is more influential for the model output (here the SWI flux of NO₃; compare to Fig. 5).

Table 11 Range of model parameters used for the sensitivity analysis of model-predicted output.

Sources: (1) Arndt et al. (2013); (2) Van Cappellen and Wang (1996); (3) Krom and Berner (1980); (4) Gypens et al. (2008); (5) Slomp et al. (1996); (6) Van Cappellen and Berner (1988).

The PAWN method, as implemented within the Sensitivity Analysis for Everyone (SAFE) MATLAB toolbox (Pianosi et al., 2015), is used to investigate M=11 model parameters for ranges as specified in Table 11. Sensitivity indices for all resulting SWI fluxes for two idealised sediment conditions (i.e. anoxic at 400 m and oxic at 4000 m; see Table 12) are calculated. We use NU = 200 samples to estimate the unconditional CDF, NC = 100 samples to estimate the conditional CDFs and n=10 conditioning points. Thus, as ${\mathrm{N}}_{\mathrm{eval}}=\mathrm{200}+\mathrm{100}\cdot \mathrm{10}\cdot \mathrm{11}$, 11 200 model evaluations are performed for each sediment condition. The resulting indices are then translated into a colour code and summarised in a pattern plot to simplify comparison (Fig. 5).

Table 12 Model boundary conditions for the two idealised sediment conditions used for the sensitivity analysis (Figs. 5 and 6). All solute concentrations are in nmol cm⁻³.

3.1.2 Results

Figure 5 summarises the results of the sensitivity analysis as a colour map. Results indicate that generally the most significant parameters for all model outputs are the degradation rate constant for the labile OM pool (k₁) and the fraction of this pool to the total OM stock (f₁). Other parameters play a minor role for the SWI fluxes, with the exception of the secondary redox parameters (i.e. ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}},{\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$) in the oxic scenario. Here, NH₄, SO₄ and H₂S fluxes are very sensitive to changes in ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$ and ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$, as these parameters determine how much of the respective TEA is produced in situ via reoxidation, thus affecting the resulting SWI fluxes. For the oxic scenario, the reoxidation of H₂S produced in the sulfidic layer has also a strong influence on alkalinity (${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$; Table 8, Eq. 5) as it decreases alkalinity by 2 moles per mole of S oxidised (${\mathrm{ALK}}^{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$; Table 10). However, these high sensitivities are partially caused by the wide range of allowed values (${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}\in [\mathrm{0.5};\mathrm{1.0}]$). Yet, for oxic deep sea conditions it is more likely that reduced substances are almost completely reoxidised (e.g. Hensen et al., 2006). For the anoxic scenario the secondary redox parameters are essentially non-influential as no O₂ is available for the reoxidation of reduced substances. Especially for the oxic condition the PO₄ SWI flux appears to be insensitive to P-related parameters (i.e. $K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{ox}}$, $K_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{anox}}$, k_s, k_m, k_a) as the majority is absorbed to Fe oxides. The sensitivities change if other PO₄-related equilibrium concentrations ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{s}}$, ${{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{a}}$ and FeP^∞ are used (not shown). Overall the results of the sensitivity analysis are in line with what one would expect from a diagenetic model and thus provide grounds to confirm that OMEN-SED provides sensible results. The parameterisation of the organic matter pools (f₁) and their degradation rate constants (k₁, k₂) is critical, especially when the model is used in a global Earth system model framework, as these parameters, as well as the γ parameters, can have a very important influence on the flux of dissolved species through the SWI. At the same time these are the weakest constrained parameters. Thus, one should rather choose γ values close to 1 and consider carefully where a relaxation of the “all reoxidised” assumption is appropriate. In contrast, the importance of the OM degradation rate constants cannot be overemphasised. Therefore, much care should be given to how these are parameterised in coupled simulations and a range of different plausible scenarios should be tested to quantify uncertainty.

Figure 5 Pattern plot showing the output sensitivity for each SWI flux (i.e. the chemical compounds on the vertical axis) and each input factor (i.e. the model parameters on the horizontal axis) for two idealised sediment cores. White patterns are assigned where the SWI flux is independent of the specific parameter.

Because of the strong sensitivity of model results to OM degradation rate parameters, we further explore the sensitivity of simulated sediment–water exchange fluxes to variations in organic matter degradation parameters by varying k₁, f₁ and $\stackrel{\mathrm{̃}}{k_{\mathrm{2}}}$ while all other model parameters are set to their default values (Tables 9 and 10). Minimum and maximum values for k₁, $\stackrel{\mathrm{̃}}{k_{\mathrm{2}}}$ and f₁ in the shallow ocean are as in Table 11. For the deep sea condition we account for the presence of more refractory OM by sampling $f_{\mathrm{1}}\in [\mathrm{0.02},\mathrm{0.3}]$, whereas the variation of k₁ and $\stackrel{\mathrm{̃}}{k_{\mathrm{2}}}$ is as in the shallow ocean. The parameter space is sampled using another Latin hypercube approach with sample sizes of N=3500 for each idealised sediment condition. Figure 6 summarises the results of the sensitivity study, and the ranges of observed O₂ and NO₃ sediment–water interface fluxes extracted from a global database (Stolpovsky et al., 2015) are indicated on the colour scale. The colour patterns in Fig. 6a and b reveal the complex interplay between the amount of labile OM f₁ and its degradation rate k₁ for the resulting SWI fluxes of NO₃ in anoxic sediments and O₂ in aerobic sediments. In general, a higher degradation rate in combination with more labile OM available leads to a higher SWI flux. However, higher fluxes extend over a larger range of k₁ values when the amount of labile OM f₁ is high. The absence of a colour pattern in Fig. 6c highlights the limited interaction of the two model parameters for NO₃ SWI fluxes under oxic conditions. Figure 6 shows that SWI fluxes can vary widely over the range of plausible organic matter degradation parameters and that simulated fluxes generally fall within the range of observed SWI fluxes. However, a large number of different k−f combinations can result in SWI fluxes that fall within the observed ranges reported by Stolpovsky et al. (2015), further emphasising the care that should be devoted to constraining OM degradation parameters.

3.2 Case study: simulations of sediment cores

3.2.1 Methodology

In order to illustrate the capabilities of OMEN-SED, comprehensive datasets from the Santa Barbara Basin (Reimers et al., 1996), the Iberian margin and the Nazaré Canyon (Epping et al., 2002) are modelled. Modelled profiles are compared with measured pore water data from different depths including the continental shelf (108 m) and the lower slope (2213 m) located at the Iberian margin, the upper slope (585 m) from the Santa Barbara Basin and a deep sea site (4298 m) in the Nazaré Canyon. The Santa Barbara Basin is characterised by anoxic bottom waters, high POC concentrations and varved sediments (Reimers et al., 1990), and therefore the depth of bioturbation in OMEN-SED is restricted to the upper 0.01 cm. In the uppermost sediments iron(III) hydroxides are reduced, releasing Fe²⁺, which reacts with sulfide to form iron sulfides. Thus, the Fe cycle exerts a strong control on sulfide concentrations in the sediments of this basin (Reimers et al., 1996). In addition, the sediments are generally supersaturated with respect to carbonate fluorapatite by and below 2 cm (Reimers et al., 1996). The Iberian margin, situated in the north-eastern Atlantic, generally belongs to the more productive regions of the global ocean (Longhurst et al., 1995); however, seasonal changes in upwelling creates a strong temporal variability in primary productivity and organic carbon deposition. Submarine canyons in this area (like the Nazaré Canyon) may deliver organic carbon from the shelf to the ocean interior (van Weering et al., 2002; Epping et al., 2002). For a more detailed description of the study areas and the experimental work, the interested reader is referred to the publications by Reimers et al. (1996) and Epping et al. (2002).

In OMEN-SED sediment characteristics and boundary conditions are set to the observed values where available (Table 13). Other sediment characteristics (e.g. sedimentation rate, porosity, density), stoichiometric factors and secondary reaction parameters are set to the default value (see Tables 9 and 10). Organic matter is modelled as two fractions with different first-order degradation rate constants. The POC and pore water profiles were manually fitted by optimising the POC partitioning into the fast- and slow-degrading pool and their respective first-order degradation rate constants (priority is given to reproducing the POC and O₂ profiles). For phosphorus the equilibrium concentration for authigenic P formation (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{a}}$) was adjusted to fit the PO₄ concentration at z_max.

Figure 6 Scatter plots (k₁ vs f₁) of resulting OMEN-SED SWI fluxes for the 400 m anoxic (a: NO₃) and 4000 m oxic (b: O₂, c: NO₃) scenario. Negative values represent a flux from the water column into the sediments. Ranges indicated in red on the colour scale correspond to observed benthic fluxes as reported in the global database of Stolpovsky et al. (2015).

Table 13 Boundary conditions for simulated sediment profiles at the Iberian margin (108 and 2213 m), the Santa Barbara Basin (585 m) and the Nazaré Canyon (4298 m) reported in Fig. 7. For all sites a DIC bottom water concentration of 2400 nmols cm⁻³ is assumed.

3.2.2 Results

Figure 7 Modelled (curves) and measured (filled dots) solid-phase and dissolved pore water profiles for four different sediment cores. Note that different scales are used for different stations. The blue POC curve represents the sum of the refractory (green) and labile (red) POC fraction. The horizontal dashed lines in each panel indicate the bioturbation depth (black) and the penetration depths of oxygen (blue), nitrate (green) and sulfate (red) as calculated by OMEN-SED.

Figure 7 compares modelled and observed sediment profiles for the Santa Barbara Basin and the Iberian margin. Results show that OMEN-SED is able to capture the main diagenetic features across a range of different environments without changing model parameters (other than the four we tuned, i.e. k₁, k₂, f₁ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{a}}$) to site-specific conditions. For the two open Iberian margin stations (108 and 2213 m) OMEN-SED fits all observations well. OMEN-SED does especially well at a sea-floor depth (SFD) of 2213 m by reproducing the deep O₂ penetration and the subsurface maximum in NO₃ concentration due to the nitrification of NH₄ (note that NH₄ is overestimated at this SFD). For the anoxic Santa Barbara Basin (585 m) the decrease in SO₄ and the increase in ALK concentration with sediment depth is well represented, indicating the importance of sulfate reduction as the primary pathway of OM degradation at this site (compare with Meysman et al., 2003). However, a misfit is observed for H₂S and PO₄ in the upper 20 cm of this sediment core. The discrepancy for H₂S can be explained by high iron(III) hydroxide concentrations, which is reduced to degrade organic matter (especially in the 2–4 cm depth interval), therefore placing the beginning of the sulfate reduction zone and the production of H₂S in the deeper sediments (Reimers et al., 1996). Iron processes are currently not dynamically represented in OMEN-SED. In addition, produced dissolved Fe reacts with H₂S to form iron sulfides (e.g. pyrite, FeS₂) and thus further inhibits the rise of H₂S (Reimers et al., 1990). The iron cycle also plays a critical role for phosphorus, as the reduction of iron(III) hydroxides in the surface sediments releases sorbed phosphate, leading to pore waters around and below 2 cm which are supersaturated with respect to fluorapatite, thus initiating CFA precipitation. Reimers et al. (1996) could even show that the accumulation of CFA is mainly restricted to the near-surface sediments (∼5 cm) instead of throughout the sediment column. As OMEN-SED currently does not include an iron cycle and Fe-bound P and CFA processes are highly parameterised, the model is not able to capture these complex, non-steady-state phosphorus dynamics at this specific site. For the Nazaré Canyon station (4298 m) satisfactory fits could be realised apart from NH₄. However, Epping et al. (2002) also could not obtain a better fit using the diagenetic model OMEXDIA. They suggested non-local solute exchange resulting from bio-irrigation as responsible for the higher NH₄ concentrations at this site, which is neglected in their model and in OMEN-SED. Furthermore, the fractured POC profile (indicating episodic depositional events through the canyon) could have been approximated using a different partitioning of the bulk POC into the labile and refractory pool with different degradation rate constants, thus potentially leading to a better fit of the NH₄ profile. In general, better approximations of the data could have potentially been acquired by applying a sensitivity study using different NC ratios (e.g. Epping et al., 2002, report different ratios from Redfield stoichiometry) and exploring the parameter space for the secondary reaction parameters (${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$). However, considering these generalisations and our assumption of steady state, which might not be valid, particularly for the complex Santa Barbara Basin, the shallow core and the Nazaré Canyon, which are affected by seasonality and biology, OMEN-SED generally reproduces the observed pore water trends and hence captures the main diagenetic processes.

Table 14 Sea-floor depth dependency of key model parameters and boundary conditions (adapted from Thullner et al., 2009).

^a Derived from Middelburg et al. (1997).^b Derived from Van Cappellen and Wang (1995). ^c Derived from Conkright et al. (2002). ^d Derived from Boudreau (1997).
^e Calculated with OMEN-SED from POC_flux.

3.3 Case study: stand-alone simulations of global ocean transect

3.3.1 Methodology

In this section we explore to what degree OMEN-SED is capable of capturing the dynamics of organic matter degradation pathways and related TEA fluxes as simulated with a more complete and complex numerical diagenetic model (Thullner et al., 2009). Therefore, we reproduce the simulations of typical conditions along a global ocean hypsometry of Thullner et al. (2009) and compare our modelled TEA fluxes with the results of the complex model and with observations from Middelburg et al. (1996). To explore the global degradation of OM in the sea floor, Thullner et al. (2009) quantified various diagenetic processes using the Biogeochemical Reaction Network Simulator (BRNS; Aguilera et al., 2005), a flexible simulation environment suitable for reactive transport simulations of complex biogeochemical problems (e.g. Jourabchi et al., 2005; Thullner et al., 2005). Thullner et al. (2009) used sea-floor depth (SFD) as the master variable and calculated model parameters, such as w, D_bio and ϕ, from existing empirical relationships (e.g. Van Cappellen and Wang, 1995; Middelburg et al., 1997). Organic matter degradation was described with a 1G approach, thus assuming a single pool of organic matter of uniform reactivity. The first-order rate constant was related to the sediment accumulation rate, w (cm yr⁻¹), following the empirical relationship of Boudreau (1997):

$$\begin{array}{}\text{(50)}& {\displaystyle }k=\mathrm{0.38}\cdot w^{\mathrm{0.59}}.\end{array}$$

This rate constant can be assumed as the mean reactivity of the organic matter fractions which are degraded in the upper, bioturbated 10–20 cm of the sediments. Thus, more reactive fractions (degraded during days or weeks close to the SWI) and more refractory fractions (degraded on longer timescales deeper in the sediments) are not captured by this relationship (Boudreau, 1997). BRNS simulations were performed using boundary conditions and parameters for depths representative of shelf, slope and deep sea sediments (i.e. SFD of 100, 200, 500, 1000, 2000, 3500 and 5000 m). In order to reproduce these results, OMEN-SED is configured here as a 1G model and boundary conditions and model parameters are defined as in Thullner et al. (2009, see Table 14). As OMEN-SED assumes a fixed fraction (i.e. ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}$, ${\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$) of reduced substances to be reoxidised, which exerts a large impact on the resulting SWI fluxes (compare to Sect. 3.1), two sets of simulations are performed in order to show the range of possible model outputs. In the first set-up 95 % of the reduced substances are reoxidised (i.e. ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}={\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\mathrm{0.95}$), and in the second less realistic case, only 5 % are reoxidised (all other model parameters and boundary conditions are equal).

3.3.2 Results

Figure 8 compares simulated SWI fluxes of TEAs (i.e. O₂, NO₃ and SO₄) along the global hypsometry using OMEN-SED (black lines) with the results of Thullner et al. (2009) (red lines). Observations for O₂ and NO₃ fluxes are taken from Middelburg et al. (1996). Due to the applied empirical relations, organic matter flux to the sea floor decreases by 2 orders of magnitude from 100 to 5000 m and its degradation rate constant by 1 order of magnitude (Table 14). Therefore, the rate of organic matter degradation is about 50 times greater at 100 m than at 5000 m (compare to Thullner et al., 2009), thus resulting in a decrease in TEA fluxes along the hypsometry (Fig. 8). The 95 % reoxidation experiments (dots) show proportionally higher O₂ influxes than the 5 % reoxidation experiments (triangles) because more O₂ is utilised for the in situ production of NO₃ and SO₄ in the sediments. This is also mirrored by the increased NO₃ outflux and decreased SO₄ influx for shallower SFDs. This is in line with the results of Thullner et al. (2009), which showed that in situ production is an important pathway of SO₄ supply in the sediment responsible for ∼ 80 % of the total OM degradation at depths between 100 and 2000 m (in our results SO₄ is not used for OM degradation in OMEN-SED below 2000 m). In general, Fig. 8 shows that OMEN-SED captures the main trends in observed and numerically simulated TEA fluxes well. Results also confirm that higher γ values better represent SWI fluxes for most of the global hypsometry. A slight overestimation of shallow ocean SWI fluxes (SFD < 200 m) for the high γ scenario indicates that slightly lower γ values would better capture SWI fluxes in these areas where rapid oxygen consumption favours the escape of reduced species across the SWI.

Figure 8 Fluxes of O₂, NO₃ and SO₄ to the sediment along the global hypsometry. Red lines (with open symbols) are modelled fluxes from Thullner et al. (2009) using BRNS; black lines are results from OMEN-SED (• : ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}={\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\mathrm{0.95}$; ▾: ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}}}={\mathit{\gamma }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\mathrm{0.05}$). Observations of TEA fluxes are taken from Middelburg et al. (1996) (◇: Atlantic, □: Pacific, ×: Arctic–Indian Ocean). Also plotted in the figure (far left panel) are the total oxygen uptake (TOU) estimates of Thullner et al. (2009) (filled red symbols). The green line indicates OMEN-SED results for low oxygen and high nitrate levels and the lower NC ratio. Positive fluxes are directed from the ocean into the sediments.

In addition, observed O₂ fluxes in the upper 2000 m are generally encompassed by our total range in predicted OMEN-SED fluxes. Oxygen fluxes for the deep sea sediments, however, are slightly underestimated. These deviations can presumably be related to the assumed 1G description of organic matter degradation, which neglects the more labile OM pool. This highly reactive pool is degraded close to the sediment surface, thus promoting higher aerobic degradation rates and higher O₂ fluxes. Nitrate fluxes in the upper 500 m of the Atlantic Ocean are well predicted. However, as in Middelburg et al. (1996) the direction of calculated nitrate fluxes in the upper 1000 m of the Pacific Ocean differ from the observations. Middelburg et al. (1996) related these discrepancies to the globally averaged model parameters and the applied boundary conditions. They could reduce the disagreements significantly by using more representative bottom water concentrations for the eastern Pacific and a higher flux of labile organic matter for their 2G model. By changing the boundary conditions and the N : C elemental ratio of organic matter for the whole hypsometry, it is possible to obtain a better model–data fit with OMEN-SED for the shallow Pacific Ocean (green line in Fig. 8b). Bohlen et al. (2012) report that the elemental N : C ratio strongly deviates from Redfield stoichiometry (0.151) with specifically lower values for the eastern Pacific Ocean. The use of their globally averaged value of 0.067 allows the modelled and observed values to be reconciled provided that bottom water conditions are also changed to the low oxygen and high nitrate levels more likely to be found in the shallow Pacific Ocean (O₂=10 nmol cm⁻³ and NO₃=80 nmol cm⁻³).

4 Coupled pre-industrial Earth system model simulations

4.1 Coupling to the cGENIE Earth system model

Figure 9 Schematic of the relationship between OMEN-SED and cGENIE. Arrows represent the information transferred between models.

In a final step, we couple OMEN-SED to the carbon-centric version of the GENIE Earth system model (cGENIE; Ridgwell et al., 2007) in order to illustrate how a fully coupled ocean–sediment system can be configured and applied. We start by providing a brief description of cGENIE and the coupling procedure (Fig. 9).

cGENIE is a model of intermediate complexity based on the efficient climate model C-GOLDSTEIN of Edwards and Marsh (2005), featuring a frictional-geostrophic 3-D ocean circulation model coupled to a fast energy–moisture balance 2-D atmosphere together with a dynamic–thermodynamic sea-ice component. The version of cGENIE used here includes the marine geochemical cycling of carbon, oxygen, phosphorus and sulfur (Ridgwell et al., 2007), the preservation of carbonates in deep sea sediments (SEDGEM; Ridgwell and Hargreaves, 2007) and terrestrial weathering (Colbourn et al., 2013). The ocean model is implemented on a 36 × 36 equal-area horizontal grid with 16 vertical levels using the pre-industrial continental configuration and bathymetry as in Archer et al. (2009). A finer grid (72 × 72) is used for the sediments (see Fig. 11c and  Ridgwell and Hargreaves, 2007) and OMEN-SED is called by SEDGEM for each wet ocean grid point.

In our Earth system model set-up, we prescribe the burial sediment fluxes of detrital material, opal and CaCO₃, while leaving OMEN-SED to calculate organic matter preservation. This assumption serves two purposes. First, the runtime of the model is minimised as steady-state conditions are reached earlier compared to the ca. 20–50 kyr adjustment time for surface sediment CaCO₃ (Ridgwell and Hargreaves, 2007). Second, invariant flux fields remove feedbacks between OMEN-SED and the calculation of CaCO₃ preservation (changes in organic matter preservation affect CaCO₃ dissolution and hence sedimentation accumulation rates, which in turn affects the weight percent of organic matter in the sediments) that would not only lengthen the sediment adjustment time but also make it impossible to carry out unbiased comparisons between different assumptions regarding organic matter reactivity in OMEN-SED.

We derive these fields form the data compilation of Archer (1996) as follows. First, we re-grid the Archer (1996) interpolated non-carbonate mass accumulation rate field (NC_flux) to the 72 × 72 cGENIE sediment grid. This field includes detrital material plus opal (plus a minor contribution from organic matter). We could then directly calculate ∑_flux (total burial flux of all components or total sediment accumulation rate) from this plus measurements of coretop wt % CaCO₃ (C_wtpct) (Archer, 1996) as ${\sum }_{\mathrm{flux}}={\mathrm{NC}}_{\mathrm{flux}}\cdot (\mathrm{1}-\frac{C_{\mathrm{wtpct}}}{\mathrm{100}}{)}^{-\mathrm{1}}$. However, some of the Archer (1996) database C_wtpct values are both close to 100 % and associated with high NC_flux and hence would lead to unrealistically high values for ∑_flux. We therefore impose a plausibility filter by also re-gridding coretop wt % opal (O_wtpct) and quartz (Q_wtpct) and for grid points in which more than one component is reported and the sum exceeds 100 wt %, normalising the individual components (for grid points with only a single solid component, no change is made). We then calculate the individual solid component burial fluxes, and sum them up. To interpolate between the grid points associated with data, we iteratively average nearest (adjoining) neighbours. The distribution of the total burial flux ∑_flux (in g cm⁻² kyr⁻¹) is shown in Fig. C1 in the Appendix.

Depending on the configuration of the overlying biogeochemical ocean model, processes can be included or excluded in OMEN-SED and stoichiometric factors (Table 10) need to be matched between models to ensure preservation of mass. As nitrogen is not modelled explicitly in the employed cGENIE configuration, NC_i, ALK^NIT and ALK^DEN in OMEN-SED are set to zero. cGENIE, however, implicitly includes the effects of NH₄ release and its complete nitrification on alkalinity but neglects the impact of P release. Therefore, alkalinity stoichiometries for aerobic degradation, sulfate reduction and methanogenesis are changed to ${\mathrm{ALK}}^{\mathrm{OX}}=-\mathrm{16}/\mathrm{106}$, ${\mathrm{ALK}}^{\mathrm{SUL}}=\mathrm{122}/\mathrm{106}$ and ${\mathrm{ALK}}^{\mathrm{MET}}=-\mathrm{16}/\mathrm{106}$, respectively (compare to default in Table 10).

Various biogeochemical tracers and parameters are transferred from SEDGEM to OMEN-SED (see Fig. 9) and are converted into the required units. Bottom water concentrations of solutes are converted from mol kg⁻¹ to mol cm⁻³ and the depositional flux of POC (POC_flux) is converted from cm³ cm⁻² yr⁻¹ to mol cm⁻² yr⁻¹ assuming an average density of POC of 1.0 g cm⁻³. Within the water column in cGENIE, POC is partitioned into two fractions with different degradation length scales of ∼ 590 m and 1 000 000 m. The labile pool thus degrades while sinking through the water column, whereas the refractory pool is assumed relatively unreactive (Ridgwell et al., 2007). Thus, depending on sea-floor depth, the partitioning of bulk POC reaching the sediments is different (Fig. 10a, b). This information is used by OMEN-SED to define the parameter f₁. Other parameters used from cGENIE are sea-floor depth and local temperature. The sediment accumulation rate (w) is taken from the previous time step of cGENIE; however, it is assured that w is not smaller than the detrital flux (Det_flux) to the sediments (e.g. w<0 can occur if initially carbonate-rich sediments are eroded during the spin-up of cGENIE). In the case (w≤Det_flux & Det_flux=0.0), all POC is remineralised at the ocean floor. Furthermore, a minimum value of w=0.4 cm kyr⁻¹ is imposed as OMEN-SED tends to be less stable for lower values. For comparison, this threshold is crossed for sea-floor depths below 7000 m when applying the relationship between the sediment accumulation rate and water depth of Middelburg et al. (1997) and below 5200 m for the Burwicz et al. (2011) parameterisation. The bulk POC_flux is separated into the labile and refractory component and the routine to find the steady-state solution for POC is called. Here, the two POC depositional fluxes are first converted into SWI concentrations (POC_i(z=0), in mol cm⁻³) by solving the flux divergence equation

$$\begin{array}{}\text{(51)}& {\displaystyle \frac{\partial F}{\partial z}}=-{\displaystyle \frac{\partial }{\partial z}}\left(-\mathit{\xi }{D}_i{\displaystyle \frac{\partial {\mathrm{POC}}_i}{\partial z}}+\mathit{\xi }w{\mathrm{POC}}_i\right)\end{array}$$

for z=0. OMEN-SED then computes the fraction of POC preserved in the sediment (f_POC, see Eq. 5) and subsequently calls the routines to find the steady-state solutions for the solute substances. Note that in this initial coupling the calculated benthic uptake and return fluxes $F_{C_i}$ of dissolved species C_i (compare to Eq. 6) are adjusted for the advective loss at the lower sediment boundary (w⋅C_i(z_max)) to ensure the conservation of mass in the coupled model:

$$\begin{array}{}\text{(52)}& {\displaystyle }{F}_{C_i}=\mathit{\varphi }\left(\mathrm{0}\right)\left(D_i{\displaystyle \frac{\partial C_i\left(z\right)}{\partial z}}{|}_{z=\mathrm{0}}-w\left[C_i\left(\mathrm{0}\right)-C_i\left(z_{\max }\right)\right]\right).\end{array}$$

In the case that OMEN-SED computes unrealistic results for POC preservation (i.e. f_POC<0.0 or f_POC>1.0) we discard the results of OMEN-SED and all POC is remineralised at the ocean floor. For the modern ocean set-up using the adjustments for w described above, this has not occurred and is just installed as a safety check. Finally, f_POC and the SWI fluxes of solutes ($F_{C_i}$, in mol cm⁻² yr⁻¹) are returned to cGENIE. In the case that no POC is deposited on the sea floor (i.e. POC_flux=0), OMEN-SED is not executed and f_POC and $F_{C_i}$ for all i are set to zero. In order to reduce memory requirements, the sediment profiles (e.g. as shown in Fig. 7) are not calculated in the FORTRAN version of OMEN-SED; however, the boundary conditions are saved and sediment profiles for specific grid cells, ocean basins and ocean transects can be plotted at the end of each experiment using the stand-alone MATLAB version of OMEN-SED.

4.2 Parameterising the OM degradation rate constants in a global model

Figure 10 Idealised relationship of organic matter decomposition during remineralisation in the water column and the sediments. (a, b) Upper panels: water column development of the two organic matter fractions as represented in cGENIE for two ocean depths (red: labile OM with degradation length scale of 589 m; green: refractory OM, which is unreactive in the water column). The values are normalised to OM export at 100 m. Age estimates for the OM since its export from the euphotic zone are calculated using a sinking velocity of 125 m day⁻¹. (a, b) Lower panels: schematic representation of the development of the two OM fractions in the sediments (normalised to OM deposited on the sea floor). For the age estimates in the sediment column a sediment accumulation rate of 0.01 and 0.001 cm yr⁻¹ is assumed. (c) Idealised distribution functions of OM reactive types during remineralisation for different OM ages assuming a reactive continuum model for OM degradation. The initial distribution (at t=0) represents fresh OM when it is exported from the euphotic zone (characterised by $a=\mathrm{3}{e}^{-\mathrm{4}}$ yr⁻¹ and ν=0.125; Boudreau et al., 2008).

As shown in our sensitivity analysis (Sect. 3.1) and discussed by Arndt et al. (2013), the degradation rate constants for OM (k_i) are the most influential parameters and exert a dominant control on the SWI flux of redox-sensitive elements and the preservation of organic matter. Yet, their spatial variability is unknown at the global scale and reported rate constants in the sediments can vary by about 10 orders of magnitude or more (Middelburg et al., 1993; Arndt et al., 2013). Furthermore, when OMEN-SED is coupled to cGENIE, very different timescales have to be considered for OM degradation in the sediments compared to the water column (Fig. 10a, b), and thus the diagenetic rate constants cannot be easily implied by the assumed water column POC flux profiles in cGENIE. To illustrate this, let us consider the degradation of fresh marine organic matter as it is transported and degraded along the ocean–sediment continuum. The bulk material is composed of a complex mixture of different organic carbon compounds that can be described by a reactivity continuum. Microbes preferentially degrade the more reactive organic matter compounds first (Emerson and Hedges, 1988; Wakeham et al., 1997; Lee et al., 2000), resulting in the preferential preservation of more unreactive compounds and rendering the remaining mixture less and less reactive with time. Thus, depending on the age of OM (or depth in the water and sediment column) the reactivity distribution of its compounds changes significantly (Fig. 10c) and the multi-G (2G in this case) approximation of this continuum has to take this shift into account. Figure 10 illustrates these changes in the original reactivity distribution within an ocean–sediment framework. The reactivity distribution t<1 year represents the organic matter mixture after it settled through the water column (Fig. 10c). Only the most reactive OM compounds are remineralised. This explains why the POC flux in the ocean can be represented with a 1G or pseudo 2G degradation model. In the sediments, however, much longer timescales have to be considered and a wider range of more unreactive compounds are degraded. As a consequence, significant changes in the reactivity distribution already take place in the upper millimetres of the sediments (t∼10 years; Fig. 10c). Therefore, a broader range of OM reactive types must be represented by the degradation model to capture the reactivity spectrum of OM in surface sediments, explaining why at least a pseudo 3G model (including two degradable and one refractory fraction; Soetaert et al., 1996; Boudreau, 1997; Stolpovsky et al., 2015) is required. To complicate the situation even further, different sediment depths can represent very different timescales. For instance, half a metre of sediment can be deposited within a year in a coastal setting, while it will represent thousands of years (if not more) of sedimentation in a deep ocean setting. Therefore, residence times and thus degradation rate timescales (or OM age) are mainly controlled by sediment accumulation rates. For instance, assuming a sediment accumulation rate of 0.01 cm yr⁻¹ for the shallow ocean, OM at 5 cm of depth has been degraded for approximately 500 years, whereas a deep ocean sediment accumulation rate of 0.001 cm yr⁻¹ allowed for OM degradation of approximately 5000 years at the same depth. As a consequence, organic matter degradation in deep ocean sediments affects a much wider range of the reactivity continuum and our simple pseudo 3G approximation of the complex OM mixture needs to reflect this by allowing for different k and f values (Fig. 10c). Furthermore, by assuming steady state in OMEN-SED we assume that deposition fluxes of OM are constant over the characteristic timescales of the reaction–transport processes.

Thus, defining appropriate OM degradation rate constants is a major challenge and source of uncertainty for diagenetic models. The rate constants in models are either determined through profile fitting for a specific site or, for global applications, they are related to a single readily available characteristic (or master variable) of the local environmental conditions. For instance, considerable effort has been expended to relate the apparent rate constant for oxic and anoxic OM degradation to sedimentation rate (w) and various empirical relations have been proposed (Toth and Lerman, 1977; Tromp et al., 1995; Boudreau, 1997). Stolpovsky et al. (2015, 2018) suggested empirically derived approaches to constrain degradation rate constants in a 2G model on a global scale. These approaches are derived from present-day observations and might help constrain parameters for present-day applications. However, the problem of constraining 2G degradation model parameters remains for largely different environmental conditions encountered in the past that could also prevail in the future (Arndt et al., 2013). We hence test several alternative schemes in the coupled OMEN–cGENIE framework. Our objective is not to perform and discuss a detailed calibration of the coupled models as this is beyond the scope of this sediment model development paper. Rather, we want to showcase the feasibility of the model coupling, illustrate the range of results and thus the information that can be generated with OMEN-SED, and verify that the model results capture the main observed global benthic biogeochemical features.

4.2.1 Methodology

In this section we compare modelled mean POC weight percentages (wt %) in the upper 5 cm of the sediments (POC_5 cm) to the global distribution pattern of POC content in surface sediments (< 5 cm sediment depth) of Seiter et al. (2004) using different parameterisations for the degradation rate constants k₁ and k₂. For our observational target we take the original POC distribution pattern in $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid resolution (interpolated from > 5500 measurements; compare to Seiter et al., 2004) and transform it onto the 72×72 SEDGEM grid (Fig. 11). The re-gridding of the original POC distribution obviously affects the resolution of the data, especially for the continental margin, as some sites with higher POC wt % are lost in the re-gridding process (compare e.g. maximum values for the eastern Pacific and upwelling waters of the Namibian shelf; Fig. 11a, b). The colour of the points in Figs. 12–13 indicates the sea-floor depth (SFD) of the respective cGENIE grid cell. As the individual data points are highly scattered and in order to see if a certain relation between k₁ and k₂ performs better for specific ocean depths, the data points are binned into six uniform depth classes of 1000 m each (respective mean POC wt % and SFD are represented by the triangles). The regression line is calculated for the six bin classes and included in the figures.

Figure 11 Observed distribution of sediment surface (< 5 cm) POC wt % (a, b) and cGENIE bathymetry (c). (a) Original global distribution of POC wt % interpolated on a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid from more than 5500 individual data points (compare to Seiter et al., 2004, for the interpolation procedure). (b) Observed POC wt % data transformed onto the 72×72 SEDGEM grid. Grid points without any observations are left blank (grey). (c) Gridded continental configuration and ocean bathymetry of the 16-level, 72×72 equal-area cGENIE grid.

To parameterise the reactivity of organic matter in OMEN-SED two different schemes are tested and compared. First, spatially uniform degradation rate constants k₁ and k₂ are assumed. By simulating two different pools of POC in the water column characterised by different degradation length scales (Ridgwell et al., 2007), cGENIE implicitly accounts for the decrease in mean POC reactivity with water depth. The rate constants for the more refractory OM pool, k₂, are systematically varied and the more labile OM component, described by k₁, is assumed to degrade multiple times faster. Values for these parameters are selected in accordance with the best fit to the global surface POC observations (Seiter et al., 2004). Although accounting for the decrease in mean POC reactivity with sea-floor depth, this approach does not take into account the change in the distribution of organic matter reactivity types caused by different sediment accumulation rates and thus different residence timescales in the sediments (Fig. 10). Therefore, the second approach uses the empirical relationship proposed by Boudreau (1997), which relates the apparent OM degradation rate constant in the upper sediments to the sediment accumulation rate, w (cm yr⁻¹; see also Sect. 3.3):

$$\begin{array}{}\text{(53)}& {\displaystyle }{k}_{\mathrm{app}}=\mathrm{0.38}\cdot w^{\mathrm{0.59}}.\end{array}$$

Following Boudreau (1997) and Stolpovsky et al. (2015) it can be assumed that k_app represents the mean OM reactivity within the upper 10–20 cm of the sediments. The following assumptions are made in order to calculate the two degradation rate constants for OMEN-SED:

$$\begin{array}{}\text{(54)}& {\displaystyle }{k}_{\mathrm{app}}& {\displaystyle }=f_{\mathrm{1}}\cdot k_{\mathrm{1}}+(\mathrm{1}-f_{\mathrm{1}})\cdot k_{\mathrm{2}},\text{(55)}& {\displaystyle }{k}_{\mathrm{1}}& {\displaystyle }=x\cdot k_{\mathrm{2}},\end{array}$$

where x describes the relation between k₁ and k₂ and is subject to sensitivity experiments (with values of $x\in \mathit{\{}\mathrm{2},\mathrm{5},\mathrm{8},\mathrm{10},\mathrm{12},\mathrm{15},$ 20,25}). Note that the difference between k₁ and k₂ using this approach is significantly larger as in the globally uniform approach. As the fractions of labile and refractory OM reaching the sediments f₁ is known from cGENIE, k₁ and k₂ can be calculated independently for each grid cell.

To simulate steady-state sediment composition we configure the model as a “closed” system, i.e. one in which there is no loss of CaCO₃ through burial. The redox-dependent P cycle in OMEN-SED is not used in these experiments and all organic phosphorus is returned at the sea floor. To speed up the calculation and to ensure that ocean redox changes caused by OMEN-SED do not impact the sediment composition of CaCO₃, we use the prescribed solid fields as described earlier. Apart from the prescribed fields and the 72×72 sediment grid the model is configured as in Archer et al. (2009) and atmospheric CO₂ is restored to a pre-industrial value of 278 ppmv. First a 20 000 year spin-up is performed without OMEN-SED being coupled. All presented coupled cGENIE–OMEN simulations are run for 10 000 years to steady state from this spin-up. OMEN-SED is called for each grid cell in every time step, feeding back the resulting SWI fluxes and the fraction of POC preserved in the sediments to cGENIE.

4.2.2 Results

Figure 12 Cross plots comparing modelled and observed mean POC wt % in the upper 5 cm of the sediments using the relationship of Boudreau (1997) and the assumptions of Eqs. (54) and (55) to calculate k₁ and k₂. Data points are binned into six uniform depth classes of 1000 m, and each class is represented by a triangle. Grid points with more than 4.0 POC wt % are not shown.

Figure 12 presents results for the relationship of Boudreau (1997) using the assumptions of Eqs. (54) and (55) to calculate k₁ and k₂. Here, the relation between the two degradation rate constants (Eq. 55) is changed globally and thus independent of the sea-floor depth. The cross plots show that it is not possible to achieve a solution in which all bin classes fall onto or close to the 1:1 line. Also, the slope of the regression lines is generally much larger or smaller than 1.0 (with the exception of Fig. 12c), indicating that the relationship between depth and observed POC wt % by bin class is not adequately represented by the model. When looking at the individual bin classes it can be seen that shallow ocean depths are better represented by smaller differences between k₁ and k₂ (e.g. $k_{\mathrm{1}}=\mathrm{5}\cdot k_{\mathrm{2}}$ for SFD < 1000 m; Fig. 12b) and the deep ocean by a larger spread (e.g. $k_{\mathrm{1}}=\mathrm{25}\cdot k_{\mathrm{2}}$ for SFD > 3000 m; Fig. 12h). These results reflect the preferential degradation of more reactive organic matter types (Wakeham et al., 1997; Lee et al., 2000) and thus the change in the distribution functions of OM reactive types for different OM ages (Fig. 10c). In the shallow ocean bulk POC consists of fresher organic matter types on average and is therefore generally more reactive overall (i.e. higher k_app due to higher w in the model) as in the deep ocean. In addition, OM at 5 cm of sediment depth in the deep ocean is generally older than in the shallow ocean due to lower sediment accumulation rates, and therefore more reactive types are affected by degradation and a larger spread between k values is needed to capture these dynamics (compare to Fig. 10c).

Figure 13 Mean POC concentrations in the upper 5 cm of the sediments (POC_5 cm) using the depth-dependent parameterisation $k_{\mathrm{1}}=x\left(\mathrm{SFD}\right)\cdot k_{\mathrm{2}}$ adapted from Boudreau (1997, left) and the globally uniform model (k₂=0.005, $k_{\mathrm{1}}=\mathrm{1.3}\cdot k_{\mathrm{2}}$; right). (a, b) Cross plots as shown in Fig. 12. (c, d) Histograms of the residuals of modelled minus observed POC_5 cm. (e, f) POC_5 cm as calculated with OMEN-SED. (g, h) Difference map of POC_5 cm as calculated with OMEN-SED and interpolated data from Seiter et al. (2004).

Departing from our theoretical considerations (see discussion of Fig. 10c), we use these observations to create a depth-dependent relationship between the two degradation rate constants, where x in Eq. (55) is a function of SFD and takes values of x=5 for SFD < 1000 m, x=8 for 1000 m ≤ SFD < 2000 m, x=12 for 2000 m ≤ SFD < 3000 m and x=25 for SFD ≥ 3000 m for the six SFD bin classes, respectively. Figure 13 compares this depth-dependent approach with the best model using spatially uniform degradation rate constants (k₂=0.005, $k_{\mathrm{1}}=\mathrm{1.3}\cdot k_{\mathrm{2}}$). In the depth-dependent approach all bin classes are close to the 1:1 line and the slope of the regression line (0.9662) indicates that the relationship between depth and observed POC wt % for the bin classes is well represented by the model (Fig. 13a). Using spatially uniform degradation rate constants, five of the six bin classes are located close to the 1:1 line (Fig. 13b), indicating that the simpler parameterisation also adequately captures the relationship between depth and observed POC wt % by bin class. The reason for this is that BIOGEM provides a depth-dependent POC flux and partitioning between the two fractions (Fig. 10). The shallowest bin class (between 0 and 1000 m) represents an exception, as OMEN-SED tends to overestimate POC preservation for this depth class. However, this could also be related to the re-gridding of the original POC distribution pattern of Seiter et al. (2004) onto the SEDGEM grid, as some data grid cells with higher POC wt % on the continental margin are lost due to the restricted SEDGEM resolution (compare to Sect. 4.2). The histograms (Fig. 13c, d) visualise the difference between modelled and observed mean POC concentrations and demonstrate the high density of data points close to the 1:1 line. For the depth-dependent approach, 92.5 % of the cGENIE grid cells show a difference between modelled and observed POC concentration of less than 1.0 POC wt %; in 79.9 % of the grid cells, the difference is less than 0.5 POC wt % (for the globally uniform approach the percentages are 95.37 and 70.95 %).

Both experiments reproduce minimal POC concentrations in the subtropical gyres and generally higher concentrations along the continental margins (Fig. 13e, f). Both experiments, however, underestimate mean POC wt % in the surface sediments of the equatorial eastern Pacific and overestimate POC concentrations in the North Pacific and Southern Ocean (Fig. 13g, h). The depth-dependent approach of Boudreau (1997) shows more spatial variability in POC preservation than the other parameterisation. In general, implementing lower anaerobic degradation rate constants when bottom water oxygen concentrations fall below a threshold value could potentially improve the simulation of higher POC concentrations in areas with high POC input to the sediments (Palastanga et al., 2011).

4.3 Modelled fluxes and sediment characteristics

Figure 14 Sediment characteristics related to POC degradation and oxygen consumption for the depth-dependent parameterisation after Boudreau (1997) with $k_{\mathrm{1}}=x\left(\mathrm{SFD}\right)\cdot k_{\mathrm{2}}$. Total POC_degr rate and fraction of aerobic POC_degr are the respective values for the first 5 cm in the sediments.

For the depth-dependent Boudreau (1997) approach modelled SWI fluxes and sediment characteristics are shown in Fig. 14. Modelled total POC degradation (POC_degr) rates in the upper sediments decrease from the shelves to the deep sea by up to 2 orders of magnitude (Fig. 14b). This is in agreement with data from the literature (e.g. Middelburg et al., 1993, 1997; Burdige, 2007) and other model results (e.g. Thullner et al., 2009) which indicate that the highest degradation rates in marine sediments are found in the coastal ocean (SFD < 200 m). Oxygen fluxes into the sediments (Fig. 14c) range from 0.0 for the deep ocean and sites without OM deposition to values of about 300 µmol cm⁻² yr⁻¹ for the shallow ocean with the highest POC degradation rates. Influx of SO₄ into the sediments is rather low (between 0.0 and 23.9 µmol cm⁻² yr⁻¹) because in OMEN-SED 95 % of produced H₂S is reoxidised to SO₄, and therefore sulfate reduction is mainly driven by in situ sulfide oxidation. However, in general the coupled model fluxes fall well within the ranges predicted by the stand-alone global hypsometry experiments (O₂ between 0.0 and 800 µmol cm⁻² yr⁻¹ and SO₄ between 0.0 and about 300 µmol cm⁻² yr⁻¹; compare to Sect. 3.3). In accordance with the total POC degradation rates the release of PO₄ shows a maximum value of 8.12 nmol cm⁻² yr⁻¹ on the shelves (Fig. 14d). The relative contribution of aerobic POC degradation in the upper sediments increases from the shelves to the deep sea (Fig. 14g), which is also consistent with estimates from Thullner et al. (2009) who found that oxygen is responsible for less than 10 % of POC_degr at 100 m SFD and for more than 80 % in the deep sea. The oxygen penetration depth in OMEN-SED increases from values below 1 cm at the shelves to more than 10 cm in the deep ocean (Figs. 14h and 15). Small oxygen penetration depths of a few millimetres are typical for bioturbated sediments in the coastal ocean (e.g. Wenzhöfer and Glud, 2002) and the oxygen penetration depth has been shown to increase rapidly with SFD to more than 10 cm in the deep sea (Meile and Van Cappellen, 2003; Glud, 2008). Fischer et al. (2009) and D'Hondt et al. (2015) even found cores along a transect in the South Pacific gyre being oxygenated over their entire length (up to 8 m or even 75 m, respectively), which is consistent with our model results (not shown). Simulated mean oxygen penetration depths for the six depth bin classes also agree well with observations compiled by Glud (2008) and Meile and Van Cappellen (2003, Fig. 15).

Figure 15 Sea-floor depth versus O₂ penetration depth for the depth-dependent parameterisation after Boudreau (1997) with $k_{\mathrm{1}}=x\left(\mathrm{SFD}\right)\cdot k_{\mathrm{2}}$. Diamonds represent observations compiled by Meile and Van Cappellen (2003) and squares observations from Glud (2008). Circles are the mean model results for the six SFD bin classes (with standard deviations). Grid cells in which the entire sediment column is oxic (i.e. z_ox=100 cm) are not considered in these statistics (17, 32, 102, 300, 477 and 307 cells for the six bin classes, respectively).

5 Scope of applicability and model limitations

Because of the high computational cost associated with resolving benthic dynamics, most Earth system models of intermediate complexity (EMICs) and also some of the higher-resolution Earth system models either completely neglect or merely include a highly simplified representation of benthic–pelagic exchange processes (Hülse et al., 2017). However, benthic–pelagic coupling plays an important role for carbon cycling and the lack of its representation in EMICs compromises our ability to assess the response and recovery of the Earth system to major past, present and future carbon cycle and climate perturbations. As a consequence, there is a need for benthic biogeochemical models that are able to capture the main features of benthic biogeochemical dynamics, but that are at the same time computationally efficient enough to allow for a direct, dynamic coupling to an ocean biogeochemical model. Therefore, we have developed OMEN-SED, a one-dimensional analytical early diagenetic model that offers a predictive ability similar to complex, numerical diagenetic models at a significantly reduced computational cost.

OMEN-SED is thus not problem specific. Its reaction network resolves the most pertinent benthic biogeochemical species and the most important processes that control their cycling and burial in marine sediments. OMEN-SED can thus be coupled to a wide range of regional to global ocean biogeochemical models, EMICs and higher-resolution Earth system models to investigate a wide range of research questions associated with past, present or future carbon and macronutrient cycling. For instance, OMEN-SED can be used to (i) quantify benthic macronutrient recycling from the shallow coastal to the deep open ocean, (ii) investigate the role of benthic–pelagic coupling in the development of past or future ocean anoxia–euxinia, or to (iii) estimate global organic carbon burial in marine sediments. In theory, its scope of applicability thus ranges from the regional to the global and from the seasonal to the millennial timescale. In order to simulate organic matter preservation in the deeper sediments and thus address questions concerning long-term, geological carbon burial the degradation rate constant for the refractory OM pool has to be scaled down. The resulting larger difference between degradation rate constants can be interpreted as being needed to capture the broader range of OM reactivities degraded over the entire sediment column (see Fig. 10). Instead, more collapsed degradation rate constants are needed to model OM degradation in the upper sediments, such as the first 5 cm as shown in Sect. 4.2.2. In addition, the computational efficiency of OMEN-SED allows for the calculation of quantitative sensitivity indices requiring large sample sizes such as variance- or density-based approaches. Therefore, OMEN-SED can also help quantitatively investigate the sensitivity of benthic model output to systematic variations in model parameters when the model has been tuned to a site-specific problem.

However, OMEN-SED is inevitably associated with a certain degree of simplification that may compromise the applicability of the model in its current version under certain circumstances. First, we have assumed steady-state conditions to allow for an analytical solution to the coupled diagenetic equations. This steady-state assumption is only valid if the variability in boundary conditions and fluxes is generally longer than the characteristic timescales of the reaction–transport processes. As a consequence, OMEN-SED is well suited for coupling to EMICs and the investigation of long-term dynamics in sediment–water exchange fluxes, for instance during past extreme climate events. Yet, in its current version, OMEN-SED is not able to predict the transient response of benthic process rates and fluxes to short-term or seasonal variations in boundary conditions. Future versions of OMEN-SED could approximate non-steady-state conditions by incorporating a time-step-dependent relaxation between different steady states, similar to the schemes used in Ruardij and Van Raaphorst (1995) and Arndt and Regnier (2007). Such a pseudo-transient approach would enable the application of OMEN-SED to systems characterised by high-frequency fluctuations in boundary conditions, such as the coastal ocean or estuaries. Furthermore, by their very nature, analytical models do not allow for overlapping biogeochemical zones or depth-dependent porosity, which introduces a certain error to simulation results. However, the energy-yield-dependent sequence of oxidants is generally valid (e.g. Hensen et al., 2006) and the good agreement between OMEN-SED and the results obtained with a fully formulated numerical RTM (allowing for overlapping TEA use and depth-dependent porosity; Sect. 3.3) shows that these are not critical limitations of OMEN-SED, even for shallow sediments.

Although the model explicitly simulates DIC and alkalinity production and thus has the potential to predict pH profiles within the sediment, a major limitation at this stage is the lack of an explicit description for CaCO₃ precipitation and dissolution coupled to OM decomposition, which also controls the inorganic carbon system (Krumins et al., 2013). In addition, the current version of OMEN-SED does not yet explicitly resolve iron and manganese dynamics (although note that iron is implicitly accounted for in the PO₄ equation). This lack currently limits the applicability of OMEN-SED to iron- and manganese-rich environments, such as coastal marine environments, upwelling regions or ferruginous oceans. In addition, it also compromises the ability of OMEN-SED to predict H₂S fluxes in Fe-rich anoxic environments, where high iron pore water concentrations can deplete pore water H₂S by iron-sulfide mineral precipitation (e.g. Meyers, 2007). Therefore, already planned future extensions of OMEN-SED include an explicit description of carbonate dissolution and iron. Also note that our 1-D diffusion–bioturbation model might not be appropriate to simulate non-accumulating permeable sands of the coastal ocean.

Finally, just as all global models, the global application of OMEN-SED is complicated by the lack of an objective global framework for biogeochemical process parameterisation. The sensitivity study presented here shows that this lack is particularly critical for OM degradation rate parameters (k_i, f_i) and the γ values describing the completeness of secondary redox reactions. A comparison between simulated OM contents and observations indicates that a depth-dependent k−f relationship provides the best fit (Sect. 4.2.2). These results confirm that reducing the continuous distribution of organic matter reactivities into two distinct reactivity classes (2G model) requires different k−f values for shallow vs. deep ocean sediments because of the largely different reaction timescales involved (also see Fig. 10). With respect to γ values, model simulations along the global hypsometry (Sect. 3.3) have shown that high γ values generally capture the main SWI flux features, but have also highlighted that slightly lower γ values would result in a better fit of SWI fluxes to observations of the shallow ocean.

6 Conclusions

Here, we have described in detail and tested OMEN-SED, a new, analytical early diagenetic model resolving organic matter cycling and the associated biogeochemical dynamics. OMEN-SED has been explicitly designed for coupling to EMICs and combines biogeochemical complexity with computational efficiency. It is the first analytical diagenetic model to implicitly represent methanogenesis and explicitly represent oxic degradation, denitrification and sulfate reduction, as well as the reoxidation of reduced substances, adsorption and desorption, and mineral precipitation and dissolution. Explicitly resolved pore water species include O₂, NO₃, NH₄, SO₄, H₂S, DIC and ALK and the solid phase includes two degradable fractions of organic matter, Fe-bound P and authigenic Ca–P minerals.

An extensive sensitive analysis based on the density-based PAWN method (Pianosi and Wagener, 2015) emphasises the importance of OM degradation rate parameters (k_i, f_i) and thus highlights the need for the development of an objective global framework to parameterise OM degradation rate parameters. We have shown that the performance of OMEN-SED at the system scale is similar to that of a fully formulated, multi-component numerical model. The new analytical model is able to reproduce observed pore water profiles across a wide range of depositional environments and captures observed global patterns of SWI fluxes, oxygen penetration depths, biogeochemical reaction rates and surface sediment organic matter contents. Coupled to EMICs or higher-resolution Earth system models, OMEN-SED is thus well suited to examine the role of sediments in global biogeochemical cycles in response to a wide range of past or future carbon cycle or climate perturbations over various timescales.

Appendix A: Reaction network

Table A1 Primary pathways of organic matter degradation, secondary redox reactions and stoichiometries implemented in the reaction network.

Appendix B: Sensitivity analysis

Figure B1 Box plot of parameter sensitivities for the calculated SWI fluxes for the 4000 m oxic condition. Average sensitivities (black lines) and 90 % confidence intervals using N=11 200 model evaluations and Nboot = 100 bootstrap resamples.

Appendix C: Prescribed burial flux fields

Figure C1 Distribution of prescribed total burial fluxes of detrital material, opal and CaCO₃ (in g cm⁻² kyr⁻¹) re-gridded from the data compilation of Archer (1996) using a method explained in the text. Note that latitude and longitude are shown in cGENIE grid cells and not in degrees.