A new look at the multi-G model for organic carbon degradation in surface marine sediments for coupled benthic-pelagic simulations of the global ocean

The kinetics of particulate organic carbon (POC) mineralization in marine surface sediments is not well constrained. This creates considerable uncertainties when benthic processes are considered in global biogeochemical or Earth system circulation models to simulate climate-ocean interactions and biogeochemical tracer distributions in the ocean. In an attempt to improve our understanding of the rate and depth distribution of organic carbon mineralization in bioturbated (0 – 10 20 cm) sediments at the global scale, we parameterized a 1-D diagenetic model that simulates the mineralization of three discrete POC pools (a ‘multi-G’ model). The rate constants of the three reactive classes (highly reactive, reactive, refractory) are fixed and determined to be 70 yr -1 , 0.5 yr -1 , and ~0.001 yr -1 , respectively, based on the Martin curve model for pelagic POC degradation. In contrast to previous approaches, however, the reactivity of the organic material degraded in the seafloor is continuous with, and set by, the apparent reactivity of material sinking through the water column. Despite the 15 simplifications of describing POC remineralization using G-type approaches, the model is able to simulate a global database (185 stations) of benthic oxygen and nitrate fluxes across the sediment-water interface in addition to porewater oxygen and nitrate distributions and organic carbon burial efficiencies. It is further consistent with degradation experiments using fresh phytoplankton reported in a previous study. We propose that an important yet mostly overlooked consideration in upscaling approaches is the proportion of the relative reactive POC classes reaching the seafloor in addition to their reactivity. The 20 approach presented is applicable to both steady-state and non-steady state scenarios, and links POC degradation kinetics in sedimentary environments to water depth and the POC rain rate to the seafloor.


Introduction
Mineralization and burial of particulate organic carbon (POC) in marine sediments is a core component of the Earth's carbon cycle, helping to regulate atmospheric oxygen and carbon dioxide levels over glacial timescales and longer (Berner and Canfield, 1989;Mackenzie et al., 2004).On timescales relevant to society and nutrient oceanic residence times (10 3 yr), the recycling of POC in the seafloor leads to a loss of fixed N and P from the ocean, thus modulating basinwide nutrient inventories (Van Cappellen, 2003).Global biogeochemical box or circulation models (collectively termed here Earth system models, ESMs) that are designed to study e.g., climate driven environmental change or proxies in sediment archives, thus rely on knowledge of the fraction of POC reaching the seafloor that is mineralized, and the fraction that is buried.
Earth system models deal with the recycling of biogenic material in seafloor sediment in a variety of ways, from the very simple (ignoring it completely) to the more complex yet computationally expensive coupling with 1-D vertically resolved or layered diagenetic models (e.g., Munhoven, 2007).Intermediate approaches make use of computationally efficient transfer functions (vertically integrated models) to simulate benthic feedbacks of the N, P, Fe and O 2 cycles (Middelburg et al., 1997;Soetaert et al., 2000;Wallmann, 2010;Tschumi et al., 2011;Bohlen et al., 2012;Somes et al., 2013;Dale et al., 2015;Capet et al., 2016;Yang and Gruber, 2016).Nutrient feedbacks computed using these functions may treat POC degradation implicitly and be independent of one another.Whilst transfer functions undoubtedly offer a powerful and pragmatic approach, 1-D diagenetic models provide a Published by Copernicus Publications on behalf of the European Geosciences Union.more complete coupling between carbon and nutrient cycles and their fluxes across the seafloor (Soetaert et al., 2000).
In 1-D sediment models, the depth in the seafloor over which POC degrades is usually determined by the reactivity or "freshness" of the material and its net downward movement due to burial and bioturbation (mixing by animals).It is important to gain an improved understanding of this interplay because the mineralization length scale strongly determines the flux of redox-sensitive elements (e.g., oxygen, nitrate, phosphate, iron etc.) across the sediment surface (Stolpovsky et al., 2015).Yet, it is notoriously difficult to prescribe the reactivity of POC once it arrives at the seafloor due the complex interaction of many ecosystem properties (Bianchi et al., 2016).These include mineral interactions, geopolymerization, thermodynamic factors and the activity of micro-and macro-biota (e.g., de Leeuw and Largeau, 1993;Aller, 1994;Mayer, 1995;LaRowe and Van Cappellen, 2011).The net result is an apparent reactivity of POC that decreases with time, often parameterized using reactive continuum models that assume a continuous distribution of organic matter reactivity as opposed to discrete compound classes with fixed reactivity (Middelburg, 1989;Boudreau and Ruddick, 1991).Implementation of continuum POC dynamics in bioturbated sediments is possible, yet cumbersome (Dale et al., 2015).Recently, Stolpovsky et al. (2015) successfully parameterized a continuum model for bioturbated sediments as a function of the POC rain rate to the seafloor using a global database of in situ flux measurements.In that model, POC degradation rate responds instantaneously to changes in rain rate and is thus applicable under quasi-steady state conditions, i.e., in settings where changes in rain rate are slow compared to the mean lifetime of POC in sediments.
Benthic POC mineralization in some ESMs relies on constitutive equations that describe linear degradation kinetics of one or more sedimentary POC fractions, or "G", as a function of e.g., burial rate (reviewed by Arndt et al., 2013).An advantage of so-called multi-G models (more than one POC class; Jørgensen, 1978) is that by modeling POC content explicitly, they can be employed in situations with marked temporal variability in POC rain rate.Whilst multi-G models work well for individual sites where POC degradation can be constrained against empirical data, there is no evidence that the transferability of these models is guaranteed at the global scale.In fact, building upon previous work (Toth and Lerman, 1977;Emerson et al., 1985;Boudreau and Ruddick, 1991;Tromp et al., 1995;Boudreau, 1997), Arndt et al. (2013) showed that no statistically significant relationship exists between the rate constants for POC degradation and major controlling factors such as burial rate and water depth.This severely limits the predictive capacity of ESMs to faithfully reproduce, and draw conclusions about, benthic feedbacks and the role of sediments in global biogeochemical cycles of carbon, nutrients and other non-conservative tracers.
In this study, we present a novel methodology to constrain the rate constants of discrete organic matter reactive classes in sediments using the multi-G approach.Our motivation is to provide a more realistic POC degradation model that can be coupled to ESMs.We employ the same global database of benthic O 2 and NO − 3 fluxes to constrain the model that we previously used to derive the continuum model (Stolpovsky et al., 2015).The model is further consistent with global POC burial rates and porewater O 2 and NO − 3 distributions.In contrast to all previous approaches (to the best of our knowledge), the apparent reactivity of the organic material degraded in the seafloor described here is continuous with, and set by, the apparent reactivity of material sinking through the water column.This is achieved by assuming that the efficiency of carbon transfer to the ocean interior and to the seafloor (i.e., the biological pump) can be approximated using a power law, that is, the Martin curve (Martin et al., 1987).We find that the relative proportion of the reactive classes reaching the seafloor allows for transferability, upscaling and extrapolation of the model to the global scale.The model also differs from our previous approach since it can more realistically deal with changes in carbon rain rate (Stolpovsky et al., 2015).The coherence between POC degradation in the water column and sediments should improve predictions of benthic-pelagic coupling of redox sensitive elements (e.g., O 2 , NO − 3 , PO 3− 4 , Fe 2+ ) in ESMs that explicitly account for seafloor mineralization processes.

