New insights into large-scale trends of apparent organic matter reactivity in marine sediments and patterns of benthic carbon transformation

Constraining the mechanisms controlling organic matter (OM) reactivity and, thus, degradation, preservation, and burial in marine sediments across spatial and temporal scales is key to understanding carbon cycling in the past, present, and future. However, we still lack a detailed quantitative understanding of what controls OM reactivity in marine sediments and, consequently, a general framework that would allow model parametrization in data-poor areas. To fill this gap, we quantify apparent OM reactivity (i.e. OM degradation rate constants) by extracting reactive continuum model (RCM) parameters (a and v, which define the shape and scale of OM reactivity profiles, respectively) from observed benthic organic carbon and sulfate dynamics across 14 contrasting depositional settings distributed over five distinct benthic provinces. We further complement the newly derived parameter set with a compilation of 37 previously published RCM a and v estimates to explore large-scale trends in OM reactivity. Our analysis shows that the large-scale variability in apparent OM reactivity is largely driven by differences in parameter a (10⁻³–10⁷) with a high frequency of values in the range 10⁰–10⁴ years. In contrast, and in broad agreement with previous findings, inversely determined v values fall within a narrow range (0.1–0.2). Results also show that the variability in parameter a and, thus, in apparent OM reactivity is a function of the whole depositional environment, rather than traditionally proposed, single environmental controls (e.g. water depth, sedimentation rate, OM fluxes). Thus, we caution against the simplifying use of a single environmental control for predicting apparent OM reactivity beyond a specific local environmental context (i.e. well-defined geographic scale). Additionally, model results indicate that, while OM fluxes exert a dominant control on depth-integrated OM degradation rates across most depositional environments, apparent OM reactivity becomes a dominant control in depositional environments that receive exceptionally reactive OM. Furthermore, model results show that apparent OM reactivity exerts a key control on the relative significance of OM degradation pathways, the redox zonation of the sediment, and rates of anaerobic oxidation of methane. In summary, our large-scale assessment (i) further supports the notion of apparent OM reactivity as a dynamic ecosystem property, (ii) consolidates the distributions of RCM parameters, and (iii) provides quantitative constraints on how OM reactivity governs benthic biogeochemical cycling and exchange. Therefore, it provides important global constraints on the most plausible range of RCM parameters a and v and largely alleviates the difficulty of determining OM reactivity in RCM by constraining it to only one variable, i.e. the parameter a. It thus represents an important advance for model parameterization in data-poor areas.

1 Introduction

Organic matter (OM) buried in marine sediments represents the largest reactive reservoir of reduced carbon on Earth (e.g. Hedges, 1992). The majority of OM that reaches modern global ocean seafloors originates from contemporary primary productivity (PP) in terrestrial (global net primary productivity, NPP = ∼ 59 Pg C yr⁻¹) or marine environments (global NPP = ∼ 49 Pg C yr⁻¹) (Hedges et al., 1997; Hedges and Oades, 1997; Kharbush et al., 2020). Additionally, it comprises OM derived from chemoautotrophy (Lengger et al., 2019; Vasquez-Cardenas et al., 2020) and ancient PP, the latter delivered through rock weathering, soil erosion, and thawing permafrost (Hedges et al., 1997; Kusch et al., 2021; Grotheer et al., 2020). Since OM production slightly exceeds degradation, a small fraction of the photosynthetically produced organic carbon escapes heterotrophic degradation and is buried in sediments (Berner, 2003; Burdige, 2006). This relatively small imbalance connects the short- and long-term carbon cycles, controls the long-term evolution of atmospheric CO₂, and has enabled the accumulation of oxygen in the atmosphere (Berner, 2003; Hedges, 1992). Nevertheless, several data compilations reveal that rates of OM degradation and preservation vary significantly in space (Arndt et al., 2013; Burdige, 2006; Emerson and Hedges, 1988; Middelburg, 2019) and time (Berner, 2003). Significant progress has been made over the past decades in identifying the potential environmental controls of this spatio-temporal variability (Arndt et al., 2013; Bogus et al., 2012; LaRowe et al., 2020b; Zonneveld et al., 2010). It is becoming increasingly clear that the susceptibility of OM to heterotrophic degradation is not solely an inherent characteristic of the OM itself, but also depends on its environmental context (Kayler et al., 2019; LaRowe et al., 2020b; Mayer, 1995; Schmidt et al., 2011; Zonneveld et al., 2010). In fact, OM reactivity is determined by the dynamic interplay of a plethora of different factors (Arndt et al., 2013; Burdige, 2006; Middelburg, 2019; LaRowe et al., 2020b) (Table 1 and references therein). However, the relative significance of each of these factors remains poorly quantified. Consequently, we currently lack a general framework that allows quantification of the reactivity of OM in a given environmental context. This knowledge gap seriously compromises our ability to constrain OM reactivity in both space and time (Arndt et al., 2013; Lessin et al., 2018). The current lack of such a general framework ultimately limits our ability to quantify pathways and rates of biogeochemical processes over different temporal and spatial scales and, thus, predict the impact and feedbacks of past, present, and future climate change on global biogeochemical cycles and climate (Hülse et al., 2018).

Table 1 Environmental factors influencing organic matter reactivity in marine sediments.

Traditionally, apparent (i.e. estimated representation of) OM reactivity has been determined through laboratory-based incubation experiments in shallow (few centimetres deep) and mixed sediments (e.g. Dai et al., 2009; Grossi et al., 2003; Westrich and Berner, 1984) or integrated model–data approaches in deeply buried sediments (hundreds of centimetres to metres deep and hundreds to thousands and millions of years old) (e.g. Boudreau and Ruddick, 1991; Middelburg, 1989; Wehrmann et al., 2013), using reaction–transport models (RTM). The latter are powerful tools since they allow the extraction of quantitative measures of OM reactivity from comprehensive multi-species porewater datasets integrated over greater sediment depths and ages (e.g. Middelburg, 2019). However, due to our limited quantitative understanding of the environmental controls on OM reactivity, the application of RTM approaches to data-poor settings is still complicated by the difficulty in constraining model parameters.

The need for predictive frameworks that constrain OM degradation model parameters has been recognized since the early days of RTM. Based on the rationale that the reactivity of OM settling onto the sediment is controlled by the degree of degradation, and thus its residence time in the oxygenated water column, early efforts have explored correlations between OM reactivity and water depth, sedimentation rate, or OM fluxes (Boudreau, 1997; Boudreau and Ruddick, 1991; Müller and Suess, 1979; Tromp et al., 1995). However, the derived relationships were based on limited datasets, and a more recent compilation of global data did not reveal statistically significant correlations between OM reactivity and single depositional environmental descriptors (Arndt et al., 2013). However, the latter combined OM degradation model parameters derived from structurally different models and used observational data of different complexity, thus somewhat compromising the comparability of their findings. Nevertheless, it emphasized the role of a complex interplay of environmental controls in reactivity (Kayler et al., 2019; LaRowe et al., 2020b; Mayer, 1995; Schmidt et al., 2011; Zonneveld et al., 2010).

The efforts by Arndt et al. (2013) to develop a general predictive framework highlighted two critical issues: OM reactivity is variable within and across settings with no significative global trends, and there is a crucial need to resolve the lack of comparative studies. Therefore, we here use a common, integrated model–data approach to inversely determine the model parameters that control apparent OM reactivity from contemporary sediment depth profiles (Sect. 3) at 14 different sites from five different depositional environments (Sect. 2). OM degradation is formulated according to the reactive continuum model (RCM) (Boudreau and Ruddick, 1991), which represents the most suitable degradation model type to explore links between reactivity parameters and environmental controls (see Sect. 3.2 for details). Derived RCM parameters are discussed in the environmental context of depositional settings. We further complement the newly derived parameter set with a compilation of previously published RCM parameter values to explore the distribution of model parameters and provide first-order constraints on model parametrization in data-poor areas. Finally, we evaluate the control of OM reactivity on benthic cycling and fluxes by linking RCM parameters to estimated OM degradation rates and relative contributions of metabolic pathways, anaerobic oxidation of methane (AOM) coupled to sulfate reduction, and fluxes of dissolved species across the sediment–water interface (SWI) at each of the 14 studied sites.

