Advancing on large-scale trends of apparent organic matter reactivity in marine sediments and patterns of benthic carbon transformation

Abstract. Constraining the mechanisms that control 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 quantitative understanding of what controls OM reactivity in marine sediments and, as a result, how to constrain it in global models. To fill this gap, we quantify apparent OM reactivity (i.e., model-derived estimates) by extracting reactive continuum model parameters (a and v) from observed benthic organic carbon and sulfate dynamics across 14 contrasting depositional settings distributed over five distinct benthic provinces. Our analysis shows that the large-scale range in apparent OM reactivity is largely driven by the wide variability in parameter a (10−3 



Introduction 50
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 y -1 ) or marine environments (global NPP = ~49 Pg C y -1 ). In addition, it also comprises OM derived from ancient PP, which is delivered through rock weathering, soil erosion, and thawing permafrost, as well as chemoautotrophy, and secondary production by microbes and animals (e.g., Hedges 55 et al., 1997;Hedges and Oades, 1997). 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). This small imbalance between OM production and degradation connects the short-and long-term carbon cycles, controls the longterm evolution of atmospheric CO2 and has enabled the accumulation of oxygen in the atmosphere (Berner, 2003;Hedges, 1992) . Yet, compilations of field data reveal that OM degradation, and preservation rates vary significantly in space (Arndt et 60 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 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 an inherent characteristic of the OM itself, but depends on its environmental context (e.g., Kayler et al., 2019;LaRowe et al., 2020b;Mayer, 1995;Zonneveld et al., 2010). In fact, OM 65 reactivity is determined by the dynamic interplay of a plethora of different factors that include, but are not limited to OM sources and composition (Burdige, 2005;Meyers, 1997); oxygen exposure time during settling and burial 3 (Hartnett et al., 1998;Hoefs et al., 2002;Huguet et al., 2008;Sinninghe Damsté et al., 2002;Sun and Wakeham, 1994;Wakeham et al., 1997); terminal electron acceptor (TEA) availability (Aller, 1994;Aller and Aller, 1998); microbial activity (LaRowe et al., 2020b); sediment biological mixing (Aller, 1980(Aller, , 1994Aller and Aller, 1998;Boudreau, 1994;Grossi et al., 70 2003;Middelburg, 2018); OM ageing and transport history (Bianchi et al., 2018;Cathalot et al., 2010;Griffith et al., 2010;Mollenhauer et al., 2003Mollenhauer et al., , 2007Ohkouchi et al., 2002); mineral protection through OM-sediment association (Bianchi et al., 2018;Hemingway et al., 2019;Keil and Cowie, 1999;Mayer, 1994aMayer, , 1994b; adsorption onto reactive iron minerals (e.g. Faust et al., 2020;Ma et al., 2018;Salvadó et al., 2015;Shields et al., 2016;Wang et al., 2019); or hydrogen sulfide exposure time (i.e., sulfurization) (Sinninghe Damsté et al., 1988, 1998. However, the relative significance of each of these factors 75 remains poorly quantified. As a consequence, we currently lack a general framework that allows quantifying 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 Lessin et al., 2018), quantifying both pathways and rates of biogeochemical processes over different temporal and spatial scales (Hülse et al., 2018) and, thus, ultimately predicting the impact and feedbacks of past, present, and future climate change feedbacks on global biogeochemical cycles and climate. 80 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 meters deep and hundreds to thousands and million years old) (e.g. Boudreau and Ruddick, 1991;Middelburg, 1989;Wehrmann et 85 al., 2013), using Reaction-Transport Models (RTM). Especially 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/ages. RTMs can be broadly subdivided into two classes regarding OM conceptualisation: discrete models and continuum models (e.g., Arndt et al., 2013). Discrete or G-type models (Berner, 1980), and particularly multi-G models (Jørgensen, 1978;Soetaert et al., 1996;Westrich and Berner, 1984) compartmentalise OM into discrete pools each having 90 distinct reactivity. Although mathematically simple, multi-G models require the definition of both the size and reactivity of its discrete OM pools. Observational data generally allows unequivocal identification of a maximum number of three distinct compound classes (Jørgensen, 1978;Middelburg, 1989). Continuum models (Boudreau and Ruddick, 1991;Middelburg, 1989) on the other hand are mathematically more complex models, which assume a continuous distribution of OM compounds over the entire reactivity spectrum, thus avoiding the partitioning of OM into a discrete compound classes. Unlike Multi-G models 95 that converge towards a constant reactivity once more reactive compound classes are depleted at depth, continuum models account for the widely observed continuous decrease in OM reactivity with sediment depth. Consequently, the application of multi-G models is limited to shallow, contemporary sediments (Ait Ballagh et al., 2020;Dale et al., 2008aDale et al., , 2008bDale et al., , 2015Luff et al., 2000;Paraska et al., 2014;Soetaert et al., 1996;Stolpovsky et al., 2018), whereas continuum models are successfully employed in both shallow, contemporary environments, as well as deeply buried, old sediments (Arndt et al., 2009;Boudreau 100 and Ruddick, 1991;Bradley et al., 2020;Contreras et al., 2013;Henkel et al., 2012;Jørgensen et al., 2019a;LaRowe et al., 2020a;Mogollón et al., 2012Mogollón et al., , 2016Wehrmann et al., 2013;Ye et al., 2016). However, due to our limited quantitative understanding of the environmental controls on OM reactivity, the application of both model types to data-poor settings is still complicated by the difficulty in constraining model parameters, inducing important uncertainties at the global scale as well as past and future applications. 105 In a first attempt to address this knowledge gap, Arndt et al. (2013) compiled a comprehensive global dataset of previously published, model-derived apparent OM reactivities across a wide range of depositional environments. Their analysis revealed that reactivity is largely variable within each sedimentary setting with no clear global trends. Nevertheless, broad regional reactivity trends were highlighted. For instance, the highly productive upwelling region of the Arabian Sea is characterized by 110 the burial of highly reactive OM. Monsoon-induced high PP rates and an associated pronounced oxygen minimum zone (OMZ) promote an efficient export and burial of fresh, reactive OM (Cowie, 2005;Koho et al., 2013;Luff et al., 2000;Rixen et al., 2019;Vandewiele et al., 2009). In coastal settings and continental margins, a high variability in OM reactivity is attributed to variable inputs of terrestrial and marine OM with complex degrees of alteration prior to burial Oni et al., 2015b). In open ocean sediments, a large range in the reactivity spectra arises from several processes operating simultaneously 115 at great water depths; e.g., efficient degradation during settling (Emerson and Hedges, 1988;Sun and Wakeham, 1994;Wakeham et al., 1997), extensive exposure time to oxygen (Hartnett et al., 1998;Mewes et al., 2014Mewes et al., , 2016Mogollón et al., 2016;Volz et al., 2018), and ageing (Griffith et al., 2010). The global compilation of OM reactivity  also revealed a lack of significant global correlation between OM reactivity and single characteristics of the depositional environment (water depth, sedimentation rate, OM fluxes), thus confirming cautious statements regarding the global validity 120 of previously published relationships (Boudreau, 1997;Müller and Suess, 1979;Tromp et al., 1995) and emphasising the role of a complex interplay of environmental controls on reactivity (e.g., Kayler et al., 2019;LaRowe et al., 2020b;Mayer, 1995;Zonneveld et al., 2010).
Crucially, the analysis by Arndt et al. (2013) revealed that a lack of comparative studies was a roadblock towards developing 125 a robust global scale framework for OM reactivity. The large variability in model formulation and approaches adopted by different studies limits the comparability and transferability of results , which hampers our ability to develop a systematic understanding of processes controlling OM reactivity on a global scale and OM cycling in the past, present, and future (Hülse et al., 2018). Therefore, we here use a consistent, integrated model-data approach to inversely determine the parameters that control apparent OM reactivity in the Reactive Continuum Model (RCM), and (Boudreau and Ruddick, 130 1991) from contemporary sediment depth profiles (Sect. 3) at 14 different sites from five different depositional environments (Sect. 2). Briefly, those parameters, alongside burial age, determine the depth-evolution of apparent OM reactivity; and consequently, OM degradation dynamics (see Sect. 3.2 for details). The large-scale variability in inversely determined parameters and is then assessed and first order constraints on model parametrization are provided. Additionally, following the approach of Arndt et al. (2013), we compile and values derived from previously published data-model approaches 135 5 (Arndt et al., 2009;Boudreau and Ruddick, 1991;Contreras et al., 2013;Henkel et al., 2012;Mewes et al., 2016;Wehrmann et al., 2013;Ye et al., 2016). Inversely determined parameters are compared to previously published values and the extended dataset is then used to explore the links between OM reactivity and master characteristics of depositional environments. Furthermore, we quantify key processes associated with OM transformation within sediments to assess the role of OM reactivity on benthic biogeochemical cycling. From our model simulations, we estimate OM degradation 140 rates and relative contributions of metabolic pathways, anaerobic oxidation of methane coupled to sulfate reduction (AOM) rates, as well as dissolved TEA and nutrient fluxes across the SWI. With this comprehensive and global-scale approach, we aim to explore the links between OM reactivity and characteristics of depositional environments, and to understand the controls of reactivity on OM transformation processes.

Study sites 145
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 1). 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 characterised by predominantly temperate climate regimes, except for the Arabian Sea region (Arid) and the Bering Sea area (Snow/Polar) 150 (Chen and Chen, 2013). OM characteristics also cover a broad spectrum, ranging from organic rich ( > 2 wt%) to lean ( < 1 wt%) sediments (Premuzic et al., 1982;Seiter et al., 2004), and from predominantly terrestrial (coast), mixed terrestrial and marine (coast/shelf) to predominantly marine (shelf/slope) OM 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, which are triggered by various mechanisms (Table 1). This is also reflected by the 155 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 (Table 1) offers a unique opportunity to quantify OM reactivity and explore its multidisciplinary links with environmental drivers.