Model structure
In this section, we first describe how the apparent reactivity (k app in yr −1 ) of POC sinking through the water column is calculated.This information is then used to define the rate constants and fractions of three reactive POC pools in bioturbated surface sediments (upper 20 cm) by simulating global databases of O 2 and NO − 3 fluxes using a reaction-transport sediment model.This latter model is modified from our previous study where POC degradation rates were defined using a reactive continuum approach (Stolpovsky et al., 2015).

Apparent reactivity of organic matter sinking through the water column
The apparent reactivity of POC raining to the seafloor was calculated by first considering the sinking flux of POC.This was achieved by simulating the Martin curve (Martin et al., 1987) that has been widely used to describe the rain rate of POC to the sediment (RRPOC), as a function of water depth (wd, m): where RRPOC(100) is the flux at 100 m and the flux attenuation coefficient of b = 0.86 is considered to be a globally averaged value (we shall return to this point later).We fitted an analytical model for POC advection (sinking) and degradation to this power function to derive the apparent reactivity of organic matter raining to the seafloor as a function of water depth, k app (wd).Three POC reactivity classes, i (highly reactive, reactive and poorly reactive), and two particle size classes, j (small and large) characterized by different sinking speeds were included.We considered that a minimum of three reactivity classes was required to be able to simulate the Martin curve accurately.
The degradation of each POC fraction in the water column was defined as follows: where k i,j represents a first-order degradation constant and w j is the sinking speed.The model was solved for wd ≥ 100 m with the following boundary conditions: where POC i,j (100) is the concentration of each reactivity/size class at 100 m water depth, POC tot ( 100) is the total POC concentration at 100 m, and f i and f j represent the reactivity and size class fractions, respectively.The contribution of large particles to total POC at 100 m was set to 20 % (f large = 0.2) based on results from the global 3-D biogeochemical model NEMO-PISCES (Aumont et al., 2017).This model was constrained against a global database of particle size distributions and dissolved tracers including iron, phosphate, silicate, nitrate and ammonium.Large and small particles were assigned a sinking speed of 50 and 2 m day −1 , respectively, based on the NEMO-PISCES model.The nominal cutoff size between these two size classes was assumed to be 100 µm (Aumont et al., 2015).These simplifications are necessary, considering that sinking speeds range over several orders of magnitude from 10 1 to 10 3 m day −1 (Kriest and Oschlies, 2008).Moreover, particles can undergo a variety of aggregation, disaggregation and repackaging processes during settling.The uncertainties in our model created by these additional factors will be discussed later.
In order to reduce the number of adjustable parameters, the attenuation of small and large organic particles was the same.The rationale for this particular choice is based on (i) the assumption that small particles originate at the top of the water column and not from disintegration of large particles at depth, and on (ii) observations that the fraction of POC bound in the fine particles to total POC concentration in the water column remains at ∼ 20 % below ca.1000 m (Aumont et al., 2017).To maintain this proportion, POC in small particles must be less reactive than POC in large particles, otherwise the reduced sinking speed of small particles would lead to a systematic increase in the f large /f small ratio with water depth.Therefore, the mineralization constants of small particles were normalized to the ratio of sinking speeds: With this assumption, the reduced sinking speed of small particles does not affect the POC concentration profile.The solution of Eq. ( 2) has the following form: The apparent reactivity, k app (wd), is then equal to The model (Eq.5) was fit to the Martin curve for water depths ≥100 m by adjusting the three decay constants (k large, highly_reactive , k large, reactive and k large, poorly_reactive ) and the relative contribution of the two reactive POC fractions (f highly_reactive and f poorly_reactive ).The reactive fraction was calculated as 1 − f highly_reactive -f poorly_reactive .The five unknown parameter values were constrained using Monte Carlo analysis by varying k large, highly_reactive from 50 to 200 yr −1 , k large, reactive from 5 to 50 yr −1 and k large, poorly_reactive from 0.1 to 5 yr −1 .These ranges span values determined for fresh plankton detritus sinking through the water column (McDonnell et al., 2015).The values for f highly_reactive and f poorly_reactive were 0.1 to 0.8 and 0.1 to 0.4, respectively, based on observations that the bulk of fresh POC is labile (Arndt et al., 2013;Belcher et al., 2016).A total of 127 combinations of k i,j and f i that provided a fit to the Martin curve (to within a nominal R 2 of 0.98) is shown as the cloud of red lines in Fig. 1, which were then averaged to determine k app (wd).The continuous decrease in k app (wd) with depth reflects the gradual consumption of labile organic matter moieties in the sinking material and the increasing proportion of poorly reactive material (Wakeham et al., 1997;Dauwe et al., 1999).We found that the ability of the model to simulate the Martin curve with only one particle size was very much reduced, with k app (wd) being very sensitive to its prescribed sinking speed (result not shown).Consideration of at least two particle types with different sinking speeds seems to be required to provide robust estimates of k app (wd).As a first approximation, the simulated k app (wd) can be approximated with the following power law (Fig. 1): This function is the basis for determining the rate constants of POC degradation in the sediment model, as described in the following section.
Figure 1.Apparent reactivity of total POC raining to the seafloor (degraded and buried fractions), k app (wd).The cloud of red lines shows k app (wd) from 127 combinations of values for k i,j that resulted in a fit to the Martin curve to within R 2 = 0.98.The solid line shows the power function fit to these combinations (Eq.8).The dashed curve shows Eq. ( 14) which corrects k app (wd) for the buried fraction from which k 0 and k 1 are estimated (see text).