Figure 1 Depositional environment geographic locations. (a) Global scale; (b) European margin; (c) Arabian Sea region. Map produced with Ocean Data View (Schlitzer, 2019).

2 Study sites

Our data–model study covers a wide range of depositional environments, from coastal to shelf and slope settings characterized by widely different environmental conditions (Fig. 1, Table 2). It comprises the following benthic provinces (Seiter et al., 2004): northern European margin (EUR1), NW Mediterranean margin off Rhone delta (EUR2), western Arabian Sea (WARAB), Bering Sea (NWPAC), and SW Atlantic off La Plata River (RIOPLATA). Most of those regions are characterized by predominantly temperate climate regimes, except for the Arabian Sea region (arid) and the Bering Sea area (snow/polar) (Chen and Chen, 2013). OM characteristics also cover a broad spectrum, ranging from organic-rich (TOC > 2 wt  %) to lean (TOC < 1 wt  %) sediments (Premuzic et al., 1982; Seiter et al., 2004) and from predominantly terrestrial (coast), to mixed terrestrial and marine (coast/shelf), to predominantly marine (shelf/slope) OM (Hedges et al., 1997; Meyers, 1997). The pelagic ecosystem structure reflects local and regional characteristics. Overall, as a common characteristic they display the occurrence of spring–summer phytoplankton blooms (e.g. Cowie, 2005; Rixen et al., 2019; Fleming-Lehtinen and Laamanen, 2012; Jensen et al., 1990; Coyle et al., 2008). This is also reflected by the composition of the phytoplankton community and thus OM production at the photic zone and vertical transport efficiency to sediments (e.g. Bach et al., 2019). This wide spectrum of depositional environments (Tables 2 and S1) offers a unique opportunity to quantify OM reactivity and explore its multidisciplinary links with environmental drivers.

Table 2 Oceanographic context of depositional environments studied in the large-scale organic matter reactivity quantification.

NA: not available.

3 Methods

Here, we adopt an integrated model–data approach to determine inversely OM reactivity parameters from contemporary total organic carbon (TOC) and sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) profiles measured at 14 sites across the five different depositional environments described in Sect. 2. The following sections provide an overview of the observational data (Sect. 3.1), a brief description of the applied RTM (Sect. 3.2.1; for further details see Supplements, Sect. S2), a detailed description of the inverse model approach used to extract RCM parameters from the observed depth profiles (Sect. 3.2.4), and a summary of the analysed model output (Sect. 3.2.5).

3.1 Observational data

We benefit from previously published datasets to develop our large-scale OM reactivity assessments. TOC and porewater ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles as well as site-specific physical characteristics and supporting data are available for all sites and were compiled from the PANGAEA database (https://www.pangaea.de, last access: 20 November 2020) (except from Rhone delta and Severn estuary datasets), as well as corresponding published literature (see below). Specific details about sampling strategies and analytical methods can be found in the respective publications. The dataset for the Severn estuary (Table 2; row A) is from Thomas (2014). Rhone delta data (Table 2; rows B and C) are from the DICASE project (Rassmann et al., 2016). Data for the northern European margin sites Aarhus Bay, Arkona Basin, and the Skagerrak transect (Table 2; rows D, E, G, H, and I) are part of the Europe-funded project METROL (Aquilina et al., 2010; Dale et al., 2008a, b; Knab et al., 2008; Mogollón et al., 2012). Data from the Helgoland mud area, North Sea (Table 2; row F) (Henkel and Kulkarni, 2015), the Arabian Sea oxygen minimum zone (OMZ) transect (Table 2; rows J, K, and L) (Bohrmann and cruise participants, 2010a, b, c), Bering Sea (Table 2; row M) (Gersonde, 2009), and Argentine Basin (Table 2; row N) (Henkel et al., 2011) are presented in the PANGAEA database (see above). We use those site-specific datasets to inform our data–model analysis and to constrain OM reactivity and related benthic processes.

3.2 Model description

3.2.1 General RTM framework

The Biogeochemical Reaction Network Simulator (BRNS) (Aguilera et al., 2005; Regnier et al., 2002) is an adaptive simulation environment that has been successfully employed to reproduce and quantify diagenetic processes in marine sediments across a wide range of depositional environments and timescales (Dale et al., 2008a; Thullner et al., 2009; Wehrmann et al., 2013). The BRNS is suitable for large, mixed kinetic–equilibrium reaction networks (Dale et al., 2009; Jourabchi, 2005; Thullner et al., 2009). The concentration depth profiles of solid and dissolved species in marine sediments are calculated according to the vertically resolved mass conservation equation of solid and dissolved species in porous media (Berner, 1980; Boudreau, 1997):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\sigma }{C}_i}{\partial t}}& ={\displaystyle \frac{\partial }{\partial z}}\left(D_{\mathrm{bio}}\mathit{\sigma }{\displaystyle \frac{\partial C_i}{\partial z}}+D_i\mathit{\sigma }{\displaystyle \frac{\partial C_i}{\partial z}}\right)-{\displaystyle \frac{\partial \mathit{\sigma }\mathit{\omega }{C}_i}{\partial z}}\\ & +{\mathit{\alpha }}_i\mathit{\sigma }\left(C_i\left(\mathrm{0}\right)-C_i\right)+{\sum }_n{s}_i^n{r}_n.\end{array}\end{array}$$

The first three terms on the right-hand side represent the transport process (bioturbation and molecular diffusion, advection, and bioirrigation), whereas the last term denotes the sum of all reactions (production and consumption) affecting species i. Table 2 summarizes all symbols employed here.

Briefly, the implemented reaction network encompasses the most pertinent primary and secondary redox reactions found in the upper layers of marine sediments (e.g. Aguilera et al., 2005; Thullner et al., 2009; Van Cappellen and Wang, 1996; Wang and Van Cappellen, 1996). It explicitly accounts for the heterotrophic degradation of OM coupled to the consumption of oxygen (aerobic OM degradation), nitrate (denitrification), sulfate (organoclastic sulfate reduction), and methanogenesis. Additionally, the reaction network accounts for secondary redox reactions, i.e. re-oxidation of reduced species produced during primary redox reactions. It explicitly resolves nitrification, sulfide re-oxidation by O₂, anaerobic oxidation of methane (AOM) coupled to sulfate reduction, and CH₄ reoxidation by O₂.

The model transport and reaction equations were solved sequentially according to Regnier et al. (1998). The size of the model domain was fixed at 1000 cm for all sites, except for the Bering Sea, in which the model domain was extended to 1500 cm, due to the low sedimentation rate assumed for this site (Table 2). This choice is based on initial tests and ensures that the model domain covers the diagenetically most active zone, thus reducing the influence of biogeochemical dynamics in underlying sediments on biogeochemical dynamics within the model domain. The model equations were solved on an uneven grid and with a time-adapting time step (Δt < 0.01) run until steady state was reached. The Supplement includes a detailed description of the general RTM framework used here, its parametrization, and its numerical solution (Sect. S2).

3.2.2 OM degradation model

OM is composed of a complex and dynamic mixture of compounds that are distributed over a wide, continuous spectrum of reactivities. Thus, OM degradation is described by the reactive continuum model (RCM) (Boudreau and Ruddick, 1991), which assumes a continuous distribution of OM compounds over the entire reactivity spectrum. The RCM assumes a gamma distribution to describe the probability density function of OM distribution, om(k,t) (Aris, 1968; Boudreau and Ruddick, 1991; Ho and Aris, 1987). The overall rate of OM degradation R_TOC is given by the integral

$$\begin{array}{}\text{(2)}& r_n=R_{\mathrm{TOC}}=-{\int }_{\mathrm{0}}^{\mathrm{\infty }}k\cdot om\left(k,t\right)\mathrm{d}k.\end{array}$$

Here, om(k,t) determines the concentration of OM having a degradability between k and k+dk at time t, with k being equivalent to the first-order degradation rate constant. The initial distribution (t=0) of OM is given by

$$\begin{array}{}\text{(3)}& om(k,\mathrm{0})={\displaystyle \frac{\mathrm{TOC}\left(\mathrm{0}\right)\cdot a^v\cdot k^{v-\mathrm{1}}\cdot e^{-a\cdot k}}{\mathrm{\Gamma }\left(v\right)}},\end{array}$$