Methods
Here, we adopt an integrated model-data approach to inversely determine OM reactivity parameters and from 160 contemporary and 4 2− profiles measured at 14 sites across the five different depositional environments described in  Thomas (2014). Rhone delta data (Table 1; 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 1; rows D, E, G, H, and I) are part of the European funded project METROL (Aquilina et al., 2010;Dale et al., 2008aDale et al., , 2008bKnab et al., 2008;Mogollón et al., 2012). Data from the Helgoland Mud Area, North Sea (Table 1; row F) (Henkel and Kulkarni, 2015), the Arabian Sea transect (Table 1; rows J, K, 175 and L) (Bohrmann and cruise participants, 2010a(Bohrmann and cruise participants, , 2010b(Bohrmann and cruise participants, , 2010c, Bering Sea (Table 1; row M) (Gersonde, 2009), and Argentine Basin (Table 1; row N) (Henkel et al., 2011) are presented on PANGAEA database (see above). We use those sitespecific datasets to inform our data-model analysis and to constrain OM reactivity and related benthic processes (Sect. 3.2.2 and 3.2.3).

Model description
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 185 Jourabchi, 2005;Thullner et al., 2009). Here, we design a BRNS set-up that explicitly resolves the fluxes and transformations of OM ( 2 ), the most pertinent TEAs (oxygen -2 ; nitrate -3 − ; and sulfate -4 2− ) and reduced species (methane -4 ; sulfide -− ; and ammonium -4 + ). Due to the limited availability of data to constrain manganese oxide ( 2 ) and iron hydroxide ( ( ) 3 ) depositional fluxes, the model does not account for metal oxides reduction pathways. At certain depositional settings (e.g., Skagerrak), metal oxides can be relatively important TEAs (e.g., Canfield et al., 1993;Rysgaard et 190 al., 2001), particularly in continental margins that receive considerable inputs of iron (Beckler et al., 2016). Nevertheless, a previously published global assessment of the importance of metabolic pathways in marine sediments has found their contributions to the overall heterotrophic OM degradation to be negligible at a global scale (Thullner et al., 2009). For each of the species involved in heterotrophic OM degradation, the BRNS simulates changes in both solid phase and porewater concentrations as result of transport processes (advection, molecular diffusion, bioirrigation and bioturbation), as well 195 production/consumption due to reactions (Aguilera et al., 2005;Thullner et al., 2009;Wehrmann et al., 2013). In the following sections we provide a detailed description of the model approach, including parameters values (Table 2), boundary conditions (Table 3), rates (see Table S1 for stoichiometric equations) and fluxes calculations.

Model formulation 200
The concentration depth profiles of solid and dissolved species in marine sediments are calculated according to the verticallyresolved mass conservation equation of solid and dissolved species in porous media (Berner, 1980;Boudreau, 1997): where is the concentration of specie , denotes time, and is the sediment depth. For solid species, the porosity term is given by = (1 − ), whereas for dissolved species porosity assumes = . The effective molecular diffusion coefficient of dissolved species is given by ( = 0 for solid species), represents the bioturbation coefficient, the sedimentation rate, and denotes the bioirrigation coefficient ( = 0 for solid species). The sum of consumption/production process rates is given by ∑ , where the stoichiometric coefficient of specie is given by for the kinetically controlled 210 reaction , with rate .

Transport processes and parameters
The RTM accounts for sediment accumulation and compaction, molecular diffusion, bioturbation, and bioirrigation (see Eq. 1). Sediment porosity is assumed to decrease exponentially with depth due to sediment compaction: 215 where is the porosity attenuation coefficient, is the porosity at depth and ∞ is the porosity at greater depth.
Consequently, the burial velocity is corrected for the effect of compaction assuming steady state compaction (e.g. Berner, 220 1980): The diffusive fluxes are commonly quantified by means of Fick's first law of diffusion, which depends on the molecular 225 diffusive coefficient ( = dissolved species) and is corrected for sediment tortuosity. is dependent on the charge and size of the diffusing species, as well as temperature and viscosity of the medium (Boudreau, 1997;Burdige, 2006). Here, are corrected for temperature, salinity, and tortuosity. The model also accounts for the effect of sediment reworking by infaunal 8 organisms in the bioturbated upper sediment layer ( < ). The process is generally described by a dispersive term with constant bioturbation diffusion coefficient (Boudreau, 1986). Bioirrigation, i.e. the mixing by benthic macrofaunal 230 organisms that build burrows or tubes in the sediment for feeding is described as a nonlocal transport process with a nonlocal bioirrigation coefficient , which describes the exchange rates between SWI and porewater at depth in bioirrigated zone of sediments (Aller, 1994;Aller and Aller, 1998;Burdige, 2006). Implemented transport parameters can be subdivided into two classes: (i) global model parameters (Table 2), which are 235 universal for all depositional environments, and thus assure model transferability across sites; and (ii) site-specific transport parameters (Table 3), which are constrained for each depositional environment investigated here, either by direct observations or derived from empirical relationships.

Reaction network 240
The reaction network implemented here encompasses the most pertinent primary and secondary redox reactions found in the upper layers of marine sediments. Its formulation and parametrization builds on a number of previous studies that investigate diagenetic dynamics across several depositional environments and scales (Aguilera et al., 2005;Thullner et al., 2009;Van Cappellen and Wang, 1996;Wang and Van Cappellen, 1996;Wehrmann et al., 2013). It explicitly accounts for the heterotrophic degradation of OM coupled to the consumption of oxygen (aerobic OM degradation), nitrate (denitrification), 245 sulfate (organoclastic sulfate reduction), as well as methanogenesis. Additionally, it accounts for nitrification, sulfide reoxidation by 2 (in the absence of ( ) 3 , which would most likely be responsible for re-oxidation of sulfides; e.g., Findley et al. (2020)), anaerobic oxidation of methane coupled to sulfate reduction (AOM) and 4 reoxidation by 2 . Stoichiometric reactions and kinetic rate formulations for each of the primary and secondary redox reactions involved in OM heterotrophic degradation are summarized in the Supplements (Table S1 and Table S2). 250 OM degradation rates are calculated assuming that the bulk OM is continuously distributed over a range of reactivities , i.e. a reactive continuum model (RCM). OM is degraded according to first-order kinetics with respect to the electron donor, (i.e. particulate organic matter (POM) which is assumed to be ( 2 ) 106 ( 3 ) 16 ( 3 4 )). A recent study on the sensitivity of OM stoichiometry towards degradation dynamics (Arning et al., 2016) supports the use of 2 simplification for 255 describing complex OM composition. As such, the overall rate of OM degradation for a continuous distribution of organic compounds is given by the integral: 9 where ( , ) denotes the probability density function that determines the concentration of OM having a degradability between and + at time , with being equivalent to the first-order degradation rate constant. At = 0, the initial OM distribution ( , 0) can assume distinct mathematical formulations, however it cannot be constrained from observations.
Here, a Gamma function is adopted (Aris, 1968;Boudreau and Ruddick, 1991;Ho and Aris, 1987). Assuming first-order degradation kinetics, the initial distribution ( = 0) of OM is given by: 265 where (0) denotes the initial OM content, is the Gamma function, is the average lifetime of the more reactive components of bulk OM, and represents the dimensionless scaling parameter of the distribution near = 0. The free, 270 positive parameters and delineate the shape of the initial distribution of OM compounds along the range of , 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 .
Due to the rapid depletion of the most reactive compounds, the reactivity of the bulk material decreases during degradation 275 and thus reflects the widely observed reactivity decrease with burial time/depth/age (Boudreau and Ruddick, 1991;Middelburg, 1989). In the RCM, OM compounds are continuously and dynamically distributed over a range of reactivities that captures the decrease in apparent reactivity with burial age/depth as the most reactive compounds are successively degraded.
The interplay of and drives the bulk OM reactivity depth profiles and consequently the yielded profiles from the reaction network implemented in the BRNS. The Gamma distribution captures the OM degradation dynamics in nature. As such, the 280 evolution of OM concentration as a function of depth ( ( )) is given by: where and are the free RCM parameters Boudreau and Ruddick, 1991) and ( ) denotes the age of 285 the sediment layer at depth . For non-bioturbated sediments the burial ( = − ) can be calculated as a function of the burial velocity and is given by: 290 where 0 is the initial age at the bottom of the bioturbated zone, is the porosity attenuation coefficient, is the porosity at depth and ∞ is the porosity at greater depth. However, within the bioturbated upper sediment layers, the age distribution of reactive species is controlled by both sedimentation, bioturbation, and the reactivity of reactive species (Meile and Van Cappellen, 2005). Thus, we here apply a multi-G approximation for the RCM in the bioturbated sediments (Sect. 3.2.1.3).

295
The sequential utilisation of TEA is described by a combination of Michalis-Menten and inhibition terms. The rate of heterotrophic OM degradation is independent of the TEA concentration for TEA > ( = TEA). However, for low TEA concentrations, the specific OM degradation pathway rate becomes limited by TEA availability, expressed as a first-order dependency of the rate with respect to the TEA (Table S2). TEAs are consumed sequentially following the decrease in Gibbs free energy yield of the respective metabolic pathway reactions. The resulting classic redox sequence is described by a sequence 300 of inhibition terms , where denotes the TEA consumed by the respective metabolic pathway: After all TEAs are depleted, OM degradation proceeds via methanogenesis. In addition to primary redox reactions, the reaction network also accounts for secondary redox reactions, i.e. re-oxidation of reduced species produced during primary redox reactions. Following the classical approach, the rates of secondary redox reactions are described by biomolecular rate laws with the rate constant (Table 2 and Table S1) . 310

Multi-G approximation of RCM in the bioturbated layer
Within the bioturbated layer ( < ), the RCM is approximated by a multi-G model (200 fractions) to circumvent the difficulty of quantifying OM ages within bioturbated sediments (Dale et al., 2015;Meile and Van Cappellen, 2005). The multi-G model approach divides the bulk OM into multiple discrete compound classes , each degrading according to first-order 315 kinetics with a degradation rate constant (Jørgensen, 1978). It thus accounts for the heterogeneity of bulk OM by assigning different reactivities to each pool of organic matter ( ) . The RCM is approximated by dividing the initial, continuous distribution of OM compounds over the reactivity spectrum (Eq. 4) into discrete compound classes. The fraction of total OM in compound class can be directly calculated by integrating the initial probability density function (Eq. 4), which determines the concentration of OM having a degradability between and + at time zero with: 320 = +1 2 (10).

11
In bioturbated sediments, the RCM is then approximated by dividing the defined reactivity range = [10 −15 , 10 − log( )+2 ] into 200 equally spaced reactivity pools . The initial fraction of total in compound class defined by the reactivity 325 bin −1 and in the 200-G model can be calculated as: The least and the most reactive fractions 1 and 200 with reactivity 1 = 10 −15 and 200 = 10 − log( )+2 y -1 , respectively, are 330 calculated based on the upper incomplete Gamma function: . 335 Once and are determined, the steady-state solution of the diffusion-advection-reaction equation (Boudreau, 1997) for OM in the bioturbated zone is then given by the general analytical solution: A and B denote integration constants that can be determined by defining appropriate boundary conditions (Boudreau, 1997) 350 for OM at the upper and lower boundaries. Below the bioturbated zone, the depth evolution of is determined by the RCM formulation (Eq. 6).

Additional model parameters
The bioirrigation rate was calculated from the bioirrigation coefficient at sediment surface 0 and the bioirrigation 355 attenuation depth length (Thullner et al., 2009;Wehrmann et al., 2013): The effective molecular diffusion coefficients for each of the species (Table 2) are derived from Van Cappelen and Wang 360 (1996). values are corrected for temperature ( ), salinity ( ) and tortuosity ( 2 ). For solid species = 0, whereas for dissolved species the corrected * is given by Boudreau (1997) Sediment porosity 0 , ∞ , , sedimentation rate , and temperature were either measured at site or derived from literature when actual measurements were not available. Salinity was defined as 35 at all sites. Bioturbation depth was fixed at 10 cm for most bioturbated sediments (Table 3), based on a compilation of mixed layer depths (Boudreau, 1994(Boudreau, , 1998. At anoxic depositional environments (i.e., where 2 concentrations are zero), was set to zero (Table 3). Similarly, the bioturbation 370 diffusion coefficient was set to zero at anoxic environments. For bioturbated sites, was constrained based on an empirically derived relationship proposed by Middelburg et al. (1997):

375
where ℎ denotes the water depth. For the Rhone delta, values were derived from Pastor et al. (2011). was kept constant within the bioturbated zone ( < ) then set to zero below .
The reaction parameters (except for the RCM parameters and ) implemented in the BRNS are universal for all depositional environments (Table 2). This approach ensures that model set-ups are consistent and model-derived RCM parameters are 380 independent of model formulation and fully comparable across depositional settings. The half-saturation concentrations of the respective Michaelis-Menten terms were adopted from Van Cappellen and . Secondary redox reaction constants for each explicitly resolved reaction are based on a compilation of published values . The RCM parameters and reflect the site-specific OM reactivity and are free parameters that are determined by inverse modelling (Sect. 3.3). 385

Boundary conditions
Boundary conditions place the BRNS in the environmental context of each of the study sites (Table 3). The boundary conditions are constrained based on either site-specific measurements or alternatively on published data if direct observations were not available. Here, we assume a fixed boundary concentration ( 0 ) of OM at the sediment water interface. 390 Concentrations of 2 are based on measurements in Skagerrak (Canfield et al., 1993) and Rhone delta (Rassmann et al., 2016) sediments. Similarly, concentrations of 3 − are based on measurements in Skagerrak sediments (Canfield et al., 1993).
Because of the lack of site specific information, these 2 and

Model solution
Transport and reaction equations were solved sequentially. Firstly, the diffusion term was discretized at each time-step of the numerical integration using the semi-implicit Crank-Nicholson scheme. This was followed by the calculation of the advective transport, using a 3 rd order accurate total variation diminishing algorithm with flux limiters (Regnier et al., 1998). The reaction network was subsequently solved. The mass-conservation equation (Eq. 1) was discretized on an uneven grid (Boudreau, 405 1997): where ( ) is the depth of the ℎ grid point, denotes the length of the model domain, is a point in a hypothetical grid, 410 and is depth relative to which ( ) is quadratically distributed for ≪ and linearly distributed for ≫ . and were chosen so that the grid size, ∆ , increases downcore from SWI to a maximum of . The size of model domain was fixed at 1,000 cm for all sites, except for the Bering Sea, in which the model domain is extended to 1,500 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 415 on biogeochemical dynamics within the model domain. BRNS was run until steady state (Δt < 0.01) was reached: The assumption of steady-state condition for estuarine sites may be impacted by the high dynamic nature of such settings. At 420 the Severn estuary a strong tidal dynamics results in continuous cycles of sediment resuspension (Manning et al., 2010). The Rhone delta experiences strong pulses of fresh water and sediments associated to flood events, as well as variable sedimentation rates (Antonelli et al., 2008;Cathalot et al., 2010;Zebracki et al., 2015). Despite that (and somewhat surprisingly), previous investigations at the Rhone delta found similar biogeochemical trends on an interannual basis, i.e. data collected along the years in this system are consistent Pastor et al., 2011;Rassmann et al., 2016Rassmann et al., , 2020. As such, it is likely 425 that the majority of sediment load and OM are delivered during late fall and early winter flood events, and the diagenetic system close to the river mouth reaches a mature state during early spring when the typical porewater and benthic-pelagic fluxes have established. A surprisingly similar steady benthic cycling has been observed along a sea-sea-ice gradient at the Barents Sea shelf, which is a very seasonally dominated environment (Freitas et al., 2020). Thus, for the purpose of our investigations, we assume that on longer time/depth scales those sediments reach steady state and the limitations imposed 430 should be minimal.

Inverse modelling
Here, we apply the BRNS to determine RCM parameters ( and ) from present-day downcore observations by inverse modelling of observed depth profiles (e.g., Arndt et al., 2009;Wehrmann et al., 2013). The first efforts to quantify and 435 parameter values, and thus reactivity, on the basis of observations solely focused on downcore profiles (Boudreau and Ruddick, 1991). While this approach has provided useful insights, it is nevertheless compromised by several factors. First, it is generally difficult to quantify OM contents at the sediment-water interface (SWI) due to the common loss of the upper few centimetres during sampling with gravity or piston corers. In addition, a steady-state OM deposition over the entire time span of sediment deposition is unlikely, especially in slope and coastal settings. Finally, multiple and pairs could potentially fit 440 a given depth profile (Meister et al., 2013). Parameter estimates based on profiles alone are thus associated with a generally unknown degree of uncertainty.
Adding further constraints increases the robustness of predictions. OM heterotrophic degradation is the driver behind biogeochemical dynamics in marine sediments, and downcore profiles of TEAs (e.g., 4 2− ) (e.g., Bowles et al., 2014;445 Jørgensen et al., 2019b) or reduced products of OM degradation (such as 4 and determine a minimum set of observational data that is widely available and comparably easy to measure, we tested the suitability of different artificial porewater data sets for the inverse determination of apparent OM reactivity. Using an artificially generated dataset (see Supplements, Table 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 is considered as a single constraint, 450 multiple pairs of and produce equivalent depth-profiles (Supplements, Fig. S1), as previously suggested (Meister et al., 2013). Model results show that by adding 4 2− depth profiles as an additional constraint (Supplements, Fig. S2), the determination of and becomes more robust (Supplements, Fig. S3). 4 profiles were used as an additional constraint when the data were available, although measured ex situ 4 profiles are often not reliable due to sampling and measurement uncertainties in deeper sediment layers where concentrations exceed mM (Dale et al., 2008a;Hensen et al., 2003;Hilligsøe et 455 al., 2018). The dynamics of 4 2− and 4 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 be also affected by bioirrigation in particularly shallow sulfate-methane transition zones (SMTZs) (e.g., Dale et al., 2019). In anoxic settings and deeply buried sediments, 4 2− is the dominant TEA and 4 is the most common reduced species (Jørgensen et al., 2019a(Jørgensen et al., , 2019b. Thus, a combination of , 4 2− , and 4 (if available to verify the depths of the SMTZ) depth profiles incorporates the information contained in the 460 observed benthic sulfur and carbon dynamics and is sufficient to extract apparent OM reactivity and its evolution from the sediment-water interface down to the SMTZ. Thus, here we perform a site-specific data-model fit based on and 4 2− (and 4 ) depth profiles. The RCM parameters were inversely determined by running the model for each site over the entire range of previously published ( = 10 -3 -10 7 years) and ( = 10 -2 -10 0 ) values. First, the model was run over a coarsely resolved ensemble of and combinations to explore the order of magnitude in which the RCM parameters fall. Once the 465 order of magnitude was established, a new model ensemble was run over a finely resolved grid of and pairs falling in the previously defined range.

Quantification of organic matter reactivity and associated benthic processes
The apparent reactivity of bulk OM for each site is calculated based on the RCM reactivity parameters and (Boudreau 470 and Ruddick, 1991) that yield the best data-model fit and is given by: .
The apparent reactivity of bulk OM at the sediment water interface ( = 0) is thus given by: 475 The BRNS simulates downcore concentration profiles for each species that is explicitly resolved in the reaction network, as well as reaction rates . From those reaction rates, the total rate of organic carbon oxidation, i.e. the depth integral of for 480 each primary redox reaction can be calculated by integrating the rate over the entire model domain (e.g., Thullner et al., 2009):

485
The OM fluxes ( ) at the SWI are calculated based on the sum of ∑ and burial fluxes (e.g. Burdige, 2006): Additionally, the fluxes of species at the SWI, i.e. the benthic-pelagic fluxes, are calculated from the model-derived depth profiles of species according to:

Results and discussion
The integrated data-model analysis yielded a comprehensive picture of OM reactivity parameters ( , , and ( )), as well as 500 OM degradation dynamics and benthic-pelagic coupling for our large-scale compilation of depositional environments (Table   4). 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. In the following sections, we explore the inverse modelling results in the context of their depositional environment, (e.g. water depth, sedimentation rates, and OM fluxes, bottom water and sediment redox conditions, OM quality, and mineral protection), explore emerging patterns and 505 provide important, first order constraints on model parameterization.

Inverse Modelling
We inversely determined the set of RCM parameters and that produces the best model fit to the observed (Fig. 2) and 510 4 2− (Fig. 3) depth-profiles. Overall, profiles display a satisfying fit between model and observations. However, the quality of fit is often compromised by a lack of observations/data at the SWI due to core loss during sediment sampling by gravity and piston corers, as well as a poor downcore sampling resolution. Therefore, determining the reactivity parameters based solely on observations, an approach adopted in the past (Boudreau and Ruddick, 1991), would not be feasible for those sites and adding additional constraints relieves such limitation. The 4 2− depth-profiles (Fig. 3) display a good 515 agreement between present-day observations and model simulations, except for the Arabian Sea site below the OMZ (Fig. 3l) and the Bering Sea (Fig. 3m). At the Bering Sea, the data-model discrepancy is likely associated to kinks in both (Fig. 2m) and 4 2− (Fig. 3m) profiles. They likely represent changes in the depositional regime (non-steady state) 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 4 2− recycling through 2 oxidation by 2 and ( ) 3 (e.g., Rassmann et al., 520 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.

Apparent OM reactivity parameters and
The inversely determined OM reactivity parameters and (Eq. 6; Eq. 22) derived from the best fit model solutions (Table  525 4) reveal important differences in apparent OM reactivity and its distribution on a global scale. Alongside and 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., 2012Mogollón et al., , 2016Wehrmann et al., 2013;Ye et al., 2016), they provide useful constraints on the plausible range of and , as well as their control on the global variability in apparent OM reactivity. Thus, they provide important first order constraints for model parameterization in data poor areas. 530 Parameter exerts a scaling influence on OM reactivity distribution. High -values result in a high initial OM reactivity, whereas low -values produce low initial reactivity . For instance, an extremely low -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 535 published values for contemporary ocean sediments fall within the range 0.1 and 1.08 Boudreau and Ruddick, 1991). Inversely determined values predominantly fall within the 0.1-0.2 interval (apparent 6 th -11 th order of reaction), although lower and higher values were also determined (Fig. 4a). They thus confirm the previously observed dominance of values between 0.1 and 0.2 (Boudreau and Ruddick, 1991) and suggest that is relatively constant across a 18 wide range of different depositional environments encountered in the contemporary ocean. These findings are also in 540 agreement with the parametrization of the empirically derived power law description of OM degradation (Middelburg, 1989) another continuum model, which applies a globally constant parameter of 0.16:

545
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: Thus, the parameter range of the RCM is in good agreement with the 0.16 value imposed in the power law model.
Nevertheless, our results also reveal exceptions from this overall pattern. Particularly high -values ( > 0.25) were determined for the Arabian Sea and Bering Sea. Similarly, Boudreau and Ruddick (1991) also report exceptionally high values ( ≈ 1) for the Peru Margin and North Philippine Sea but attribute them to non-steady state depositional conditions. In our study, > 0.25 555 also coincides in part with the poorest data-model fits, such that the robustness of such values is uncertain. Particularly in the Bering Sea, the high -value (Table 4)   . Therefore, the somewhat atypical high -values determined 565 here (together with comparably low -values; see below) for the Arabian Sea transect are realistic and indicate that those sediments receive an OM mixture containing exceptionally large fractions of highly reactive compounds (Koho et al., 2013;Luff et al., 2000;Vandewiele et al., 2009), which are quickly consumed at the SWI (i.e., rapid decrease in reactivity with burial). This is further supported by the high PP rates in the photic zone during Monsoon periods and efficient vertical transport of diatom and haptophyte detritus (Barlow et al., 1999;Cowie, 2005;Latasa and Bidigare, 1998;Shalapyonok et al., 2001) 570 and deep carbon fixation by microorganisms (Lengger et al., 2019). Luff et al. (2000) showed that highly reactive OM in Arabian Sea sediments represents around 50-70% of the total OM fluxes at the SWI close to the continental margin, and it is composed mostly of fresh plankton, large aggregates, and faecal pellets. Further, the most reactive fraction is only present in the few uppermost millimetres of the sediment column. Thus, anoxic water conditions, associated with enhanced vertical OM fluxes within the OMZ enhance the burial of more reactive OM (e.g., Maerz et al., 2019;Zonneveld et al., 2010). In summary, 575 our inversely determined -values, alongside those compiled globally (Fig. 4a) and in addition to the empirically fitted power law parameter, reinforce the similarity of values across a wide range of depositional environments, thereby providing an important first order constraint on model parametrization on the global scale.
Parameter determines both the overall apparent reactivity, as well as the shape of the reactivity decrease with depth. High 580 -values result in an overall low OM reactivity, but a slow reactivity decrease with depth. In contrast, low -values result in a higher reactivity, but also a more rapid decrease in reactivity with sediment depth . Parameter is often related to the degree of degradation Middelburg et al., 1993) and very fresh OM generally reveals very low a values ( = 3.1•10 −4 y) (Boudreau and Ruddick, 1991;Westrich and Berner, 1984), while the heavily sulfurized OM deposited during Cretaceous Oceanic Anoxic Event 2 in the strongly euxinic subtropical Atlantic Ocean is characterized by 585 very high values ( = 10 4 y) (Arndt et al., 2009). Inversely derived -values span several orders of magnitude (Table 4; Fig.   4b). They thus corroborate the high variability observed in previously published -values compilations (e.g., Arndt et al., 2013;Boudreau and Ruddick, 1991) and suggest that parameter exert the major control on the spatial variability in OM reactivity.
Inversely determined values also broadly agree with previously observed global patterns. Low values ( < 10 2 y; i.e. high apparent OM reactivity) are generally determined for depositional environments that are characterized by the rapid deposition 590 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 ( > 10 3 y). Between these endmembers, the inversely determined parameter reveals a large variability that is further explored in the context of the respective depositional environments in Sect. 4.1.4. However, despite the large total range of inversely derived -values which spans eleven orders of magnitude, most values fall into the range 10 0 < < 10 4 years, thus narrowing the most plausible 595 interval of -values. This is a valuable finding that can help better inform model parameterisation in data poor regions and timescales (e.g., Hülse et al., 2018).

OM reactivity distributions
The inversely determined parameters and not only provide insights into the initial apparent reactivity of the depositing 600 OM but also allow the assessment of its evolution during burial. Here, we calculate the probability distribution of OM over at the time of deposition on the SWI, as well as after burial at 10, 100, and 1,000 cmbsf (centimetres below seafloor). The distribution profiles (Fig. 5) reflect both and . Yet, it is important to note that changes in OM-distributions with ongoing burial are not solely controlled by parameters and . Sedimentation rates also exert an impact on the degree of degradation since controls the residence time of OM at a distinct depth (and redox) horizon (i.e. the burial age; Eq. 7). Finally, 605 20 degradation rates also depend on benthic OM contents (donor control) . They thus also dictate changes in OM-distribution with burial. Overall, high OM contents result in higher degradation rates. Therefore, the resulting changes in distribution reflect different exposure to heterotrophic degradation at similar burial depth. Since -values are broadly similar across depositional settings, the differences in distributions are mainly driven by the range of inversely determined -values.

610
The Arabian Sea region (Fig. 5 jl) is characterized by the deposition of highly reactive OM (highest -values; Table 4).
However, OM reactivity rapidly decreases with sediment depth (low -values). The OM reactivity distributions show a clear depletion of the most reactive OM compounds in the uppermost 10 cm, resulting in a marked decrease in the overall reactivity over this depth, i.e., the distributions shift towards lower values, in agreement with previous findings (Luff et al., 2000).
Below 100 cmbsf the most reactive compounds have been degraded, and the remaining bulk OM is dominated by less reactive 615 compounds ( ( ) < 10 −5 y −1 ). The inversely determined rapid decrease in OM reactivity can be explained with the often rapid deposition of predominantly marine OM (Cowie, 2005;Seiter et al., 2005), whose degradation during settling is limited by short residence times and/or locally anoxic conditions (Luff et al., 2000). The shallow, sub-oxic/anoxic sediments of the Northern European margin (Aarhus Bay - Fig. 5d; Arkona Basin - Fig. 5e) reveal a similar OM-distribution and depth profile. However, the loss of the most reactive fractions within the uppermost 10 cm is less pronounced. This slowed decrease 620 in OM reactivity with sediment depth could be due to the burial of terrestrial, yet reactive OM (Aquilina et al., 2010). In contrast, Skagerrak sediments (S10 - Fig. 5g; S13 - Fig. 5i) are characterised by the most homogeneous, albeit comparably unreactive OM mixture (high -values; Table 4). The OM-distributions reveal small changes from the SWI down to 1,000 cmbsf, reflecting a low overall reactivity, and thus low degradation rates. This somewhat surprising finding can be explained with the atypical character of the Skagerraka shallow environment that connects the brackish Baltic Sea with the saline 625 North Sea and is characterized by an estuarine type 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). As a consequence, more than 90% of the OM that is 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(Dauwe et al., , 1999b.

630
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-distributions within the uppermost 10 cm. At all sites, below 100 cmbsf -values have converged to < 10− 4 y− 1 , since ( ) >> (see Eq. 23), and thus dominates -values at great depths.
As such, only the least reactive fractions are preserved, leading to an overall low OM reactivity in deeper sediments that is largely independent of the initial apparent OM reactivity (0). Sediments from the Bering Sea (Fig. 5m) are a notable 635 exception. OM reactivity decreases gradually throughout the sediment column, reflecting both unusually high and values (Table 4). 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 21 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). The heterogeneous distribution of the Bering Sea OM is, to a certain extent, 640 comparable to that observed in the Arabian Sea region (Fig. 5 jl). However, there is a fundamental distinction in the downcore evolution of OM-distributions at these contrasting regions that is driven by the magnitude of -values (Table 4), as well as timescales (Table 1).

Emerging environmental patterns in OM reactivity 645
Inverse model results reveal that parameter exerts the main control on the global variability in apparent OM reactivity, but also highlight the large range of plausible values spanning several orders of magnitude (Fig. 2b, Table 4). Our ability to quantify and predict benthic OM dynamics strongly depends on a general framework that allows the constraint of parameter when observational data is scarce or unavailable (e.g. global, future and past model scenarios). The need for such a general framework has been recognized since the early days of RTM. Based on the rationale that the reactivity of OM settling onto 650 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 . That global compilation, however, combined 655 OM degradation model parameters derived from structurally different models and using observational data of different complexity, thus somewhat compromising their comparability. Therefore, we further explore these links combining our newly generated dataset with those compiled from previously published RCM parameters (Fig. 6).
In agreement with the previous global assessment ( Arndt et al., 2013), we do not find any significant correlation between the 660 reactivity parameters ( and (0)) and master characteristics of depositional environments, as proposed for − (Boudreau, 1997;Tromp et al., 1995), − (Boudreau, 1997;Müller and Suess, 1979), and − (Boudreau and Ruddick, 1991). These results, therefore, further confirm that the spatial variability in and thus apparent OM reactivity is not determined by one single characteristic of the depositional environment but is the result of a complex and dynamic interplay of multiple environmental controls. However, despite the large scatter, very broad trends that generally agree with previously proposed 665 relationships emerge (Fig. 6). Depositional environments that are characterised by low sedimentation rates, low OM fluxes, or great water depths generally display high -values (i.e. low apparent reactivity), whereas settings characterized by high sedimentation rates, high OM fluxes, or shallow water depth generally display lower -values (high apparent reactivity). While these trends are extremely weak on a global scale, they are relatively more pronounced within a given depositional environment (i.e. on a local/regional scale). Despite the lack of clear global relationships between apparent OM reactivity and one 670 22 environmental master variable, general but, not significant correlations between OM reactivity and single environmental controls can be identified within a given environment.
Inverse model results show that water depth seems to be a useful and easily accessible first-order proxy for the complex interplay of environmental controls on OM reactivity. At the global scale, inversely determined apparent OM reactivity for the 675 deep, well-oxygenated sites in the Bering Sea and Argentine Basin fall in the lower end of the reactivity spectrum. Here, the prolonged exposure to oxygen (Hartnett et al., 1998) and pelagic degradation (Cram et al., 2018;Maerz et al., 2019;Müller and Suess, 1979;Weber et al., 2016), as well as the continuous ageing during settling (Griffith et al., 2010), contribute to the determined low apparent OM reactivity in the Argentine Basin (e.g., Riedinger et al., 2005Riedinger et al., , 2014Riedinger et al., , 2017 and Bering Sea. In contrast, shallow coastal environments, such as the Northern European region, reveal generally higher apparent OM 680 reactivities. In addition to this weak global trend, regional -ℎ reactivity trends emerge. For instance, the Northern European region, inversely determined apparent OM reactivity generally decreases (i.e. increase in ,) with increasing water depth ( Fig.   6 b,e). A notable exception to this regional trend is the comparably low OM reactivity determined for the Severn estuary. The Severn estuary 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 not only mix along the gradient 685 but also interact with high loads of fine, suspended particular material. In addition, intense, tidally driven depositionresuspension cycles (Dyer, 1984;Jonas and Millward, 2010;Manning et al., 2010;Uncles, 2010) result in a continuous reexposure of benthic OM, which might cause the lower apparent reactivity of the buried OM. Similar to the Northern European region, the Rhone delta transect, also reveals a decrease in apparent OM reactivity with an increase in water depths from the proximal zone to the shelf (Fig. 6e). This general trend is supported by the previously observed decrease in OM quality 690 recorded by individual OM fractions (Pruski et al., 2015), as well as the ageing of bulk OM across the Rhone shelf (Cathalot et al., 2010). However, the weak trends and numerous exceptions at both global and regional level caution against the uncritical use of such simplistic relationships. Our findings suggest that further important additional controls often complicate and sometimes inhibit the use of water depth as a predictor for apparent OM reactivity. For instance, the Arabian Sea is a notable exception from the general global trend observed between OM reactivity and water depth due to the presence of the pronounced 695 OMZ (e.g. Bogus et al., 2012;Fischer et al., 2012). Here, reactivity is 2-3 orders of magnitude higher than those derived for other deep-sea sites. Reactivity is rather comparable to reactivities determined for shallow sites in the Northern European margin and the Rhone delta. Across the Arabian Sea OMZ transect, only a modest regional depth trend is observed, with a decrease in reactivity (i.e. increase in ) from the OMZ to increasing water depths and bottom water oxygen concentrations (Cowie, 2005;Luff et al., 2000). While distinct from the overall -ℎ discussed above, the Arabian Sea OMZ region exhibit 700 relatively similar patterns to other high productive, anoxic water column systems, e.g., Peru margin (Boudreau and Ruddick, 1991), Blake ridge (Marquardt et al., 2010;Wallmann et al., 2006), as well as 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.

705
Inversely determined apparent OM reactivities reveal that the use of sedimentation rates are poor proxies for the complex interplay of environmental controls on OM reactivity (Fig. 6 c,f). For instance, sediments from the Northern European margin experience similar depositional rates but reactivity spans several orders of magnitude, suggesting the weak control imposed by sedimentation rates on OM reactivity (e.g., Arndt et al., 2013). This could partly be due to strong contribution of remobilized, pre-degraded OM and mineral protection in the Skagerrak sediments (Mayer, 1994a(Mayer, , 1994b, as well as the higher 710 OM fluxes onto shallow, sub-oxic/anoxic sediments of the Northern European margin, like Aarhus Bay (Dale et al., 2008b) and Arkona Basin . Similarly, in the Arabian Sea the absence of relationships between reactivity and sedimentation rates is likely to be associated to high fluxes of unprocessed OM within the OMZ (Cowie, 2005;Luff et al., 2000) and enhanced OM mineral protection outside of the OMZ (Keil and Cowie, 1999). In contrast, across the Rhone seaward transect a decrease in sedimentation rates is followed by an increase in , and thus a decrease in apparent OM reactivity . At 715 the proximal zone, the sedimentary OM deposition is characterised by terrestrial sources and aquatic PP, which is highly reactive. In contrast, at the shelf area, the major source of OM is marine PP . Here, slow settling through the water column allows longer exposure of OM to heterotrophic degradation prior to burial, resulting in an overall decrease in OM reactivity (e.g., Middelburg et al., 1997), whereas at shallow depths with great sedimentation rates the burial of reactive OM is favoured. Additionally, there is evidence of strong OM-mineral association, and thus mineral protection in the Rhone 720 shelf sediments (Mayer, 1994a(Mayer, , 1994b, which results in decreased reactivity. Our findings thus suggest that sedimentation rates might be useful predictors of reactivity across strong environmental gradients like those generally observed in land-ocean transition systems (e.g., Epping et al., 2002). In other environments, additional and more dominant controls on OM reactivity blur the use of sedimentation rates in other environments.

Influence of OM reactivity on benthic-pelagic coupling
OM degradation is the engine behind the complex and dynamic interplay of diagenetic processes in marine sediments.
Consequently, the quantity and quality of the deposited OM exerts an important control on benthic-pelagic coupling processes with potential implications on global biogeochemical cycles and climate. While the influence of OM quantity on benthicpelagic exchange across a range of depositional environments has been previously assessed (Krumins et al., 2013;Soetaert et 730 al., 1996Soetaert et 730 al., , 1998Thullner et al., 2009), these assessments made simplified assumptions about the variability of apparent OM reactivity. Consequently, the role of reactivity on benthic-pelagic fluxes remains poorly quantified. Therefore, we further explore the influence of differences in OM reactivity on OM degradation rates and their respective metabolic pathways, AOM dynamics, as well as benthic fluxes across all study sites.

Control of OM reactivity on OM degradation rates and respiration pathways
Depth-integrated (0-10,000 cmbsf) OM degradation rates ( ; Eq. 24) mirror the large spatial variability in apparent OM reactivity (i.e. parameter ) across depositional environments. However, results show a strong correlation of with OM fluxes (Fig. 7a, r 2 = 0.97; p < 0.01; n = 14), thus confirming previous findings by Henrich (1992) who showed that high depthintegrated heterotrophic OM degradation rates strongly correlate with high (typically reactive OM). Together with the 740 fact that does not show any systematic relationship with (Fig. 6g), as previously suggested (Boudreau and Ruddick, 1991;Müller and Suess, 1979) and that apparent reactivity often converges towards similar values in deeper sediments (see Sect. 4.1.3 OM reactivity distributions), these findings indicate 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. The Arabian Sea region follows that general trend with high and high ; however, these high depth-integrated rates are driven by exceptionally high OM 745 reactivities (i.e. high ) rather than high . As previously discussed, the combination of high PP rates and the presence of a well-developed OMZ delivers and preserves fresh, reactive OM to and within the sediments, where intense heterotrophic degradation takes place, consuming up to 70% of the most reactive compounds in the uppermost few millimetres of the sediment (Cowie, 2005;Koho et al., 2013;Luff et al., 2000;Rixen et al., 2019;Vandewiele et al., 2009 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 depositional 755 environments (Rhone proximal zone, Skagerrak, Bering Sea, and Argentine Basin; Fig. 7b) due to their overall dynamic nature that favours the dominance of organoclastic sulfate reduction (Bowles et al., 2014;Jørgensen and Kasten, 2006;Thullner et al., 2009), and the length of sediment column (1,000 cmbsf) integrated for those estimates. The contribution of the aerobic OM degradation pathway further decreases with a decrease in OM apparent reactivity (i.e. increase in ), whereas the contribution of denitrification increases. This pattern can be explained by the fact that highly reactive (i.e. low ) OM is 760 predominantly degraded in the uppermost sediment layers, thus increasing the relative contribution of oxygen, if oxygen is the dominant electron acceptor at these sediment depths (Freitas et al., 2020;Glud, 2008). In contrast, a decrease in apparent reactivity leads to a higher burial of OM to sub-oxic or anoxic sediments (e.g., Meister et al., 2013). In addition, in this case, 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 reduced species 765 (e.g., ammonium and hydrogen sulfide) that diffuse up to shallower sediment layers (Glud, 2008). Here, we do not consider the degradation pathways mediated by metal oxides due to a lack of SWI data that would allow constraining boundary conditions of metal oxides. While this might result in a slightly overrepresentation of sulfate reduction (e.g., at the Skagerrak; 25 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 (Freitas et al., 2020;Thullner et al., 2009), albeit in continental margin 770 sediments iron hydroxide reduction may represents an important OM degradation pathway (Beckler et al., 2016).
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 characterised by high apparent OM reactivity ( < 10 2 y) and intermediate OM fluxes ( ~ 10 −4 -10 −3 mol C cm −2 775 y −1 ), e.g., the Arabian Sea region (e.g., Cowie, 2005;Luff et al., 2000), the Arkona Basin and Aarhus Bay (e.g., Dale et al., 2008b;Jørgensen et al., 2019a;Mogollón et al., 2012). Here, the rapid degradation of highly reactive OM compounds takes place in shallower sediment layers where sulfate can be efficiently replenished by exchange with bottom waters. In addition, high relative contributions of sulfate reduction are also simulated for environments that are characterized by a reduced supply of less reactive ( > 10 3 y) OM deposition fluxes , such as the Argentine Basin and Bering Sea sites (e.g., Hensen et al., 2003). 780 Here, the lower quantity ( < 10 −5 mol C cm −2 y −1 ) and quality ( < 10 −5 y −1 ) of OM supplied to the sediment allows for a deeper sulfate penetration (Fig. 3 m-n), thus sustaining sulfate reduction over a greater sediment volume (e.g., Bowles et al., 2014;Meister et al., 2013;Riedinger et al., 2005Riedinger et al., , 2014Riedinger et al., , 2017. In contrast, methanogenesis is an important degradation pathway in environments that are characterized by a high supply ( > 10 −3 -10 −2 mol C cm −2 y −1 ) of both highly reactive ( = 10 0 -10 1 y) OM (Rhone proximal zone and Helgoland Mud Area, North Sea), as well as less reactive ( > 10 3 y) OM (Skagerrak). 785 In the first case, fast burial of reactive OM drives the development of a shallow methanogenic zone in the sediment (shallow sulfate penetration; Fig. 3 b,f) (e.g. Aromokeye et al., 2020;Oni et al., 2015aOni et al., , 2015bRassmann et al., 2020), whereas in the latter case low reactive 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). As a result, the sulfate available at depth is mostly channelled to AOM 790 instead of organoclastic sulfate reduction (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., 2005Riedinger et al., , 2014Riedinger et al., , 2017. Our model results highlight that OM reactivity exerts a dominant control on the relative significance of OM degradation pathways. Additionally, OM deposition fluxes control the amounts of material delivered to the respective redox layers, as well 795 as the timescale of burial which influences the apparent OM reactivity at depth (Eq. 26). This has important implications for carbon burial on geological timescales since depositional regimes with fast burial rates favour the escape of OM from thermodynamically more favourable metabolic pathways to methanogenesis. In these settings, sulfate and oxygen (if available) are predominantly channelled to the re-oxidation of reduced species produced during OM breakdown, with important implications for authigenic mineral precipitations (e.g. Arndt et al., 2009;Bahr et al., 2009;Henkel et al., 2012;Nöthen and 800 Kasten, 2011;Pierre and Fouquet, 2007;Sun and Turchyn, 2014;Wehrmann et al., 2013), isotopic signatures (e.g. Ussler and Paull, 2008;Zhou et al., 2016) and benthic alkalinity fluxes (e.g. Łukawska-Matuszewska and Graca, 2018;Rassmann et al., 2020;Snyder et al., 2007;Soetaert et al., 2007).
The depth of the SMTZ can vary from cm (e.g. Dale et al., 2008a;Jørgensen et al., 2019a;Knab et al., 2008;Mogollón et al., 2012;Oni et al., 2015aOni et al., , 2015b to hundreds of meters (e.g. Arndt et al., 2009;Wehrmann et al., 2013) as a function of different 815 environmental controls such as, among others, 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. 8 a-c). Generally, we observe a deepening of the SMTZ with decreasing sedimentation rates, which agrees with global patterns (Egger et al., 2018). Not surprisingly, the controls on both AOM depths and overall rates emerge from a combination of environmental forcing. 820 Confirming previous findings (Egger et al., 2018;Regnier et al., 2011) and in line with the dominant influence of OM reactivity on benthic redox zonation, OM reactivity exerts a significant influence on the depth of the SMTZ (Fig. 8b) and the depthintegrated AOM rates (Fig. 8d). The lowest Σ and deepest SMTZ depths are determined for the low apparent OM reactivity sites in the deep-sea, well oxygenated bottom waters regions (Table 4). In those regions, low OM reactivity results 825 in low rates of heterotrophic degradation (Fig. 7a) and, thus, deep 4 2− penetration (Fig. 3 m-n) or no 4 2− depletion at all.
In contrast, sediments from the Helgoland Mud area, North Sea, are characterized by a shallow SMTZ and intense Σ (Oni et al., 2015a(Oni et al., , 2015b. High OM reactivity in combination with high depositional fluxes (Fig. 8 b-e) result in the complete depletion of TEAs in the upper sediment ( < 50 cm), and thus allows the establishment of methanogenesis at relatively shallow sediment depth. However, high OM reactivity not always results in shallow SMTZ and intense AOM rates. Associated with 830 moderate depositional fluxes, high apparent OM reactivity reduces the amount of reactive OM that reaches the methanogenic zone, and thus causes a deepening of the SMTZ and a decrease in Σ . Such a pattern can be observed at the shallow, suboxic/anoxic sediments from the Aarhus Bay (Flury et al., 2016) and Arkona Basin  sites. Here, high OM reactivity ( < 10 1 y) supports high rates of degradation close to the SWI (Fig. 7a) where 4 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;835 Regnier et al., 2011). For moderate depositional fluxes, associated with lower OM reactivity ( > 10 2 y) such as determined for the Skagerrak sites a shoaling of the SMTZ and higher Σ are simulated (Dale et al., 2008a;Knab et al., 2008). The pre-degraded OM buried in Skagerrak sediments (Aquilina et al., 2010) is predominantly degraded via methanogenesis (Fig.   7b). Our results thus show that OM reactivity exerts an important control on SMTZ depths and AOM rates. Depositional environments that are characterized by the deposition of OM of intermediate reactivity, as well as those that receive high fluxes 840 of highly reactive OM reveal a shallow SMTZ and high AOM rates, whereas AOM is substrate limited in environments that receive low amounts of highly reactive material or that are characterised by the deposition of less reactive OM.

OM reactivity controls on sediment-water interface fluxes patterns
Model results show that apparent OM reactivity exerts an important control on OM degradation and benthic redox zonation 845 and, consequently, nutrient recycling at the seafloor. Here, we estimate benthic-pelagic fluxes based on our model-derived concentration depth profiles of dissolved species ( 2 , 3 − , 4 2− , 4 + , and 4 3− ). We emphasise that these flux calculations are based on simulated depth profiles at steady-state and are thus only weakly affected by sampling limitations. Uncertainties in SWI concentrations and upper boundary conditions are induced by the common loss of the uppermost sediment layers during sampling with gravity and piston corers. Multi-corers and benthic landers are more suitable sampling devices for SWI 850 assessments since such methods usually recover the undisturbed uppermost sediments, allowing for a better quantification of concentration depth profiles near the SWI (e.g. Hensen et al., 2006). Despite of such limitations, our model-derived flux estimates across the SWI allow us to investigate the links between nutrient recycling and OM reactivity and to assess qualitative differences across sites.

855
We observe that OM reactivity broadly controls the spatial patterns of benthic-pelagic fluxes, as well as the relative importance of different transport processes (Fig. 9). Overall, molecular diffusion is the main transport pathway (Fig. 9 f-j). However, the relative importance of bioirrigation increases at environments characterized by low deposition of less reactive OM. This occurs because of relative the 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 (Middelburg et al., 1997). Advective 860 fluxes become an important transport mechanism driving downward sulfate fluxes, especially in rapidly accumulating sediments, as well as at sites characterized by low OM reactivity and moderate depositional rate.
Overall, a high apparent OM reactivity, and thus high rates of OM degradation (Fig. 7a), especially at shallow sediment depth, drive high uptake fluxes of dissolved TEA (Fig. 9 a-c). Sulfate uptake dominates the benthic TEA uptake, particularly at sites 865 where oxygen and nitrate are unavailable TEA (Fig. 9c), due to the dominance of sulfate reduction across sites. Sulfate uptake fluxes increase with an increase in apparent OM reactivity. Oxygen and nitrate fluxes follow a similar trend. For instance, at 28 the Rhone proximal zone unusually high sedimentation rates (Pastor et al., 2011), which delivers fresh and highly reactive OM (Cathalot et al., 2010;Pruski et al., 2015) support high and Σ (Fig. 7a; Table 4), thus promoting intense oxygen consumption at the SWI (Rassmann et al., 2016) that are driven by both the aerobic degradation of OM, as well as the re-870 oxidation of reduced species (see Table S1 for a reaction network review). Similarly, high degradation rates result in an increased nitrate demand for denitrification. Thus, the high input of highly reactive OM ( > 10 −2 y −1 ; Table 4) yields high rates of aerobic degradation of OM and denitrification, and consequently intense downward oxygen and nitrate fluxes.
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 875 to 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 characterised by deposition of less reactive OM ( < 10 −4 y −1 ; Skagerrak, Argentine Basin, and Bering Sea), benthic-pelagic fluxes of oxygen, nitrate, and sulfate are generally low. The low degradation rates in these regions result in low TEA consumption at the SWI. 880 Benthic fluxes exert a crucial role in nutrient recycling. In the Arabian Sea region, large phytoplankton blooms are associated to monsoon conditions, which result in upwelling of nutrient-rich bottom water (SW monsoon) and deepening of the mixedlayer (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., 885 1996;Lomstein et al., 1990;Wiltshire and Manly, 2004) and maintained by benthic nutrient fluxes. Additionally, bottom water upwelling and nutrient recycling is an important mechanism 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. 9 d-e). The largest fluxes are determined for the Arabian Sea region, as well as the shallow Aarhus Bay and Arkona Basin, where the deposition of highly reactive OM drives intense degradation close to the SWI (Table 4). In contrast, the lowest nutrient 890 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 (e.g., Colombo et al., 1997;Griffith et al., 2010;Hartnett et al., 1998;Hensen et al., 2000;Sun and Wakeham, 1994;Wakeham et al., 1997).

Conclusions
We developed and applied an inverse modelling approach to quantify apparent OM reactivity (i.e. parameters and ) based 895 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.

29
Our results provide important first-order constraints on model parameterization that can inform on the choice of parameter values in data poor areas, as well as global model applications and past and future model scenarios. They corroborate previous 900 findings Boudreau and Ruddick, 1991) that the RCM parameter is globally relatively constant ( = 0.1-0.2). Exceptionally high > 0.2 are associated to particular cases, like the deposition and burial of highly reactive phytoplankton debris in high productivity regions associated with well-established OMZs, e.g., the Arabian Sea and the Peruvian shelf. In contrast, we show that the RCM parameter can span several orders of magnitude at a global scale. Such a large variability is linked to a combination of multiple environmental drivers that control the quality of OM delivered to 905 sediments, as well as the timescale of settling and burial, and thus exposure to favourable oxic conditions and ageing. These findings suggest that the parameter exerts the main control on the variability of apparent OM reactivity, and that future modelling efforts to quantify OM reactivity on a global scale could thus be reduced to one main reactivity variable. Based on inverse model results, we narrowed the most plausible range of to 10 0 -10 4 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 -range. 910 In addition, inverse model results further confirm 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 , , and . They caution against the use of single environmental master variable such as water depth, sedimentation rate or organic matter deposition fluxes beyond the local scale. Instead, they call for a more holistic, interdisciplinary exploration of OM reactivity in its entire environmental context. Such integrative assessments are key to mechanistically understanding benthic carbon and nutrient 915 dynamics.
Model results reveal that OM deposition fluxes exert a first-order control on total depth-integrated rates of OM degradation, except at settings that receive exceptionally highly reactive OM fluxes (e.g., the Arabian Sea OMZ region), and OM degradation rates are most intense in the uppermost sediment depths ( < 50 cmbsf). Yet, our model results highlight the 920 important role of apparent OM reactivity in controlling the relative significance of different metabolic pathways, redox zonation of the sediment, the rate of AOM, as well as benthic TEA uptake and nutrient recycling fluxes. The fluxes of dissolved species across the SWI are strongly coupled to OM turnover, which is driven by reactivity. Particularly, the release of macronutrients (e.g., ammonium and phosphate) back to pelagic systems seems to be closely linked to reactivity patterns. The benthic-pelagic fluxes have a crucial impact on sustaining PP in regions characterised by the occurrence of spring and summer 925 phytoplankton blooms, e.g., the Northern European margin, the Arabian Sea, and the Bering Sea. Alongside TEA availability and burial fluxes, OM reactivity regulates oxidative pathways and the predominance of anaerobic degradation of OM even at shallow sediment depths. In particular, the balance between organoclastic sulfate reduction and methanogenesis, and consequently the development of AOM, is strongly influenced by OM reactivity patterns.

930
Our findings have direct implications on the constraints for OM burial and preservation over long timescales. OM buried deeper than 100 cmbsf is expected to be predominantly of extremely low reactivity, which progressively becomes 30 thermodynamically unfavourable as an energy source for microbial processing (e.g., Arndt et al., 2013;LaRowe et al., 2020b).
Thus, the overall low reactivity character of OM in deeply buried sediments results in carbon sequestration at geological timescales. The processes leading to the preservation of highly unreactive OM are mostly attributed to selective degradation 935 during settling and burial (Zonneveld et al., 2010) and mineral protection (Hemingway et al., 2019), and possibly sulfurization (Hülse et al., 2019;Sinninghe Damsté et al., 1988, 1998. Therefore, systematically understanding the processes driving OM burial and carbon removal from active pools (e.g., atmosphere, soils, pelagic systems, and surface sediments) can help us to better constrain uncertainties associated to perturbations to the carbon cycle and climate. Thus, our findings demonstrate that OM reactivity and its underlying environmental drivers are crucial components of carbon immobilisation over long/geological 940 timescales.
Altogether, our inverse model results not only provide first, important constraints on model parametrization and a first step towards a general framework for model parameterization, but also help elucidate how changes in environmental conditions can affect benthic biogeochemical cycling and the role of marine sediments as a long-term carbon sink, and consequently their 945 significance in the past, present and future feedback of the climate system.

Code/data availability
The model code is available at the GitHub repository (https://github.com/felipesalesdefreitas/BG_2020_435). The data is available on PANGAEA database (www.pangaea.de) where indicated, or at cited literature along the manuscript.

Competing interests 985
The authors declare no conflicts of interest.    Hydrogen sulfide molecular diffusion coefficient 2 cm 2 y -1 331.61 Van Cappellen and  Methane molecular diffusion coefficient 4( ) cm 2 y -1 263.93 Van Cappellen and