Multi-G model of organic matter degradation in the sediment
In this section, we describe how the kinetics of POC degradation in bioturbated surface sediments (upper 20 cm) were calculated using a 1-D reaction transport model (Jørgensen, 1978;Berner, 1980).Benthic POC degradation was linked to POC reactivity in the water column.k app (wd) (Eq.8) provides the sediment model with the reactivity of bulk POC at the sediment surface, from which the rate constants of discrete POC reactivity fractions, i (highly reactive, reactive and refractory) were estimated.No distinction of size fractions was considered in the sediment model since all POC classes are buried and mixed through the sediment using identical transport coefficients.The rate of POC degradation at the sediment surface, RPOC(0, t), was calculated as follows: where POC tot (0,t) (in dry wt %) is the total POC content at the sediment-water interface and k app (wd) (in yr −1 ) is its apparent reactivity, derived previously.POC tot (0, t) is calculated directly by the model as a function of the rain rate of organic matter to the seafloor (see Eq. 19).In this way, the mean reactivity of POC at the sediment-water interface is continuous with k app (wd) for POC raining out of the water column for any given water depth.k app (wd) is simply treated as a boundary condition for the sediment model without mak-ing any mechanistic assumptions regarding, for example, microbial community structure and activity and POC preservation at the transition between the bottom water and the surface sediment.The nature of this transition is not well understood, but is assumed to be ultimately determined by the flux of POC deposited to the sediment.The rate of POC degradation with depth in the sediment, x, is where k 0,1,2 (in yr −1 ) represent first-order rate constants and POC 0,1,2 denote highly reactive, reactive and refractory carbon, respectively.In keeping with previous studies, k 0 refers to highly reactive POC that is degraded at the sediment water interface (Boudreau, 1997;Arndt et al., 2013).Degradation rate constants characterizing these fractions do not have to be the same as those in the water column because organic material that is buried and mixed into the sediment column becomes subject to a multitude of preservation effects induced by mineralogical interactions that reduce its overall bioavailability (Hedges and Keil, 1995).Although the firstorder constants in our approach are fixed, bulk POC reactivity decreases through the sediment due to diminishing contents of the more labile fractions, which approximates a continuous decrease of organic matter reactivity with sediment depth (Middelburg, 1989).
From Eq. ( 9), the k i values are related to k app (wd) via the relative abundance of each reactive class: (11) where f 0 , f 1 and f 2 are the dimensionless fractions of highly reactive, reactive and refractory POC at the sediment-water interface, respectively.f 2 is considered to be largely buried below the bioturbated zone and equal to the carbon burial efficiency (CBE) (Eq.21 below).This fraction does not communicate with the ocean on ∼ 10 3 yr timescales.From the outset, therefore, the model is consistent with global rates of POC burial, essentially leaving two reactive fractions to simulate diagenesis within the bioturbated layer.
The next step is to estimate k 0 and k 1 by first solving Eq. ( 11) for the ratio f 1 /f 0 as a function of k app (wd) and CBE: Obviously, f 1 /f 0 should be positive, which is the case when the numerator and denominator of Eq. ( 12) are positive (if they are both negative, then the ratio f 1 /f 0 is still positive, but k 0 < k 1 ): Solving these two inequalities gives an expression that provides limits on k 0 and k 1 according to k app (wd) and CBE: The rate constant of the refractory POC class (k 2 = 0.001 yr −1 ) was taken from the range of values of the least reactive fractions in 3-G model studies of surface sediments from around the world (Arndt et al., 2013).Since it is largely buried below the bioturbated zone it is classed as "refractory" for the timescales of interest here.Of course, POC remineralization does not stop at the base of the modeled sediment column, but this "deep" POC turnover has little influence on the flux of solutes at the lower model boundary as explained below.However, for the deep sea where sedimentation is very low, the residence time of POC in the bioturbated zone is long enough to allow considerable degradation of the refractory fraction.To avoid this, k 2 was scaled to the sedimentation rate: where ω acc (wd) is the sedimentation rate at a given water depth and ω acc (0) is the sedimentation rate at zero water depth calculated using an empirical function (Burwicz et al., 2011).
The function in Eq. ( 14) is plotted in Fig. 1 (dashed line) alongside k app (wd).From the end-members of this curve, k 0 has to be greater than 43 yr −1 whilst k 1 must be lower than k app (wd) in the deep sea.As baseline values in the model, we used 70 yr −1 for k 0 and 0.5 yr −1 for k 1 , noting that this k 0 value is similar to that determined for shelf sediments based on O 2 microelectrode measurements (76 yr −1 , Berg et al., 2003).The corresponding lifetimes (1/k) of the three POC fractions are ∼ 1 week, 2 years and 1000 years.The sensitivity of the model to these constants will be explored later.With knowledge of the rate constants and the CBE (Eq.21), f 1 /f 0 was calculated from Eq. ( 12).
The derived rate constants were used to simulate the degradation and content of POC (3-G model) in bioturbated surface sediments (upper 20 cm) using the 1-D reaction transport model.Solutes included O 2 , NO − 3 , NO − 2 , and NH + 4 .The reaction network is similar to that described in Stolpovsky et al. (2015) and includes the major reactions of oxygen and nitrogen in surface sediments such as heterotrophic denitrification, nitrification and anammox.Mineralization of POC is linked to aerobic respiration, nitrate and nitrite reduction, and anaerobic respiration.The latter produces oxygen demand units (ODU) that lump together reduced compounds such as sulfide, dissolved iron and manganese (Soetaert et al., 1996).The reaction network and model parameters are listed in Tables 1 and 2 in Stolpovsky et al. (2015).
POC in the sediment model was transported by burial with compaction and mixing by bioturbation.Solutes were transported by burial, molecular diffusion, bioturbation and by the non-local transport pathway of bioirrigation which represents flushing of surface sediments by burrow-dwelling animals.Mass balance equations were used to solve the concentration changes with time until a steady state was reached.For each POC fraction (in wt %), the relevant equation is where t (yr) is time, x (cm) is depth below the sedimentwater interface, φ (dimensionless) is porosity, φ(L) is porosity in compacted sediments, D B (cm 2 yr −1 ) is the bioturbation coefficient and ω acc (wd) (cm yr −1 ) is the sedimentation rate as a function of water depth (Burwicz et al., 2011).
For solutes (C i (x) in mmol cm −3 of pore fluid) where D Si (cm 2 yr −1 ) is the tortuosity-and bioturbationcorrected molecular diffusion coefficient of species i, α i (yr −1 ) is the bioirrigation coefficient, and R i is the sum of biogeochemical reactions affecting C i .Constitutive equations for φ, D B , D S , ω acc and α are given in Stolpovsky et al. (2015).
The unit conversion factor between POC and solutes due to mineralization reactions is as follows: where ρ (2.5 g cm −3 ) is the dry sediment density.Benthic fluxes of O 2 and NO − 3 across the sediment surface were simulated and compared to a global database of n = 185 observations.Porewater concentrations were further used for "ground truthing" the model, although far fewer data were available.
Fixed concentrations were imposed for solutes (Dirichlet boundary) at the sediment surface (x = 0 cm).Measured bottom water concentrations were used for O 2 , NO − 3 and NH + 4 whereas NO − 2 and ODU were set to zero since they do not usually accumulate in oceanic bottom waters.Flux continuity at the sediment surface for POC is equal to The upper boundary condition of each POC fraction was defined as a fraction of the rain rate: where f 1 and f 0 were derived from Eqs. ( 11) and ( 12) and f 2 is equal to the CBE at each station.RRPOC was derived from the depth-integrated rate of POC degradation in the bioturbated layer, RPOC B (mmol m −2 day −1 ), and the CBE after solving the following equations: where CBE is calculated according to an empirical function that depends on rain rate (Dunne et al., 2007).RPOC B was approximated from a mass balance of the measured benthic fluxes of O 2 , NO − 3 and NH + 4 (Stolpovsky et al., 2015).At the bottom of the bioturbated layer (x = 20 cm), a zero gradient (Neumann) boundary was applied for all species.In reality, deeper anaerobic mineralization will lead to nonzero concentration gradients for ODU and NH + 4 at 20 cm.Yet, the overall contribution of anaerobic degradation below 20 cm is a small fraction of total degradation.Log-log relationships of radio-tracer measurements of sulfate reduction versus depth in continental marginal settings show that over 95 % of organic matter (OM) is mineralized in the upper 20 cm (Holmkvist et al., 2011).Imposing a zero-gradient condition for these species at the lower boundary therefore introduces only a small error in the model results.Further refinements to the model could be made by extending the approach to consider mineralization kinetics in deeper sediments.This is beyond the scope of this study, which focuses on POC mineralization in surface sediments only.