where TOC(0) denotes the initial OM content, Γ(v) is the gamma distribution, a is the average lifetime of the more reactive components of bulk OM, and v represents the dimensionless scaling parameter of the distribution near k=0. The free, positive parameters a and v delineate the shape of the initial distribution of OM compounds along the range of k, and thus the overall reactivity of bulk OM. As such, the RCM approach requires the definition of two parameters that will define the shape of the OM distribution over reactivity k (Sect. 3.2.5, Eq. 5).

Due to the rapid depletion of the most reactive compounds, the reactivity of the bulk material decreases during degradation, reflecting the widely observed reactivity decrease with burial time, depth, and age (Boudreau and Ruddick, 1991; Middelburg, 1989). This indicates that degradation of OM proceeds at different rates in parallel. Interactions between different compounds or transformations of compounds can change the reactivity of a given compound. While the RCM does not explicitly account for such interactions and transformations, the overall OM profiles take these interactions implicitly into account. Thus, OM compounds are continuously and dynamically distributed over a range of reactivities that capture the decrease in apparent reactivity with burial age and depth as the most reactive compounds are successively degraded. The interplay of a and v drives the bulk OM reactivity depth profiles and consequently the yielded profiles from the reaction network implemented in the BRNS. The evolution of OM concentration as a function of depth z is given by

$$\begin{array}{}\text{(4)}& \mathrm{TOC}\left(z\right)=\mathrm{TOC}\left(\mathrm{0}\right)\cdot {\left({\displaystyle \frac{a}{a+\mathrm{age}\left(z\right)}}\right)}^v,\end{array}$$

where a and v are the free RCM parameters (Arndt et al., 2013; Boudreau and Ruddick, 1991) and age(z) denotes the age of the sediment layer at depth z. For non-bioturbated sediments (z> z_bio) the burial age(z) can be calculated as a function of the burial velocity (Sect. S2.3, Eq. S8). However, within the bioturbated upper sediment layers, the age distribution of reactive species is controlled by both sedimentation, bioturbation, and the reactivity k of reactive species (Meile and Van Cappellen, 2005). In such cases, we apply a multi-G approximation for the RCM in the bioturbated sediments (Sect. S2.3, Eqs. S9–S18).

3.2.3 Boundary conditions

Boundary conditions place the BRNS in the environmental context of each of the study sites (Table 4). The boundary conditions are constrained based on either site-specific measurements or alternatively published data if direct observations were not available. Here, we assume a fixed boundary concentration of OM at the SWI. Concentrations of O₂ are based on measurements in Skagerrak (Canfield et al., 1993) and the Rhone delta (Rassmann et al., 2016) sediments. Similarly, concentrations of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ are based on measurements in Skagerrak sediments (Canfield et al., 1993). Because of the lack of site-specific information, these O₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ measurements are adopted at all oxic sites (where ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\mathit{>}\mathrm{28}$ mM) since they are representative for coastal and shallow slope sites, and their exact value exerts no influence on deep ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and CH₄ concentrations. For anoxic sites (where ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ < 28 mM), O₂ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ bottom water concentrations are set to zero. ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations were constrained based on measurements for each simulated site. Concentrations of CH₄ and HS⁻ were set to zero at SWI, as those species are rapidly oxidized in the overlying bottom water. A no-flux $\left(\frac{\partial C}{\partial z}=\mathrm{0}\right)$ condition was implemented for all species at the lower boundary, assuming that biogeochemical dynamics in underlying sediments exert no influence on diagenetic processes in the model domain (e.g. Thullner et al., 2009).

3.2.4 Inverse modelling

We use an inverse model approach to extract the optimal OM reactivity parameter set a and v (i.e. the apparent OM reactivity) by assuming that the rank of a parameter set depends on the similarity between simulated and measured data. In general, all inverse modelling approaches are ultimately affected by a few limitations. Among others, the quality of the inverse modelling results is directly affected by the quality of the observational data. In particular, the lack of observations and data at the SWI due to core top loss during sediment sampling by gravity and piston corers can impact results. Additionally, processes that are not explicitly described in the applied model can impact inverse model results because reaction rate parameters implicitly account for those processes. Furthermore, downcore changes in concentration–depth profiles may also be influenced by transient features in bottom water conditions and/or fluxes, and the validity of the steady-state assumption thus needs to be carefully evaluated. Finally, multiple parameter sets might fit the observations equally well. Because both a and v exert an influence on the apparent OM reactivity and its evolution with depth (see below, Eq. 5), they are not completely independent parameters. For instance, a decrease in v (decreasing reactivity) can be compensated for by a simultaneous decrease in parameter a (increasing reactivity). Consequently, different a and v parameter couples might result in equally statistically satisfying fits between the simulated and observed profiles.

These limitations can be alleviated by using comprehensive, multi-component observational datasets, because the inclusion of additional information adds further constraints. However, such datasets are often not available. For instance, the first efforts to quantify a and v parameter values, and thus reactivity, on the basis of observations solely focused on downcore TOC profiles (Boudreau and Ruddick, 1991). Such an approach was instrumental for the development of the RCM and for producing the first large-scale dataset of a and v. However, the increasing availability and accessibility of multi-species datasets allow for the use of additional information that is available across different sites. Because OM heterotrophic degradation is the driver behind biogeochemical dynamics in marine sediments, downcore profiles of TEAs (e.g. ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) (e.g. Bowles et al., 2014; Jørgensen et al., 2019b) or reduced products of OM degradation (such as CH₄ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$) could be used as further constraints. However, data availability is often still limited; therefore, it is important to identify a minimum set of observational data that are widely available and comparably easy to measure. We tested the suitability of different artificial porewater datasets for the inverse determination of apparent OM reactivity. Using an artificially generated dataset (see Supplement Sect. S3), we ran the BRNS in the same fashion as we did for the site-specific data–model approach (see below). Our sensitivity analysis confirms that when TOC is considered as a single constraint, multiple pairs of a and v produce equivalent TOC depth profiles (Supplement Fig. S1). Model results show that including ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles as an additional constraint (Supplement Fig. S2) facilitates the identification of a best-fit a and v parameter set (Supplement Fig. S3). CH₄ profiles were used as an additional qualitative constraint when the data were available, although measured ex situ CH₄ profiles are often not reliable due to sampling and measurement uncertainties in deeper sediment layers where concentrations exceed millimolar levels (Dale et al., 2008a; Hensen et al., 2003; Hilligsøe et al., 2018). The dynamics of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and CH₄ are solely controlled by OM degradation and AOM (e.g. Bowles et al., 2014; Egger et al., 2018; Regnier et al., 2011), although their distributions can also be affected by bioirrigation in particularly shallow sulfate–methane transition zones (SMTZs) (e.g. Dale et al., 2019). In anoxic settings and deeply buried sediments, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is the dominant TEA and CH₄ is the most common reduced species (Jørgensen et al., 2019a, b). Thus, a combination of TOC, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, and CH₄ (if available to verify the depths of the SMTZ) depth profiles incorporates the information contained in the observed benthic sulfur and carbon dynamics and is sufficient to extract robust estimates of apparent OM reactivity and its evolution from the sediment–water interface down to the SMTZ. In addition, because changes in a–v exert different effects on TOC and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles (see Figs. S1 and S2), including these two species reduces the impact of the a–v correlation on the uniqueness of fit. Thus, here we perform a site-specific data–model fit based on TOC and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (and CH₄) depth profiles.

The optimal parameter set was then determined for each site by assuming that the rank of a parameter set depends on the similarity between simulated and measured data. Best-fit parameter couples were found by minimizing the misfit between model results and simulations. We quantified the similarity between the simulated and observed TOC and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles in terms of correlation coefficient (r), centred root-mean-square difference (RMSD), and standard deviations (SD). All these measures can be graphically represented in a single Taylor diagram (Taylor, 2001) that statistically quantifies the misfit between the observed and simulated data and visualizes how closely the simulated depth profile resembles the observed one.

