Harnessing Marine Sediments: The Dance of Sorption and Transformation in Organic Carbon Preservation

Methodology

Summary of the modeling approach

An RTM is constructed to account for the cycling and retention of DOC within marine sediment. This RTM is subsequently simulated using artificial intelligence, which serves as a mechanism for performing a comprehensive process importance analysis. The flowchart illustrating the model implementation can be found in Supplementary Fig. 3. The artificial intelligence method employed in this research is an ANN, and the RTM is simulated using ANN through a Monte Carlo strategy (Supplementary Fig. 3, stage 1), where the RTM is executed numerous times (for instance, 1,000–2,000) with randomly varied input parameters to yield the model outputs (for example, the preservation rates for various MOC fractions). These input-output datasets are then utilized to train an ANN (Supplementary Fig. 3, stage 2)³⁵, facilitating the evaluation of the sensitivities of each output to diverse input parameters categorized into numerous processes.

From the 68 unknown input parameters in the model, we determined the ranges and statistical distributions for 6 of these parameters, including water depth⁴⁵, sediment accumulation rate⁴⁶, sediment surface porosity⁴⁷, and sediment-water interface concentrations of POC, which we assume to be equal to total OC⁴⁷, NO₃ (refs. ^48,49) and O₂ (refs. ^48,49) (see Supplementary Tables 1 and 8 and Supplementary Figs. 1 and 2) derived from global gridded datasets. For other typical unspecified parameters, we gathered information from ten previous studies that had carried out reactive transport modeling based on field data (summarized in Supplementary Tables 1 and 2). The ranges for newly identified parameters were also sourced either from literature²⁴ or by adjusting our model to fit field-modeling data from existing research (Burdige et al.^14,50). Although we have validated the RTM against field datasets and selected ranges for some of the new model parameters by considering the calibrated parameter values (Supplementary Section 2.2), such model fittings to field datasets are beneficial but not obligatory. This is because the RTM, which is founded on established or theoretical concepts, inherently incorporates process-based knowledge in its design, and the additional implementation of ANN and Monte Carlo enables the use of parameter range distributions rather than fixed values. A more comprehensive explanation of the overall modeling procedure can be found in Supplementary Section 1.1.

Formulation of the RTM

The RTM accounts for all prevalent early diagenesis reactions for varying compounds including dissolved species (O₂, SO₄, NH₄, NO₃, dissolved inorganic carbon (DIC), H₂S, CH₄, Fe²⁺, Mn²⁺, DOC₁ to DOC_m, GPS₁ to GPS_p, and lrDOC, where the subscript m indicates the utmost quantity of DOC pools and subscript p signifies the utmost quantity of GPS pools) and particulate species (highly reactive iron oxide, Fe(OH)₃^HR, moderately reactive iron oxide, Fe(OH)₃^MR, non-reactive iron oxide, Fe(OH)₃^UR, MnO₂, FeS, FeS₂, S⁰, and POC₁ to POC_n, where subscript n represents the maximum quantity of POC pools) as detailed in Supplementary Table 10. Here, we considered seven species for POC, four for DOC, and two for GPS in addition to lrDOC. The selection of distinct carbon species is primarily based on consistency with existing literature²⁴ and consideration of various aspects of modeling, including alignment with the conceptual framework, smooth transitions of rates across different species for numerical solutions, and minimizing the quantity of unknown model parameters. It is important to note that, as illustrated in Supplementary Section 1.2.2, the decision regarding the number of these species, such as DOC pools, does not significantly influence the model outputs.

The three governing equations of the RTM for dissolved species, particulate species, and sorbed species are enumerated as follows; full details regarding model development and validation can be found in Supplementary Sections 1 and 2:

1. (1)

   The governing equation for dissolved species:

   $$begin{array}{l}left(varphi +{rho }_{{mathrm{s}}}varepsilon {{K}_{{mathrm{d}}}}_{i}right)frac{partial {{C}_{{mathrm{d}}}}_{i}}{partial t}=frac{partial }{partial z}left((varphi D+{rho }_{{mathrm{s}}}varepsilon {{K}_{{mathrm{d}}}}_{i}{D}_{{mathrm{b}}})frac{{partial {C}_{{mathrm{d}}}}_{i}}{partial z}right)qquadqquadqquadqquad-frac{partial }{partial z}left(left(varphi {v}_{{mathrm{d}}}+{rho }_{{mathrm{s}}}varepsilon {{K}_{{mathrm{d}}}}_{i}{v}_{{mathrm{p}}}right){{C}_{{mathrm{d}}}}_{i}right)+varphi alpha left({{C}_{{mathrm{d}}}}_{i}left(0right)-{{C}_{{mathrm{d}}}}_{i}left(Zright)right)qquadqquadqquadqquad+varphi sum _{j=1}{{R}_{{mathrm{d}}}}_{i,,j}-{{k}_{{{mathrm{sorp}}}}}_{i}{{C}_{{mathrm{d}}}}_{i}+frac{{{k}_{{{mathrm{sorp}}}}_{i}}{{{{{mathrm{Kd}}}}_{{{mathrm{sorp}}}}}_{i}}{{S}_{{mathrm{d}}}}_{i}end{array}$$

   (1)

2. (2)

   The governing equation for particulate species:

   $${rho }_{{mathrm{s}}}varepsilon frac{partial {{C}_{{mathrm{p}}}}_{i}}{partial {{t}}}={rho }_{{mathrm{s}}}frac{partial }{partial z}left(varepsilon {D}_{{mathrm{b}}}frac{{partial {C}_{{mathrm{p}}}}_{i}}{partial z}right)-{rho }_{{mathrm{s}}}frac{partial }{partial z}left(varepsilon {v}_{{mathrm{p}}}{{C}_{{mathrm{p}}}}_{i}right)+{rho }_{{mathrm{s}}}varepsilon {{R}_{{mathrm{p}}}}_{i}$$

   (2)

3. (3)

   The governing equation for mineral phase, MOC, resulting from the kinetically sorbed fraction of dissolved species:

where C_di refers to the concentration of dissolved entities i (mM or µmol cm⁻³ of pore liquid), C_pi is the concentration of particulate entities i (g g⁻¹), S_di pertains to the concentration of dissolved entities i that are kinetically bound to sediment minerals (µmol g⁻¹ of solid sediments), φ denotes porosity, ε represents the solid fraction of sediments, defined as 1 − φ, v_d and v_p indicate the burial velocities of pore fluid and particulate entities (cm yr⁻¹), ρ_s is the dry density of sediments (g cm⁻³), D_i is the effective diffusion coefficient of dissolved entities i (cm² yr⁻¹), α represents the bio-irrigation coefficient (cm² yr⁻¹), D_b is the bioturbation coefficient (cm² yr⁻¹), z stands for the sediment depth concerning the coordinate system situated at the sediment–water interface (cm), while R_p, R_d and R_Sd represent reaction rates of particulate, dissolved, and kinetically bound entities (yr⁻¹, µmol cm⁻³ yr⁻¹ and µmol g⁻¹ yr⁻¹, respectively), which are temporally and spatially variable. Additionally, k_sorp indicates the mass transfer rate between the dissolved and kinetically sorbed phases to minerals (MOC pools) (yr⁻¹), and Kd_sorp is the so-called distribution coefficient in the kinetic mass transfer formulation (cm³ g⁻¹).

The initial phase of hydrolysis is regarded as analogous to the established first-order multi-POC degradation framework referred to as the multi-G model⁵¹, involving a sequence of POC pools transitioning into a singular DOC pool, DOC₁, simultaneously:

where k_i signifies the hydrolysis rate constant, considered similarly to the degradation rate constants of POC in the continuum model, following prior research^9,34,52.

The succeeding phase of hydrolysis has been articulated using a consecutive first-order reaction expression^20,53:

where λ_DOCi denotes the conversion rate from DOC_i to DOC_i+1, and λ_DOCi−1 refers to the conversion rate from DOC_i−1 to DOC_i in yr⁻¹.

An identical mathematical representation is utilized to define geopolymerization^54,55 as detailed in Supplementary Section 1 alongside other particulars.

Calculation of PE

PE, also referred to as BE, traditionally assigned to POC^28,36 is specified as follows:

where L indicates a specified depth, here taken as 1 m. In the current analysis, due to a comprehensive examination of the fate of DOC in our model, we can present a more precise assessment of PE that takes into account the portion of solid phase OC that has undergone hydrolysis and attachment to minerals:

The sorption rate reflects the net DOC kinetic sorption rate (sorption rate minus desorption rate or the net MOC formation rate). Further insights concerning the rates are elaborated in Supplementary Sections 2.4 and 3. While considering depth-versus-age horizons in initial diagenesis modeling may be crucial for global forecasting, as emphasized recently²⁸, the present investigation focuses solely on a consistent depth horizon since the objective is to gain understanding into the mechanisms that influence OC preservation rather than making global estimates.

ANN for process significance examination

The ANN serves as a flexible and universal approach for function approximation challenges and is recognized for its use in intricate, nonlinear systems^35,56. The commonly adopted ANN structure consists of a three-tier configuration incorporating input, hidden, and output layers^35,57,58,59. Each of these layers is comprised of a series of nodes (neurons) where the quantity in the input and output layers correlates with the number of input and output variables, respectively. The quantity of neurons within hidden or intermediary layers should be optimized during the fitting process^57,60. The primary equation utilized for processing information (or signal) within the ANN structure is a straightforward algebraic equation in the form of y = w × x + b, applicable to each neuron in the hidden layer. The information is subsequently aggregated for all nodes; further functions termed transfer (or activation) functions are applied to the input and output information, as elaborated elsewhere^35,57. Here, x denotes the input information (or signal), y indicates output information, w represents weights, and b signifies biases. Weights and biases are the hyperparameters of the ANN, established post-fitting and, once defined, create an empirical network utilized for new forecasts. In this study’s scope, we employ ANN solely for process significance investigation, not for predicting outcomes. We utilize the partial derivative method^35,57,58,61 for process significance analysis. In summary, this method employs the derivatives of the equations within ANN’s structure to illustrate the contribution of each ANN input in influencing the ANN output; for instance, for the principal equation y = w× x + b, the derivative corresponds to w. Consequently, the w value for each neuron indicates the intensity of the signals that flow through that neuron^35,57,58,59. This indicates that during the process significance analysis, the values of input parameters do not have a significant influence; instead, it is their variations that are critical and manifest in the architecture of the ANN. Specifics of the ANN model employed here have been chosen based on ref. ³⁵ and are presented in Supplementary Section 1.5.

Validation of the Model

We authenticate our model through several approaches. We validate our established governing equations of the RTM based on an analytical method. In this method, for conditions where equilibrium adsorption and kinetic sorption are anticipated to behave similarly, i.e., at elevated exchange rates, we initially execute the model by switching off equilibrium adsorption and subsequently run the model again by disabling the kinetic sorption expression. We then compare the model outputs from these two types of simulations. Existing field-modelling data (primarily from Santa Barbara Basin, given the extensive dataset available)^14,62 are utilized as illustrated in Supplementary Figs. 5–11 and Supplementary Table 6. Additionally, we validated the model based on mass budgets. The employment of the ANN is corroborated by its capability to describe the data (measured using goodness-of-fit metrics detailed in Supplementary Section 2.2) and through the uncertainties it generates in the process significance analysis. Lastly, the validation of the complete RTM modelling procedure was conducted using mass budgets from the averaged model outputs over numerous runs of the Monte Carlo method. This was executed employing the concept of mass flow in the model as depicted in Supplementary Fig. 14.

The outcomes of the model output comparison between instances when kinetic sorption is active and when equilibrium adsorption is active at a high mass transfer rate reveal an exceptional alignment (R² = 1.000; Supplementary Fig. 4), affirming our methodology for the development and application of sorption formulation within the governing equations.

The outcomes of the model’s fit to various field or modelling datasets, including Meysman et al.⁶² (Supplementary Fig. 5), Kraal et al.⁶³ (Supplementary Figs. 6 and 8 and Supplementary Tables 3–5) and Burdige et al.^14,50 (Supplementary Figs. 9–11 and Supplementary Tables 6) exhibit excellent correspondences between our model and existing field-model datasets for the majority of concentration-versus-depth profiles. The exceptions typically include FeS₂ and Mn²⁺ profiles, which present a less satisfactory model fit owing to the absence of carbonate species in our model. The additional intricacy of our model (with more unknown parameters) is validated in eight phases against the Santa Barbara Basin dataset^14,50 (Supplementary Table 6), demonstrating that each phase is justifiable concerning enhancements in model fits to the data relative to the complexity cost, based on model selection criteria⁶⁴ increasing from 0.626 in phase 2 to 0.843 in phase 7. Burdige et al.^14,50 further assessed δ¹³C, Δ¹⁴C and carbon-to-nitrogen ratios in their model and compared these with field measurements not undertaken in this study.

The ANN model managed to fit the data successfully in all instances, achieving the best predictive fit Nash–Sutcliffe model efficiency criterion⁶⁵ ranging from 0.923 to 0.944 (Supplementary Table 7 and Supplementary Fig. 12). The uncertainties in the ANN process significance analysis, calculated as a 95% confidence interval, were relatively small (refer to error bars in Fig. 3). These uncertainties vary between 4.6% to 29.9% (averaging at 12.6%) of the mean values for the examined cases illustrated in Fig. 3 and Supplementary Fig. 13.

The mass budgets for distinct cross-sections of the simplified conceptual model depicted in Supplementary Fig. 14 were computed based on average results from all Monte Carlo model runs (1,450) at stage 1. According to these findings, the mass budget in cross-section A–A is MB_A–A = 57.136 µmol cm⁻² yr⁻¹, in B–B is MB_B–B = 57.197 µmol cm⁻² yr⁻¹, and in C–C is MB_C–C = 57.189 µmol cm⁻² yr⁻¹, indicating an overall mass balance error of ~0.1%, which is less than the permissible mass balance error of 1% considered in our general modelling framework. It is important to mention that, despite our thorough model validation process, in this study, we utilize the model solely for process significance analysis and gaining insights into the underlying mechanisms responsible for OC preservation, rather than for making global predictions, which will be addressed in future research. Finally, based on acceptable uncertainties (95% confidence interval) associated with the total of all >1,000 RTM runs, represented as the shaded area surrounding the curves in Fig. 4 and Supplementary Fig. 16, and the uncertainties of process significance analysis obtained from the ANN phase displayed in Fig. 3, our overarching approach of random variation of input parameters is justified.

It is worth noting that limitations regarding the capacity of sorption sites, for instance, monolayer sorption^16,40, conventionally do not apply to kinetic sorption because the kinetic sorption model in this study predominantly signifies the processes that internalize DOC into the mineral matrix, such as occlusion, co-precipitation, and aggregation. Thus, limited-capacity sorption, as often considered in the literature through the monolayer surface adsorption hypothesis, does not apply to our MOC production. Additionally, kinetic sorption occurs at a slower rate than equilibrium adsorption, which is recognized as instantaneous. Therefore, kinetic sorption, which is restricted by pore water concentrations that are also regulated by hydrolysis, degradation, and other factors, is less likely to encounter a second limitation due to the capacity of sorption sites when compared to instantaneous equilibrium adsorption for which various types of isotherms, including linear, Langmuir, and Freundlich, have been established^66,67. Introducing an additional parameter to enforce a portion of the OC to be absorbed by both kinetic and equilibrium sorption sites would result in a greater number of unknown parameters and is considered unwarranted in this instance.