Results
The rain rate of POC to the seafloor for each site in the database is compared with the Martin curve (b = 0.86) in Fig. 2a.The good agreement between the two demonstrates that the independent rain rate estimates based on benthic fluxes are largely consistent with the flux of sinking material based on open-ocean sediment trap observations.From the individual rain rates, CBE (Eq.21), k app (wd) and the three rate constants for POC degradation, the proportion of reactive-to-highly reactive POC was determined for all ocean depths (Fig. 2b).With increasing water depth, the fraction of reactive POC, characterized by k 1 , increases at the expense of the highly reactive fraction, characterized by k 0 .This results from the more rapid mineralization of labile components in sinking organic detritus.f 1 /f 0 increases from ∼ 0 on the shelf to around 250 in the deep sea, meaning that < 1 % of reactive detritus is highly reactive by the time it reaches the deep ocean.This is to be expected because highly reactive particles sinking with a speed of 50 m day −1 require many weeks to reach deep ocean sediments, which compares to the turnover time of POC 0 of only ∼ 1 week (k 0 = 70 yr −1 ).
Simulated O 2 fluxes for the global database fit very well to measured data using the imposed k 0 and k 1 of 70 and 0.5 yr −1 , respectively (red symbols, Fig. 3a).Since O 2 is the ultimate electron acceptor, either for the direct oxidation of POC or for oxidation of the reducing metabolites of anaerobic respiration (i.e., ODUs), this goodness-of-fit demonstrates that the total carbon degradation rate in addition to the CBE are correctly simulated at each station.The NO − 3 result shows more scatter (Fig. 3b).The benthic N cycle is convoluted with multiple sources and sinks extending over a large range of redox potential (Table 1 in Stolpovsky et al., 2015).Furthermore, the reaction rate constants of the N cycle are unknown at the global scale.At some stations, the model strongly underestimates the observed NO − 3 flux.For example, the green ellipse corresponds to the Peruvian margin where NO − 3 uptake is dominated by nitrate-storing sulfur oxidizing bacteria such as Thioploca spp.(Fossing et al., 1995;Dale et al., 2016).The current model does not consider this process.3 uptake by sulfur oxidizing bacteria (Dale et al., 2016).The R 2 values for baseline fit are provided on each plot (p 0.05).
Both O 2 and NO − 3 fluxes are not particularly sensitive to 10-20 % variations in k 0 and k 1 (colored symbols on the plots).Although the NO − 3 fluxes are more sensitive than O 2 , no combination of these parameters is able to improve a fit at all the sites simultaneously, which can be put down to struc-tural deficiencies in the model.Even so, the measured fluxes are reproducible with the model to within ca.20 % of the prescribed variations in the rate constants.We take this to mean that the k app (wd) provides a robust constraint on the mineral- 3 concentrations in the bioturbated layer for marine settings with different water depths, POC rain rates and sedimentation rates (listed above the plots).Black curves show the model predictions with the power function (Stolpovsky et al., 2015), red curves (this study) denote the 3-G model, blue curves show results using 2-G approach, and the symbols are the measured data.The numbers on each panel indicate the corresponding fluxes in mmol cm −2 yr −1 .References for stations (see Table S1 in Supplement in Stolpovsky et al., 2015): St. NH14A = Washington margin (Devol and Christensen, 1993), St. H = Arctic shelf (Devol et al., 1997), St. 241 = Mauritanian margin (Dale et al., 2014), St. WE206 = Washington margin (Hartnett and Devol, 2003), St. PS2361-1 = Southern Ocean (Smetacek et al., 1997).
ization constants and the benthic fluxes, given the limits set by Eq. ( 14).
As a further validation of the approach, modeled and measured vertical geochemical profiles were compared for the same data set as used to verify the power law model (Stolpovsky et al., 2015).These sites include diverse settings, such as the continental shelf (41 and 114 m water depth), upper slope (241 and 1025 m) and the deep sea (3073 m) (Fig. 4).Two of the sites are located in highnitrate/low-oxygen (HNLO) areas where poorly oxygenated waters bathe the seafloor (Washington and Mauritanian margin).The 3-G model (red curves) captures the observed trends in O 2 profiles through the bioturbated layer at these sites and also reproduces the trends in NO − 3 concentrations.The prominent subsurface peaks in modeled O 2 and NO − 3 at St. NH14A and St. H are apparently absent in the field data and probably caused by too intense bioirrigation in the model.The rate of this process is described using an empirical function based on total oxygen uptake and shows considerable uncertainty when applied globally (Meile and Van Cappellen, 2003).Nevertheless, these peaks do not affect O 2 and NO − 3 gradients at the sediment surface, and the benthic fluxes are simulated accurately.Overall, the fluxes from the 3-G model are close to those from the power law model, which provides a continuous and arguably more realistic description of POC degradation.A closer correspondence with the 3-G versus 2-G model with power law is expected, because the multi-G model becomes more "continuous" as the number of POC fractions increases.We would thus expect a model with even more fractions to perform better than the 3-G model, but with the added cost of more parameter uncertainty.
As a test of our assumption that three reactive size classes are required for simulating the data, we can compare our result with a 2-G model that includes only one reactive fraction.In this case, the rate constant for the single reactive fraction, k * , can be obtained from the following equation: where k 2 is defined in the same way as in the 3-G model and values of k * are intermediate between k 0 and k 1 .The benthic fluxes derived with the 2-G model compare well with those derived with the 3-G model, although higher-end NO − 3 fluxes tend to be overestimated relative to the 3-G model (Fig. 5).However, the porewater profiles using the 2-G model are very poorly simulated, especially for NO − 3 , due to the lack of carbon mineralization below the oxic layer where denitrification and anammox take place (Fig. 4).The 2-G model is obvi- ously structurally inferior to the 3-G model for simulating N sources and sinks accurately.
An alternative way to visualize this result is to compare the simulated NO − 3 penetration depth (NPD) for an additional set of 12 stations that was used to validate the 1-D diagenetic model "Muds" (Archer et al., 2002).NPD is defined as the sediment depth where NO − 3 concentration falls to 2 % of the local bottom water level.The NPD is a useful comparative metric because the depth-profile of organic matter degradation largely determines the depth where NO − 3 is depleted.The 3-G model is able to predict the NPD at 10 from 12 stations (exceptions being no.16 and 32) to within a factor of two or 1.5 cm, whereas for the 2-G model this decreases to five stations (Fig. 6).As before, the 3-G model is also closer to the power model result.This again demonstrates that one effective reactive POC pool defined by a single rate constant is insufficient to represent the reactivity of natural organic matter.At a minimum, two reactive fractions undergoing degradation on the timescale of burial through the bioturbated layer, plus a poorly reactive fraction that is largely buried, are required to simulate the vertical structure of O 2 and NO − 3 diagenesis in addition to POC burial.