The best-fit RCM parameters were inversely determined by first running the model for each site with a set of a–v couples over the entire range of previously published a ($a={\mathrm{10}}^{-\mathrm{3}}$–10⁷ year) and v ($v={\mathrm{10}}^{-\mathrm{2}}$–10⁰) values (Arndt et al., 2013, and references therein). The approach adopted for the sampling of the parameter space is a critical step in inverse modelling since the chosen statistical measures may reveal, in addition to a global optimum, multiple local optima. Because of the comparably low computational cost of the forward modelling, the low number of unknown parameters (here, a and v), and the weak nonlinearity of the problem, we used a simple two-step nested regular sampling of the two-dimensional parameter space. First, the entire plausible a–v space (see above) was sampled on a coarse regular grid (Δ log(a) = 1, Δv=0.02), and, for each of these runs, r, RMSD, and SD were calculated for both TOC and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ profiles. Based on these measures, a new, more finely resolved a–v grid (10×10) was defined around the a–v couples that showed the smallest combined misfit between simulated and observed TOC and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ profiles. On this grid, a–v couples were again sampled regularly, and the optimized a–v couple was determined by the lowest combined misfit (Table 5).

3.2.5 Quantification of organic matter reactivity and associated benthic processes

The apparent reactivity of bulk OM k for each site is calculated based on the RCM reactivity parameters a and v (Boudreau and Ruddick, 1991) that yield the best data–model fit and is given by

$$\begin{array}{}\text{(5)}& k\left(z\right)={\displaystyle \frac{v}{\left(a+\mathrm{age}\left(z\right)\right)}}.\end{array}$$

The apparent reactivity of bulk OM at the sediment–water interface (z = 0) is thus given by

$$\begin{array}{}\text{(6)}& k\left(\mathrm{0}\right)={\displaystyle \frac{v}{a}}.\end{array}$$

The BRNS simulates downcore concentration profiles for each species i that is explicitly resolved in the reaction network, as well as reaction rates r_n (Table S4). From those reaction rates, the total rate of organic carbon oxidation, i.e. the depth integral of OM degradation for each primary redox reaction n, can be calculated by integrating the rate over the entire model domain L (e.g. Thullner et al., 2009):

$$\begin{array}{}\text{(7)}& R_{\mathrm{TOC},\mathrm{total}}={\sum }_{n=\mathrm{1}-\mathrm{4}}{\int }_{\mathrm{0}}^L{r}_n\mathrm{d}x=-{\int }_{\mathrm{0}}^{\mathrm{\infty }}k\cdot om\left(k,t\right)\mathrm{d}k,\end{array}$$

where n=1, 2, 3, and 4 denote aerobic OM degradation, denitrification, organoclastic sulfate reduction, and methanogenesis, respectively (Table S4). The OM fluxes (J_in) at the SWI are calculated as the sum of OM out of the sediment (J_OM,out) and the burial fluxes at the base of the model domain (J_OM,bur) (e.g. Burdige, 2006). Here, we assume that J_OM,out is equivalent to the depth-integrated rates of OM degradation (Eq. 7) and J_OM,bur is solely governed by advective processes at depth z=1000 cm. As such

$$\begin{array}{}\text{(8)}& {\displaystyle }{J}_{\mathrm{OM},\mathrm{out}}={\sum }_{n=\mathrm{1}-\mathrm{4}}{\int }_{\mathrm{0}}^L{r}_n\mathrm{d}x,\text{(9)}& {\displaystyle }{J}_{\mathrm{OM},\mathrm{bur}}=\left(\mathrm{1}-{\mathit{\phi }}_z\right)\cdot {\mathit{\omega }}_z\cdot {\mathrm{TOC}}_z,\text{(10)}& {\displaystyle }{J}_{\mathrm{OM},\mathrm{in}}=J_{\mathrm{OM},\mathrm{out}}+J_{\mathrm{OM},\mathrm{bur}}.\end{array}$$

Additionally, the fluxes of species i at the SWI, i.e. the benthic–pelagic fluxes, are calculated from the model-derived depth profiles of species i according to e.g. Burdige (2006):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{J}_i& =J_{i,\mathrm{Advection}}+J_{i,\mathrm{Diffusion}}\\ & +J_{i,\mathrm{Bioturbation}}+J_{i,\mathrm{Bioirrigation}}.\end{array}\end{array}$$

Positive J_i values indicate upward fluxes from the sediment, whereas negative values represent downward fluxes into the sediment.

4 Results and discussion

The integrated data–model analysis yielded a comprehensive picture of OM reactivity parameters (a, v, and k(z)), as well as OM degradation dynamics and benthic–pelagic coupling for our large-scale compilation of depositional environments (Table 6). Because model parameters implicitly account for all the processes that are not explicitly described in the model, this variability provides important insights into the environmental controls on OM reactivity.

Figure 2 Total organic carbon data–model best fit: (a) Severn estuary; (b) Rhone pro-delta; (c) Rhone shelf; (d) Aarhus Bay; (e) Arkona Basin; (f) Helgoland mud area, North Sea; (g) Skagerrak – S10; (h) Skagerrak – S11; (i) Skagerrak – S13; (j) Arabian Sea – OMZ core; (k) Arabian Sea – OMZ transition; (l) Arabian Sea – below OMZ; (m) Bering Sea; (n) Argentine Basin.

4.1 Constraining apparent OM reactivity

4.1.1 Inverse modelling

We inversely determined the set of RCM parameters a and v that produces the best model fit to the observed TOC (Fig. 2) and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Fig. 3) depth profiles. Table 5 shows the quantitative assessments of these model fits and associated uncertainties. Overall, TOC profiles display a satisfying fit between model and observations. However, the quality of fit is often compromised by a lack of observations and data at the SWI due to core loss during sediment sampling by gravity and piston corers, poor downcore sampling resolution, and possible non-steady-state conditions. Generally, the latter may be problematic when trying to reproduce sedimentary conditions of dynamic shallow systems and long-term burial in shelf and deep-sea sites. At estuarine sites, steady-state conditions may be impacted by the highly dynamic nature of such settings. The Severn estuary exhibits strong tidal dynamics that result in continuous cycles of sediment resuspension (Manning et al., 2010). The Rhone delta experiences strong pulses of fresh water and sediments associated with flood events, as well as variable sedimentation rates (Antonelli et al., 2008; Cathalot et al., 2010; Zebracki et al., 2015). In shelf and slope settings, steady-state assumptions can be impacted by fluctuations on sedimentation rates and lateral transport of sediments across shelves (e.g. Henkel et al., 2011; Riedinger et al., 2014, 2017; Hensen et al., 2003; Cowie, 2005). Further, due to the possible loss of the upper few millimetres to centimetres during sampling with gravity and piston cores, it can be challenging to quantify OM contents at the SWI. Some of our sites (e.g. northern European sites; Aquilina et al., 2010; Knab et al., 2008) are likely to be impacted by such factors. Overall, such limitations in inverse model results are highlighted by the somewhat low correlation coefficients for TOC profiles in most cases (r_TOC < 0.75; Table 3). We find that adding additional constraints to the model–data fit relieves such limitation. The ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depth profiles (Fig. 3) display a good agreement between present-day observations and model simulations, with the majority of sites exhibiting high correlation coefficients ($r_{{\mathrm{SO}}_{\mathrm{4}}}\mathit{>}\mathrm{0.95}$; Table 3). Nevertheless, Fig. 3 shows cases where the model-derived sulfate depth profiles display some mismatches relative to observations (Arabian Sea below OMZ, Fig. 3l; Bering Sea, Fig. 3m). At the Bering Sea, the data–model discrepancy is likely associated with kinks in both TOC (Fig. 2m) and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (Fig. 3m) profiles. They likely represent changes in the depositional regime over the studied timescale (e.g. Hensen et al., 2003) that are not captured by the model. The mismatch in the Arabian Sea at the site below the OMZ (Fig. 3l) could originate from ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ recycling through H₂S oxidation by MnO₂ and Fe(OH)₃ (e.g. Rassmann et al., 2020; Schulz et al., 1994); however we are unable to further explore this hypothesis due to a lack of metal oxide data to constrain those biogeochemical reactions.

Figure 3 Sulfate data–model best fit: (a) Severn estuary; (b) Rhone pro-delta; (c) Rhone shelf; (d) Aarhus Bay; (e) Arkona Basin; (f) Helgoland mud area, North Sea; (g) Skagerrak – S10; (h) Skagerrak – S11; (i) Skagerrak – S13; (j) Arabian Sea – OMZ core; (k) Arabian Sea – OMZ transition; (l) Arabian Sea – below OMZ; (m) Bering Sea; (n) Argentine Basin.