Sedimentary organic carbon degradation constants
As recently pointed out by Keil (2017), understanding the mechanisms and processes that determine the flux of organic carbon in the ocean interior (e.g., the biological pump) and preservation within sediments are among today's great oceanographic challenges (Doney and Karnauskas, 2014;Heinze et al., 2015).Inclusion of organic carbon burial improves model distributions of redox sensitive tracers in the ocean such as oxygen and nutrients (Palastanga et al., 2011;Kriest and Oschlies, 2013) and, consequently, predictions of future climate scenarios.Proper parameterization of organic carbon degradation kinetics in sediments is thus of utmost importance for the forthcoming generation of ESMs that explicitly recycle and/or preserve sinking planktonic detritus in the benthic compartment.
In this study, we derived a 3-G diagenetic model that is able to simulate a global database of benthic O 2 and NO − 3 fluxes and that is further consistent with global POC rain and burial rates (based on CBE) and porewater O 2 and NO − 3 distributions.The opportunity to further ground truth the model using fluxes of other solutes is limited by the severe lack of in situ measurements.The model does not perfectly simulate field data from site-specific locations (Figs. 3 and 4), which is probably caused by the functions describing bioirrigation, bioturbation and CBE that themselves carry considerable uncertainties.Moreover, the model performs less well on the continental shelf due to the high sediment heterogeneity there and because it is mainly tuned to the N cycle and lumps together coupled Fe-Mn-S cycling into ODUs.POC degradation coupled to the reduction of metal oxides and sulfate can be important in bioturbated sediments on the continental shelf (Canfield et al., 1993).
A novel aspect of the model is that the reactivity of the organic material degraded in the seafloor is linked to the apparent reactivity of material sinking through the water column.It performs as well as our previous approach which used a continuous power law to describe POC mineralization kinetics at the global scale (Stolpovsky et al., 2015).That study showed that the apparent rate constant for POC mineralization can decrease by orders of magnitude through the bioturbated layer, and spatially throughout the ocean.This dynamic is reproduced with the 3-G model due to the dependence of the ratio of highly reactive versus reactive organic matter on rain rate (Fig. 2).Whereas the power model is less suitable for settings with strong temporal variability in POC rain rate (see Introduction), the current multi-G model is applicable in steady-state and non-steady state situations.
The motivation for this work is to provide a methodology for parametrizing POC degradation in ESMs of the global carbon cycle.In our approach, the corresponding rate constants of the highly reactive (k 0 = 70 yr −1 ), reactive (k 1 = 0.5 yr −1 ) and refractory (k 2 = 0.001 yr −1 ) classes are fixed.However, the relative proportion of these classes varies globally as a function of the apparent reactivity of sinking POC, k app (wd).To apply this in a global model, one needs to know k app (wd) (Eq.8).The k app (wd) relationship is dependent on our choice of the b coefficient in the Martin curve.Whilst not without its detractors (see below), a coefficient of 0.86 in the Martin curve seems to provide a robust constraint on reactivity of organic matter being degraded in the Figure 6.Observed and simulated nitrate penetration depths (NPD in cm) for sites compiled by Archer et al. (2002) that are independent to those in our database.The numbering on the x axis denotes the data set number given in Table 2 in Archer et al. (2002).The power model result is from Stolpovsky et al. (2015).To simulate the data, measured RRPOC and bottom water O 2 and NO − 3 listed in Archer et al. (2002) were used as boundary conditions.sediment (Fig. 2a).The model can be applied in ESMs that do not use the Martin curve as long as the k app (wd) of POC sinking to the seafloor can be recalculated.
The rate constants for k 0 and k 1 derived here correspond to POC that is degraded over a period of weeks to years.These fractions are thus mostly responsible for driving the coupling of sediment fluxes to water column processes in global biogeochemical models.Our proposed value of the highly reactive fraction k 0 (70 yr −1 ), equivalent to a lifetime (1/k) of around 1 week, is toward the high end of those previously reported for k 0 using "stand-alone" diagenetic models.In a review on this topic, Arndt et al. (2013) reported that values of k 0 in 3-G diagenetic models tuned to field data show a maximum of 10 1 yr −1 .The highest value reported was 76 yr −1 , derived for an Arctic fjord using a carefully constrained empirical model (Berg et al., 2003).This highly reactive pool was detectable from millimeter scale dissolved O 2 gradients obtained using microelectrodes.POC reactivity in most published benthic models is not constrained by finely resolved microelectrode data or in situ fluxes that constitute a robust constraint on the mineralization of the highly reactive POC fractions (Glud et al., 2009;Bohlen et al., 2011;Dale et al., 2016).The vertical resolution of porewater sampling using traditional extraction techniques will also not capture sharp gradients in solute concentrations at the sediment-water interface resulting from the degradation of highly reactive POC.Lower rate constants are typically derived by modeling the degradation of POC over larger sediment depths, manifested in the accumulation of ammonium or alkalinity for example (Berner, 1980).
Our rate constants are in agreement with the result of a well-documented two year-long decomposition experiment of fresh phytoplankton by Westrich and Berner (1984).They extracted a 3-G mineralization model from their data with rate constants of 24 ± 4, 1.4 ± 0.7 and 0 yr −1 (Fig. 7).Our values are able to simulate the temporal decrease in phytoplankton carbon with the same degree of accuracy despite the large differences in POC attenuation in the initial stages of the experiment (R 2 = 0.98 and standard error = 0.06 g L −1 in both cases).Whilst other incubation studies cited by Westrich and Berner (1984) provide similar rate constants, we were unable to simulate the NO − 3 porewater profiles in our database using the given values.This may be because the sampling interval in experimental studies does not capture the degradation of reactive POC with very short lifetimes.
It is further noteworthy that our k 0 value, and those listed in Arndt et al. (2013), tend to be orders of magnitude higher than those used in ESMs that explicitly couple sediment and pelagic carbon cycles.For instance, rate constants for aerobic respiration in the 1-D benthic model coupled to HAMOCC range from 0.005 to 0.01 yr −1 (Palastanga et al., 2011), whilst benthic organic matter mineralization in the MEDUSA model is associated with a turnover of 0.024 yr −1 (Munhoven, 2007).Sedimentary organic matter mineralization constants in ESMs are typically tuned to global organic carbon burial rates or sediment POC content (see Arndt et al., 2013).In general, POC content in surface sediments is a poor parameter to validate POC turnover because almost the entire organic matter raining to the deep sea floor is degraded, leaving behind refractory POC that is reworked with older POC as a result of bioturbation.Furthermore, surface POC content varies by a factor of 4-5 over most of the seafloor (0.5-2 wt.%), whereas TOU fluxes vary by 2 orders of magnitude (ca.0.1-10 mmol m −2 day −1 ).