Table 3 Summary of model elements incorporated in the BRNS.

4.1.2 Apparent OM reactivity parameters v and a

The inversely determined OM reactivity parameters a and v (Sect. 3.2.4) derived from the best-fit model solutions (Table 6) reveal important differences in apparent OM reactivity and its distribution on a global scale. Alongside a and v parameters compiled from the literature (Arndt et al., 2009; Boudreau and Ruddick, 1991; Contreras et al., 2013; Henkel et al., 2012; Mewes et al., 2016; Mogollón et al., 2012, 2016; Wehrmann et al., 2013; Ye et al., 2016), they provide useful constraints on the plausible range of a and v, as well as their control on the global variability in apparent OM reactivity.

Parameter v exerts a scaling influence on OM reactivity distribution. High v values result in a high initial OM reactivity, whereas low v values produce low initial reactivity (Arndt et al., 2013). For instance, an extremely low v value of ≤ 0.01 for OM deposited during the late Cretaceous Oceanic Anoxic Event 2 at Demerara Rise (Arndt et al., 2009) has been linked to intense and rapid sulfurization of OM in the strongly euxinic subtropical Atlantic Ocean (Hülse et al., 2019). While previously published values for contemporary ocean sediments fall within the range of 0.1 to 1.08 (Arndt et al., 2013; Boudreau and Ruddick, 1991), inversely determined v values predominantly fall within the range of 0.1–0.2 (apparent 6th–11th order of reaction) – although lower and higher v values were also determined (Fig. 4a). Newly determined v parameters thus confirm the previously observed dominance of v values between 0.1 and 0.2 (Boudreau and Ruddick, 1991) and suggest that v is relatively constant across a wide range of different depositional environments encountered in the contemporary ocean. These findings also agree with the parametrization of another empirically derived continuum model description of OM degradation (Middelburg, 1989), which applies a globally constant parameter of 0.16:

$$\begin{array}{}\text{(12)}& k\left(z\right)=\mathrm{0.16}\cdot (a+\mathrm{age}(z){)}^{-\mathrm{0.95}}.\end{array}$$

Although the parameters of the so-called power model cannot be directly compared to those of the RCM due to different exponents, the RCM is nevertheless equivalent to a general power law with an exponent of −1:

$$\begin{array}{}\text{(13)}& k\left(z\right)=v\cdot (a+\mathrm{age}(z){)}^{-\mathrm{1}}.\end{array}$$

Thus, the v parameter range of the RCM is in good agreement with the 0.16 value imposed in the power-law model.

Figure 4 Frequency histogram of organic matter reactivity parameter distribution on a global scale: (a) scaling parameter – v; (b) shaping parameter – a.

Table 4 Site-specific transport parameters and upper boundary conditions implemented in the reaction–transport model.

Nevertheless, our results also reveal exceptions. Particularly high v values (v > 0.25) were determined for the Arabian Sea and Bering Sea. Similarly, Boudreau and Ruddick (1991) also report exceptionally high values (v≈1) for the Peru Margin and north Philippine Sea but invoked the possibility of non-steady-state depositional conditions as a possible explanation. In our study, v>0.25 also coincides in part with the poorest data–model fits, such that the robustness of such high v values is uncertain (Table 5). In particular, the high v value in the Bering Sea (Table 6) might result from uncertainties (see above, Sect. 4.1.1) associated with boundary conditions and transient features (e.g. Henkel et al., 2011; Hensen et al., 2003). However, the model-derived v > 0.25 for the Arabian Sea agrees with previous diagenetic model results for that region. Luff et al. (2000) developed a 3G model to investigate OM dynamics and found that bulk OM contains exceptionally reactive compounds with first-order k of 0.2–30 yr⁻¹. In the Arabian Sea, the OM flux to the sediment is supported by high primary production rates during the monsoon season (e.g. Cowie, 2005), complemented by deep carbon fixation by microorganisms (Lengger et al., 2019). Observations indicate an efficient vertical transport of diatom and haptophyte detritus to the sediment (Barlow et al., 1999; Cowie, 2005; Latasa and Bidigare, 1998; Shalapyonok et al., 2001). High transport rates and pronounced oxygen minimum zones might limit heterotrophic degradation of the sinking OM in the water column, thus delivering comparably unaltered OM to the sediment. Biomarker analysis of benthic OM indicates the deposition of highly reactive OM (Koho et al., 2013) and a reactivity decrease with extended exposure to oxic conditions in the water column (Vandewiele et al., 2009). Therefore, the somewhat atypical high v values determined (together with comparably low a values; see below) for the Arabian Sea OMZ transect are supported by interdisciplinary observations. Thus, v has a relatively narrow range of 0.1 to 0.2, with higher > 0.2 values generally restricted to depositional environments that are characterized by high primary production rates, anoxic water conditions, and enhanced vertical OM fluxes within the OMZ.

Parameter a determines both the overall apparent reactivity and the shape of the reactivity decrease with depth. High a values result in an overall low OM reactivity, but a slow decrease in reactivity with depth. In contrast, low a values result in a higher reactivity but also a more rapid decrease in reactivity with sediment depth (Arndt et al., 2013). Parameter a is often related to the degree of “pre-ageing” or degradation (e.g. Middelburg, 1989; Middelburg et al., 1993). Very fresh marine OM has been shown to be characterized by substantially low a values ($a=\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{4}}$ years) (Boudreau and Ruddick, 1991), while the heavily sulfurized OM deposited during Cretaceous Oceanic Anoxic Event 2 in the strongly euxinic subtropical Atlantic Ocean is characterized by very high a values (a=10⁴ years) (Arndt et al., 2009). Our inversely derived a values span several orders of magnitude (Table 6; Fig. 4b). They thus corroborate the high variability observed in previously published a-value compilations (e.g. Arndt et al., 2013; Boudreau and Ruddick, 1991) and suggest that the variability of parameter a exerts the major control on the spatial variability in OM reactivity. Additionally, inversely determined a values also broadly agree with previously observed global patterns. Low values (a < 10² years; i.e. high apparent OM reactivity) are generally determined for depositional environments that are characterized by the rapid deposition of predominantly fresh marine OM (e.g. Arabian Sea and shallow sites at the northern European margin), while depositional settings that receive pre-degraded OM and/or complex mixtures of OM (e.g. Skagerrak) exhibit high a (a > 10³ years). However, despite the broad range of a values, the ensemble of inversely determined a values shows that most values fall into the range of 10⁰−10⁴ years, thus narrowing the range of a values.

Table 5 Summary of statistical tests for the best-fit a–v pair at each site derived from Taylor diagrams.

4.1.3 OM reactivity distributions

Inversely determined parameters a and v also allow the assessment of the apparent OM reactivity evolution during burial. Here, we calculate the probability distribution of OM over k at the time of deposition on the SWI, as well as after burial at 10, 100, and 1000 cm b.s.f. (centimetres below seafloor). Since v values are broadly similar across depositional settings, the differences in k distributions are mainly driven by the range of a values (Fig. 5). However, it is important to note that sedimentation rates ω and OM contents also influence the evolution of these distributions as they co-control degradation rates (Eqs. 5, S8). Therefore, the changes in OM−reactivity distributions at the considered burial depths are not entirely driven by differences in reactivity, but more broadly reflect different exposure to heterotrophic degradation across sites.

Figure 5 Apparent organic matter reactivity k initial probability distributions and downcore evolution as a function of burial depths: (a) Severn estuary; (b) Rhone pro-delta; (c) Rhone shelf; (d) Aarhus Bay; (e) Arkona Basin; (f) Helgoland mud area, North Sea; (g) Skagerrak – S10; (h) Skagerrak – S11; (i) Skagerrak – S13; (j) Arabian Sea – OMZ core; (k) Arabian Sea – OMZ transition; (l) Arabian Sea – below OMZ; (m) Bering Sea; (n) Argentine Basin.