Upscaling and uncertainties
The multi-G model for benthic POC mineralization proposed here compartmentalizes organic matter into well packaged macromolecular fractions that obey first-order kinetics with rate constants that are fixed and independent of processes taking place in the sediment.This is, of course, a gross oversimplification of reality.Sedimentary organic matter is a complex array of chemical components characterized by large divergences in degradability and further altered by mineral interactions, geopolymerization, thermodynamic factors and the activity of micro-and macro-biota (e.g., de Leeuw and Largeau, 1993;Aller, 1994;Mayer, 1995;LaRowe and Van Cappellen, 2011).These, plus other unknown and poorly understood factors (Burdige, 2007), lead to a continuous decrease in POC degradation over time, that is, with an apparent time-dependent rate constant (Middelburg et al., 1989).We are still far from a mechanistic understanding of the fundamental controls on organic matter degradation dynamics in diverse marine settings.Explicit formulations capture the structural complexity of natural organic matter are currently beyond the practical capacity of diagenetic models, and the multi-G approach remains the cornerstone of most diagenetic models published to date (Arndt et al., 2013).Clearly, then, the major controls on organic matter degradation and preservation can be encapsulated within the multi-G concept with considerable success.This is despite the fact that direct observations of POC turnover would tend to suggest that organic matter degradation kinetics are best described with a power law model (Jørgensen and Parkes, 2010;Flury et al., 2016).This implies that conceptually simple approaches like the multi-G model can be coupled to ESMs to account for sediment feedbacks in simulations of climate-relevant processes.The need for benthic-pelagic coupling in ESMs has driven the search for a generalized and transferable framework for predicting sedimentary POC degradation kinetics at the regional and global scale (Stolpovsky et al., 2015).Yet, a review of the literature (> 50 data sets) by Arndt et al. (2013), failed to unearth any statistically robust relationship between the rate constants for organic matter degradation (both multi-G and continuum models) and major controlling factors including RRPOC, sedimentation rate and water depth.In fact, it would appear that the first-order rate constants are completely independent of these drivers (Arndt et al., 2013).This is consistent with our findings, where the k values are fixed.Instead, we propose that an important consideration is the changing proportion of the reactive classes reaching the seafloor in addition to their absolute reactivity.The partitioning of sinking material into discrete reactive classes is linked to the apparent reactivity of the sinking material itself, which we approximated using the Martin curve.Thus, the degradation of organic carbon in our model is consistent with, and set by, carbon degradation in the water column.With this topdown constraint, global trends in sedimentary benthic fluxes and porewater chemistry can be reproduced.Stolpovsky et al. (2015) arrived at a similar conclusion, where the reactivity of organic matter described by a power law model was directly linked to the rain rate of organic material sinking out of the water column.Inclusion of highly reactive POC fractions in diagenetic models imposes limits on the structural complexity of verwww.biogeosciences.net/15/3391/2018/Biogeosciences, 15, 3391-3407, 2018 tically resolved benthic models coupled to global biogeochemical models.Diagenetic models are typically solved over hundreds of vertical sediment layers whose thickness increases from, for example, sub-millimeter scale at the sediment surface to e.g., decimeters or meters in deeper strata.By comparison, and mainly due to computational constraints, benthic models coupled to ESMs tend to have a relatively low number of layers (e.g., 11 layers over the upper 10 cm in HAMOCC; Palastanga et al., 2011).Consider a 1-D advection-diffusion-reaction problem for the firstorder decay of POC (k = 70 yr −1 ) in a bioturbated layer (L = 10 cm) on the lower continental slope (2000 m).For these settings, D B = 0.5 cm 2 yr −1 and ω acc = 0.005 cm yr −1 are typical (see references in Stolpovsky et al., 2015).The Péclet number is < 1 for these conditions, meaning that diffusive transport (bioturbation) dominates advective transport (burial).Furthermore, for highly reactive POC, the characteristic timescale for POC degradation is much smaller than for POC transport.To maintain model accuracy, therefore, the spatial discretization in the diagenetic model ( x) must then be smaller than the relevant length scale, (D B /k) 0.5 , which in this example is less than 1 mm.This limitation is maintained throughout the ocean for typical values of D B and ω acc , requiring a stable routine for solving the diagenetic model (Berg et al., 2007).With a lower k value of 0.01 yr −1 as in the ESM mentioned above, the maximum x increases significantly.These calculations are approximate since the dynamics of redox sensitive elements is typically non-linear, and here we only intend to demonstrate the need for higher spatial resolution to capture the fine-scaled distributions of sedimentary variables as k increases.This is an important consideration for many processes, particularly for carbonate dissolution by respiratory CO 2 which takes place in the upper millimeters of the sediment (Hales, 2003;Jahnke and Jahnke, 2004).In general, for highly reactive classes, large reductions in the spatial resolution of the diagenetic model will lead to inaccurate benthic feedbacks for redox sensitive elements such as O 2 , NO − 3 and PO 3− 4 .Given these caveats and the computational penalties demanded by online coupling of diagenetic models to ESMs, it is understandable that global or regional models favor computationally efficient empirical transfer functions to simulate benthic feedbacks of N, P, Fe and O 2 cycles (Middelburg et al., 1997;Soetaert et al., 2000;Bohlen et al., 2012;Somes et al., 2013;Dale et al., 2015;Wallmann, 2010;Capet et al., 2016;Yang and Gruber, 2016).
The application of the multi-G model relies on knowledge of k app (wd) (Eq.7) and the CBE (Eq.21).Both of these are imperfectly known at the global scale.With regards to k app (wd), the depth attenuation of carbon flux in the ocean is critical to estimate the apparent reactivity of organic carbon raining to the seafloor.The Martin curve has found wide application in global models (e.g., Najjar et al., 1992;Ito and Follows, 2005).However, concerns about the accuracy of the Martin curve in capturing small-scale spatial and temporal variability of the biological pump in the upper ocean have called into question its suitability for global scale applications (Buessler et al., 2007).A rather wide range of the b coefficient between 0.4 and 1.75 has been reported and discussed at length (Keil, 2017;Guidi et al., 2015 and references therein).The sinking POC flux also varies temporally as well as spatially in the open ocean (Henson et al., 2015).Certain components of the biological pump, such a mineral ballasting or diatom coagulation, can rapidly deliver fresh organic matter to the ocean interior (Buesseler et al., 2007;Belcher et al., 2016).As a first approximation, however, the Martin curve with a b = 0.86 captures most of the independent rain rates estimates in our database (Fig. 2a).These rain rates were derived from benthic fluxes that integrate POC degradation over much longer periods of time and therefore provide a more average impression of how carbon rain varies with water depth.This supports the use of k app (wd) derived from the Martin curve.A ±20 % change in b only leads to minor changes in the benthic flux.For instance, for St. NH14A (114 m) in Fig. 4, the flux of O 2 calculated with b = 0.69 differs from that calculated using b = 1.03 by 14% and NO − 3 by 13.5 %.For St. PS236-1 (3073 m), the change is only 0.02 and 0.13 % respectively.
Trends in CBEs are generally understood throughout the ocean: low CBEs characterize the deep sea (< 5 %) versus higher CBEs on the margins (> 20 %).Yet, data are limited and subject to error due to mismatches in the timescales of input parameters needed to determine CBEs (Burdige, 2007).Additionally, the mechanistic controls on CBEs are still not resolved (Hedges and Keil, 1995).Although the Dunne relationship given in Eq. ( 21) does not perfectly describe rain rates at the global scale, it is consistent with previous studies and empirical models based on global POC burial rates (e.g., Flögel et al., 2011;Kriest and Oschlies, 2013).The good correspondence between measured and observed O 2 fluxes (Fig. 3a) indicates that the CBE formulation by Dunne et al. (2007) is appropriate to use in global models.
Finally, an additional source of uncertainty rests with the sinking speed of particulate organic matter through the water column and the distribution of particle sizes, both of which were used to derive the apparent reactivity of settling POC (see Sect. 2.1).These were parameterized according to the NEMO-PISCES model which is tuned to biogeochemical tracers in the water column as well as global benthic O 2 fluxes (Aumont et al., 2017).Whilst it is beyond the scope of the present study to discuss the uncertainty in these model aspects, it worth noting that the contribution of large POC particles to total POC in NEMO-PISCES (assumed here to be 20 %) does vary, particularly on the margins.For instance, larger particles contribute up to 50 % of the total in eastern boundary systems and areas of higher primary production (Oliver Aumont, personal communication, 2017).Since larger particles will tend to sink more quickly, the flux of POC to the seafloor in these regions will be greater than predicted by the Martin curve (i.e., b < 0.86).However, parameterizing the sinking velocities is evidently not straightforward as they span several orders of magnitude (meters to thousands of meters per day) (Kriest and Oschlies, 2008).Even small changes in the sinking speed or remineralization rates can change the e-folding depth for carbon export by hundreds of meters (DeVries et al., 2012).Ecosystem structure, and the associated packaging, disaggregation and repackaging processes during settling may create more spatial and temporal variability in POC sinking flux, hence its apparent reactivity, than assumed by a single b coefficient (De La Rocha and Passow, 2007;Wilson et al., 2008;Henson et al., 2012).Moreover, lateral transport of POC on the continental margins to deep-sea depocenters can be accelerated in benthic and intermediate nepheloid layers and by mass wasting events (Jahnke et al., 1990;Thomsen, 1999;Inthorn et al., 2006).A quantitative accounting of the contribution of lateral POC transport to the biological pump and preservation in sediments must rank among the major outstanding challenges in ocean science.