The Arabian Sea region (Fig. 5j–l) is characterized by the deposition of highly reactive OM, the concentration of which rapidly decreases with sediment depth (high v and low a; Table 6). The OM reactivity distributions show a clear depletion of the most reactive OM compounds in the uppermost 10 cm. Below 100 cm b.s.f. the remaining bulk OM is dominated by less reactive compounds (k(z) < 10⁻⁵ yr⁻¹). The rapid decrease in apparent OM reactivity reflects the often rapid deposition of predominantly marine OM (Cowie, 2005; Seiter et al., 2005), whose degradation during settling is limited by short residence times in combination with locally anoxic conditions (Luff et al., 2000). The shallow, sub-oxic and anoxic sediments of the northern European margin (Aarhus Bay – Fig. 5d; Arkona Basin – Fig. 5e) reveal a similar OM–k distribution and depth profile. However, the loss of the most reactive fractions within the uppermost 10 cm is less pronounced. Such a pattern may be attributed to the burial of terrestrial yet reactive OM (Aquilina et al., 2010). In contrast, Skagerrak sediments (S10 – Fig. 5g; S13 – Fig. 5i) are characterized by the most homogeneous, albeit comparably unreactive, OM mixture (high a values; Table 4). The OM–k distributions reveal small changes from the SWI down to 1000 cm b.s.f., reflecting a low overall reactivity, and thus low degradation rates. This somewhat surprising finding can be explained with the atypical character of the Skagerrak – a shallow environment that connects the brackish Baltic Sea with the saline North Sea and is characterized by an estuarine-like circulation. Due to this circulation pattern, the Skagerrak acts as the final sink for marine OM produced in the North Sea (Lohse et al., 1995; Van Weering et al., 1987). Consequently, more than 90 % of the OM deposited in the Skagerrak is derived from re-suspended and pre-degraded marine OM or terrestrial OM (ca. 20 % of the OM buried) (de Haas and van Weering, 1997) of generally low quality (Dauwe et al., 1999a, b).

Table 6 Benthic organic matter dynamics derived from the coupled model–data assessments across the depositional settings.

n.d.: not determined (oxygen and nitrate benthic fluxes not determined for locations where SWI concentrations of those TEA were assumed to be zero).

Between those two endmember cases (high-reactivity Arabian Sea and low-reactivity Skagerrak), sediments from other sites generally reveal small changes in the continuous OM–k distributions within the uppermost 10 cm. At all sites, k values have converged to < 10⁻⁴ yr⁻¹, below 100 cm b.s.f., where age(z) ≫ a, and changes in OM degradation rate constants with depth are thus largely controlled by age(z) (Eqs. 5, S8). As such, only the least reactive fractions are preserved, leading to an overall low apparent OM reactivity in deeper sediments that is largely independent of the initial apparent OM reactivity k(0). Sediments from the Bering Sea (Fig. 5m) are a notable exception. OM reactivity decreases gradually throughout the sediment column, reflecting both unusually high a and v values (Table 6). The Bering Sea is one of the most productive ocean regions, and OM deposition is mostly fuelled by spring and summer blooms (Coyle et al., 2008; Stabeno and Hunt, 2002; Stockwell et al., 2001). However, due to great water depths and low sedimentation rates, OM can experience extensive ageing and pre-processing prior to burial and result in the OM distribution pattern observed here (Fig. 5m).

4.1.4 Emerging environmental patterns in OM reactivity

Parameter a distributions exert the main control on the global variability in apparent OM reactivity. It also displays plausible values spanning several orders of magnitude (Fig. 4b, Table 6). Thus, our ability to quantify and predict benthic OM dynamics strongly depends on a general framework that allows constraint of parameter a when observational data are scarce or unavailable (e.g. global, future, and past model scenarios). Therefore, we further explore these links combining our newly generated dataset with those compiled from previously published RCM parameters (Fig. 6).

Figure 6 Depositional environment controls on organic matter reactivity parameters a (a–d) and k (e–g): (a) v vs. a; (b, e) water depth h; (c, f) sedimentation rate ω; (d, g) organic matter fluxes J_OM,in.

Despite the large scatter, broad trends emerge that generally agree with previously proposed relationships k–ω (Boudreau, 1997; Tromp et al., 1995), k–J_OM,in (Boudreau, 1997; Müller and Suess, 1979), and a–ω (Boudreau and Ruddick, 1991) (Fig. 6). Depositional environments characterized by low sedimentation rates, low OM fluxes, or great water depths generally display high a values (i.e. low apparent reactivity), whereas settings characterized by high sedimentation rates, high OM fluxes, or shallow water depth generally display lower a values (high apparent reactivity). While these trends are extremely weak on a global scale, they are relatively more pronounced within a given geographical setting (i.e. on a local/regional scale).

However, in agreement with the previous global assessment (Arndt et al., 2013), we do not find any statistically significant global or regional correlation between the reactivity parameters (a and k(0)) and master characteristics of depositional environments. Earlier proposed empirical relationships were based on a limited number of observations available at that time, and our results, as well as previously compiled parameter values, thus confirm cautious statements regarding the global validity of previously published relationships (Boudreau, 1997; Müller and Suess, 1979; Tromp et al., 1995). These results further emphasize that the spatial variability in a and thus apparent OM reactivity are not determined by one single characteristic of the depositional environment. Instead, it results from a complex and dynamic interplay of multiple environmental controls (e.g. Kayler et al., 2019; LaRowe et al., 2020b; Mayer, 1995; Zonneveld et al., 2010; Schmidt et al., 2011), and the ensemble of results indicates that sedimentation rates and OM deposition fluxes are poor proxies for the complex interplay of environmental controls on OM reactivity (Fig. 6c, f).

However, Fig. 6b indicates that water depth seems to be a useful and easily accessible first-order proxy for this complex interplay of environmental controls on OM reactivity. Thus, it may serve as an acceptable predictor for first-order global reactivity patterns in the absence of more suitable information. At the global scale, inversely determined apparent OM reactivity for the deep, well-oxygenated sites (e.g. Bering Sea and Argentine Basin) fall on the lower end of the reactivity spectrum, whereas shallow coastal environments, such as the northern European region, reveal generally higher apparent OM reactivities. In addition, within the northern European region and along the Rhone delta transect, parameter a generally increases (i.e. decrease in apparent reactivity) with increasing water depth (Fig. 6b, e). These trends can be explained with the broad large-scale patterns of environmental conditions along the global hypsometry. In deeper depositional environments, prolonged pelagic degradation (Cram et al., 2018; Maerz et al., 2020; Müller and Suess, 1979; Weber et al., 2016), exposure to oxygen (Hartnett et al., 1998), and the continuous ageing during settling (Griffith et al., 2010) result in a low apparent OM reactivity in underlying sediments. In contrast, several processes occurring in shallow environments, such as rapid settling, intense lateral transport, limited degradation in the water column, dynamic reworking through intense macrofaunal activity, and fluctuating redox conditions (e.g. Aller, 1994; Canfield, 1994; Hedges et al., 1997; Burdige, 2007) contribute to higher apparent OM reactivities. Local reactivity–water depth trends in shallow environments are driven by some of the same controls.

The weak water depth–reactivity trends and numerous exceptions at both global and regional level caution against the uncritical use of simplistic k–h relationships. Further important controls often complicate and sometimes inhibit the use of water depth as a predictor for apparent OM reactivity. For instance, the Severn estuary is an exception to this observed local, shallow environmental trend. This region is characterized by one of the largest tidal ranges on Earth, resulting in pronounced environmental gradients along the dynamic land–ocean transition. OM of terrestrial, anthropogenic, and marine sources mixes along the gradient and interacts with high loads of fine, suspended particular material. Additionally, intense tidally driven deposition–resuspension cycles (Dyer, 1984; Jonas and Millward, 2010; Manning et al., 2010; Uncles, 2010) result in a continuous re-exposure of benthic OM, which might cause the lower apparent reactivity of the buried OM. The Arabian Sea is another example where the k–h relationship deviates from the general global trends. Here, reactivity is 2–3 orders of magnitude higher than that derived for other deep-sea sites. Instead, reactivity is rather comparable to values determined for shallow sites in the northern European margin and the Rhone delta. In addition, across the Arabian Sea OMZ transect, only a modest regional depth trend is observed, with a decrease in reactivity (i.e. increase in a) from the OMZ to increasing water depths and bottom water oxygen concentrations (Cowie, 2005; Luff et al., 2000). While distinct from the overall k–h discussed above, the Arabian Sea OMZ region exhibits relatively similar patterns to other highly productive, anoxic water column systems, e.g. Peru margin (Boudreau and Ruddick, 1991); Blake ridge (Marquardt et al., 2010; Wallmann et al., 2006); and California, Chile, Costa Rica, and Namibia margins (Marquardt et al., 2010). This suggests a general separation from anoxic and oxic water columns shaping benthic OM reactivity.