Conclusion
Biogenic elements (C, N, P, Fe) that are packaged into organic material in the sunlit ocean are either mineralized in the ocean interior or deposited on the seafloor.Most organic material reaching the seafloor is mineralized and returned to the ocean.The rest is permanently lost from the bioavailable pool either by microbial activity (e.g., N loss to N 2 by denitrification) or through mineral precipitation and/or burial.Over long timescales, these sediment sinks begin to exert a major control on nutrient inventories, ocean fertility, and climate-ocean interactions.Realistic model simulations of ocean chemistry in past, present and future scenarios must account for the recycling of organic material at the seafloor.
In an attempt to improve our understanding of the rate and depth of organic carbon mineralization in seafloor surface sediments, we parameterized a 1-D diagenetic model for organic carbon mineralization that can be coupled to Earth system models.In contrast to previous approaches, the apparent reactivity of the organic material degraded in the seafloor is continuous with, and set by, the apparent reactivity of material sinking through the water column.We propose that an important and mostly overlooked consideration in previous upscaling approaches is the proportion of the reactive fractions classes reaching the seafloor in addition to their intrinsic reactivity.The results imply the presence of a highly reactive organic carbon class at the sediment water interface with a turnover time of days to weeks.Mineralization of this pool will have important implications for the flux of redox sensitive elements to/from the ocean as well as for carbonate dissolution by respiratory CO 2 .
Many unknowns of course remain and affect the predictive capabilities of our findings.Further studies of the dynamics of POC mineralization in the water column and its reactivity at the point of deposition on the seafloor are needed.Until the spatial heterogeneity of the vertical and horizontal flux of organic matter in the ocean is more completely understood and predictable, knowledge of the reactivity of organic matter raining to the seafloor will be limited to the very few sites worldwide where benthic studies have been undertaken.Diagenetic models and Earth system models are, by definition, simplifications of nature and will require simplified global relationships like the one proposed here for the foreseeable future.

Figure 2 .Figure 3 .
Figure2.(a) Organic carbon rain rates (RRPOC, Eq. 20) derived for each station in the database versus the Martin curve (Eq.1).The flux at 100 m used to derive the Martin curve is 12 mmol m −2 day −1 (average of 16 stations between from 90 to 114 m water depth).The benthic data are apparently consistent with a Martin coefficient of b = 0.86.(b) The dimensionless ratio f 1 /f 0 (Eq.12) describing the relative abundance of the reactive POC fraction in sediments (POC 1 ) versus the highly reactive fraction (POC 0 ).The sigmoidal distribution of the data arises from the dependency of CBE on organic carbon rain rate.

Figure 5 .
Figure 5. Normalized benthic O 2 and NO − 3 fluxes for the 3-G versus 2-G models.The solid line indicates the 1 : 1 correlation.