4.2 Influence of OM reactivity on benthic–pelagic coupling

Both quantity and quality of the deposited OM exerts an important control on benthic–pelagic coupling processes with potential implications for global biogeochemical cycles and climate. While the influence of OM quantity on benthic–pelagic exchange along the global hypsometry has been previously investigated (Krumins et al., 2013; Soetaert et al., 1996, 1998; Thullner et al., 2009), the role of reactivity on benthic–pelagic fluxes remains less constrained.

4.2.1 Control of OM reactivity on OM degradation rates and respiration pathways

In agreement with previous findings (Henrichs, 1992), depth-integrated OM degradation rates (R_TOC, total; Eq. 7) show a strong correlation with OM fluxes (J_OM,in) (Fig. 7a, r²=0.97; p<0.01; n=14). In contrast, they do not reveal a strong link with apparent OM reactivity (i.e. parameter a), thus indicating that the magnitude of the OM deposition flux rather than its quality often exerts a first-order control on depth-integrated rates of OM degradation. In our dataset, apparent OM reactivity only becomes a dominant control on R_TOC, total in depositional environments that receive exceptionally reactive OM (i.e. a < 10¹ years), such as the Arabian Sea region (Table 6).

Figure 7 Benthic organic matter degradation dynamics. (a) Organic matter fluxes vs. depth-integrated organic matter degradation rates; (b) relative contribution of each metabolic pathway to heterotrophic organic matter degradation.

OM reactivity does exert an important impact on the relative contributions of metabolic pathways (Fig. 7). Overall, rates of OM degradation coupled to consumption of oxygen (aerobic OM degradation pathway) and nitrate (denitrification) make a small relative contribution (<10 %) to total degradation rates in the investigated oxic depositional environments (Rhone, Skagerrak, Bering Sea, and Argentine Basin; Fig. 7b). The relative contributions quantified here are generally lower than previous estimates for coastal and shelf sediments (>20 %) (e.g. Long Island Sound, Mackin and Swider, 1989; Skagerrak, Canfield et al., 1993; eastern Canadian continental margin, Boudreau et al., 1998; Barents Sea, Freitas et al., 2020; and northern Gulf of Mexico, Owings et al., 2021). Additionally, estimates are consistently lower than values in the deep sea (>70 %) (e.g. Cape Hatteras continental rise; Heggie et al., 1987), as well as over the global-scale hypsometry (Thullner et al., 2009). However, the overall dynamic nature of the considered environments, which favours the dominance of organoclastic sulfate reduction (Bowles et al., 2014; Jørgensen and Kasten, 2006; Thullner et al., 2009), as well as the length of sediment column (1000 cm b.s.f.) considered here renders a direct comparison difficult. Here, the contribution of the aerobic OM degradation pathway further decreases with a decrease in OM apparent reactivity (i.e. increase in a) as a consequence of the resulting higher burial of OM to sub-oxic or anoxic sediments (e.g. Meister et al., 2013). Additionally, the enhanced production of reduced species at depth further decreases the relative contribution of the aerobic OM degradation pathway to overall OM degradation rates by channelling more oxygen consumption through the re-oxidation of upward-diffusing reduced species (e.g. ammonium and hydrogen sulfide) (Glud, 2008; Hensen et al., 2006).

It has been suggested that anoxic metabolic pathways represent significant pathways of OM degradation (e.g. Bowles et al., 2014; Jørgensen and Kasten, 2006; Thullner et al., 2009). Our results show that sulfate reduction is the main oxidative pathway in regions characterized by high apparent OM reactivity (a < 10² years) and intermediate OM fluxes (J_OM,in ∼ 10⁻⁴–10⁻³ mol C cm⁻² yr⁻¹), e.g. Arabian Sea region (Cowie, 2005; Luff et al., 2000) and Arkona Basin and Aarhus Bay (Dale et al., 2008b; Jørgensen et al., 2019a; Mogollón et al., 2012). Additionally, high relative contributions of sulfate reduction are also simulated for environments that are characterized by a reduced supply of less reactive (a > 10³ years) OM deposition fluxes, such as the Argentine Basin and Bering Sea sites (e.g. Hensen et al., 2003). While high apparent OM reactivity enhances sulfate consumption in the upper sediment, where sulfate is efficiently replenished by exchange with bottom waters, a low OM reactivity allows for a deeper sulfate penetration (Fig. 3m–n) that sustains sulfate reduction over a greater sediment volume (e.g. Bowles et al., 2014; Meister et al., 2013; Riedinger et al., 2005, 2014, 2017). In contrast, methanogenesis is an important degradation pathway in environments that are characterized by a high supply ($J_{\mathrm{OM},\mathrm{in}}\mathit{>}{\mathrm{10}}^{-\mathrm{3}}$–10⁻² mol C cm⁻² yr⁻¹) of both highly reactive (a=10⁰–10¹ years) OM (Rhone proximal zone and Helgoland mud area, North Sea) and the least reactive (a>10³ years) OM (Skagerrak). The fast burial of reactive OM drives the development of a shallow methanogenic zone in the sediment (shallow sulfate penetration; Fig. 3b, f) (e.g. Aromokeye et al., 2020; Oni et al., 2015a, b; Rassmann et al., 2020), whereas very low apparent reactivity OM escapes sulfate reduction and is degraded deeper in the sediment via methanogenesis (e.g. Dale et al., 2008a; Knab et al., 2008; Meister et al., 2013). In both cases, the produced methane diffuses up and supports an enhanced consumption of sulfate through AOM (Sect. 4.2.2). Consequently, the relative contribution of organoclastic sulfate reduction decreases due to the presence of methane (e.g. Bowles et al., 2014; Niewöhner et al., 1998; Regnier et al., 2011; Jørgensen and Kasten, 2006; Jørgensen et al., 2019a; Riedinger et al., 2005, 2014, 2017).

We do not consider the degradation pathways mediated by metal oxides due to a lack of SWI data that would allow constraint of boundary conditions of those processes. This might result in a slight overrepresentation of sulfate reduction (e.g. at the Skagerrak; Canfield et al., 1993; Rysgaard et al., 2001). It has been demonstrated that metal oxide pathways represent a minor contribution to total OM heterotrophic degradation on a global scale (Thullner et al., 2009), although in continental margin sediments iron hydroxide reduction may represent an important OM degradation pathway (Beckler et al., 2016).

4.2.2 Anaerobic oxidation of methane dynamics

AOM occurs at the sulfate–methane transition zone (SMTZ), which can vary in sediment depth from centimetres to hundreds of metres as a function of different environmental controls, such as OM quantity and quality, sedimentation rate, active fluid flow, or microbial growth dynamics (Puglini et al., 2020; Regnier et al., 2011). Reflecting the diversity of the depositional environments studied here, we observe a broad range of SMTZ depths (Fig. 8a–c). Generally and in agreement with previous findings (Egger et al., 2018; Regnier et al., 2011), we observe a deepening of the SMTZ with decreasing sedimentation rates in addition to an important control of OM reactivity on the depth of the SMTZ (Fig. 8b) and the depth-integrated AOM rates (ΣAOM, Fig. 8d).

Figure 8 Anaerobic oxidation of methane (AOM) patterns as a function of apparent organic matter reactivity and depositional rates. (a) AOM rates vs. SMTZ depth; SMTZ depths vs. (b) reactivity parameter a and (c) sedimentation rates ω; AOM rates vs. (d) reactivity parameter a and (e) sedimentation rates ω.

Figure 9 Benthic fluxes of dissolved species (a–e) and transport mechanism relative contributions (f–j): (a, f) oxygen, (b, g) nitrate, (c, h) sulfate, (d, i) ammonium, and (e, j) phosphate.

The lowest ΣAOM and deepest SMTZ depths are determined for sites that are characterized by low apparent OM reactivity (i.e. deep-sea, well oxygenated bottom waters regions; Table 6). In those regions, low OM reactivity results in low rates of heterotrophic degradation (Fig. 7a) and, thus, deep ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ penetration (Fig. 3m–n) or no ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ depletion at all. In contrast, high OM reactivity in combination with high depositional fluxes (Fig. 8b–e) result in the complete depletion of TEAs in the upper sediment (z < 50 cm), allowing the establishment of methanogenesis at relatively shallow sediment depth. For instance, sediments from the Helgoland mud area, North Sea, are characterized by a shallow SMTZ and intense ΣAOM (Oni et al., 2015a, b). Yet, high OM reactivity does not always result in shallow SMTZ and intense AOM rates. Associated with moderate depositional fluxes, high apparent OM reactivity reduces the amount of reactive OM that reaches the methanogenic zone, inducing a deepening of the SMTZ and a decrease in ΣAOM. For instance, sediments from the Aarhus Bay (Flury et al., 2016) and Arkona Basin (Mogollón et al., 2012) display high OM reactivity (a < 10¹ years), which supports high rates of degradation close to the SWI (Fig. 7a) where ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is constantly replenished from bottom waters and limits OM fluxes to the deeper methanogenic sediment (e.g. Flury et al., 2016; Meister et al., 2013; Regnier et al., 2011). For moderate depositional fluxes, associated with lower OM reactivity (a > 10² years), e.g. the Skagerrak sites, pre-aged OM buried in these sediments (Aquilina et al., 2010) is predominantly degraded via methanogenesis (Fig. 7b). Consequently, a shoaling of the SMTZ is observed (Dale et al., 2008a; Knab et al., 2008), and higher ΣAOM values are determined (Table 6).

4.2.3 OM reactivity controls on sediment–water interface fluxes patterns

Model results suggest that apparent OM reactivity also influences nutrient recycling at the seafloor. Thus, we estimate benthic–pelagic fluxes of dissolved species (O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$). We emphasize that these flux calculations are based on simulated depth profiles at steady state and are thus only weakly affected by sampling limitations (e.g. possible loss of uppermost sediment layers associated with the choice of sampling devices; Hensen et al., 2006). Our estimates allow us to investigate the links between nutrient recycling and OM reactivity and to assess qualitative differences across sites.

Apparent OM reactivity broadly controls the spatial patterns of benthic–pelagic fluxes, as well as the relative importance of different transport processes (Glud, 2008; Bourgeois et al., 2017). Overall, molecular diffusion is the main transport pathway (Fig. 9f–j). However, the relative importance of bioirrigation increases at environments characterized by low deposition of OM and less reactive OM. This occurs because of the relative decrease in diffusive fluxes as a consequence of weak concentration gradients across the SWI. Bioturbation fluxes are generally low and reflect our assumption of depth-dependent D_bio (Middelburg et al., 1997). Advective fluxes become an important transport mechanism driving downward sulfate fluxes, especially in rapidly accumulating sediments and at sites characterized by low OM reactivity and moderate depositional rate.

Overall, high apparent OM reactivity drives high sedimentary uptake fluxes of dissolved TEA (Table 6; Fig. 9a–c). Sulfate largely dominates the benthic TEA uptake, particularly at sites where oxygen and nitrate are unavailable (Fig. 9c). At the Rhone proximal zone, unusually high sedimentation rates (Pastor et al., 2011) deliver fresh and highly reactive OM (Cathalot et al., 2010; Pruski et al., 2015), supporting high OM fluxes and degradation rates (Fig. 7a; Table 6). Consequently, intense oxygen consumption occurs at the SWI (Rassmann et al., 2016) driven by both aerobic degradation of OM and re-oxidation of reduced species. Similarly, in the Arabian Sea sediments below the OMZ the high input of fresh phytoplankton debris (Cowie, 2005; Koho et al., 2013; Rixen et al., 2019; Vandewiele et al., 2009) and potentially chemoautotroph biomass (Lengger et al., 2019) associated with intense degradation rates in the uppermost sediment layers (Luff et al., 2000) result in a high relative contribution of aerobic degradation of OM and denitrification to total OM oxidation. Consequently, both benthic oxygen and nitrate uptake fluxes are comparably high. In contrast, at regions characterized by deposition of less reactive OM (k < 10⁻⁴ yr⁻¹; Skagerrak, Argentine Basin, and Bering Sea), benthic–pelagic fluxes of TEAs are generally low.

Benthic fluxes play a crucial role in nutrient recycling. In the Arabian Sea region, large phytoplankton blooms are associated with monsoon conditions, which result in upwelling of nutrient-rich bottom water (SW monsoon) and deepening of the mixed layer (NE monsoon) (Cowie, 2005; Luff et al., 2000; Rixen et al., 2019). Similarly, across the northern European margin, spring and summer diatom blooms are common (Fleming-Lehtinen and Laamanen, 2012; Jensen et al., 1990; Karlson et al., 1996; Lomstein et al., 1990; Wiltshire and Manly, 2004) and maintained by benthic nutrient fluxes. Additionally, bottom water upwelling and nutrient recycling are important mechanisms sustaining spring and summer PP at the Bering Sea (Coyle et al., 2008). Benthic ammonium and phosphate recycling fluxes mirror the reactivity and OM degradation patterns (Fig. 9d–e). The largest fluxes are determined for the Arabian Sea region, as well as the shallow Aarhus Bay and Arkona Basin (Table 6). In contrast, the lowest nutrient recycling fluxes are observed in regions characterized by the deposition of less reactive OM (phosphate fluxes are negligible at Skagerrak sites S10 and S13), where heterotrophic degradation rates are slow.

5 Conclusions

We developed and applied an inverse modelling approach to quantify apparent OM reactivity (i.e. parameters a and v) based on a common minimum set of benthic observational data, as well as a common RTM approach for 14 different sites across five different depositional environments. We provide important first-order constraints on the reactive continuum model (RCM) parameterization that can inform on the choice of parameter values in data-poor areas, including global-scale and timescale scenarios.

Our findings corroborate previous results (Arndt et al., 2013; Boudreau and Ruddick, 1991; Forney and Rothman, 2012; Middelburg, 1989) that the RCM parameter v is globally relatively constant (v=0.1–0.2). Exceptionally high v > 0.2 is often associated with deposition and burial of highly reactive phytoplankton debris in high-productivity regions associated with well-established OMZs. In contrast, in agreement with previous findings, RCM parameter a can span several orders of magnitude at a global scale, suggesting that the parameter a is the main parameter describing the variability of apparent OM reactivity. Consequently, future modelling efforts to quantify OM reactivity on a global scale could be reduced to one main reactivity variable. Based on our results and previous findings, we consolidate the range of predominant a distribution to 10⁰–10⁴ years. This is a valuable constraint when dealing with data-poor regions and timescales (e.g. Hülse et al., 2018), since it excludes extreme values at both ends of the a range.

Additionally, results indicate that the large variability in apparent OM reactivity is linked to a combination of multiple environmental drivers that control the quality of OM delivered to sediments and the timescale of settling and burial. Therefore, our findings contribute to the notion that apparent OM reactivity is controlled by a complex and dynamic interplay of environmental drivers, which are measurable and allow the quantification of a, v, and k. Our model results caution against the use of a single environmental master variable such as water depth, sedimentation rate, or OM deposition fluxes beyond well-defined geographical scales. Yet, results show that, if no other information is available, water depth can serve as a useful proxy for the complex and dynamic interplay of environmental drivers and can be used to predict first-order, large-scale OM reactivity patterns. However, more holistic, interdisciplinary exploration of OM reactivity in its entire environmental context is needed to mechanistically understand benthic carbon and nutrient dynamics.

Finally, results also show that OM reactivity exerts a dominant control on the redox zonation of the sediment, the depth of the SMTZ, depth-integrated AOM rates, and benthic–pelagic exchange fluxes, indicating that these processes could serve as predictor variables for OM reactivity. In contrast, depth-integrated OM degradation rates are largely controlled by the magnitude rather than the quality of OM deposition to the sediment.