The role of calcium in regulating marine phosphorus burial and atmospheric oxygenation

Basic 1D diagenetic model

The basic reactive-transport diagenetic model we developed includes six components, including two solid-phase species: organic matter and CFA; and four solute components: [Ca²⁺], [DIC], [DIP], and [F⁻]. The mass balance functions for this basic model are^30,31:

where C_l is the concentration of solute, C_s is the concentration of solid, ϕ is porosity, D is molecular diffusion, ʋ is the solute advection rate, ω is solid advection rate, and R_s and R_l are, respectively, the solid and solute reaction rates. The effect of molecular diffusion is calibrated to tortuosity through ({mathrm{D}} = frac{{{mathrm{D}}_m}}{{1 – 2ln left( phi right)}}), where D_m is the intensity of molecular diffusion. The effect of compaction on the solute and solid advection rates was also considered using the method described in ref. ³⁰. The depth dependence of porosity was represented by ϕ(x) = ϕ_∞ + (ϕ₀ − ϕ_∞) exp(−x/λ), where λ is the porosity attenuation length, and ϕ₀ and ϕ_∞ represent porosity at the sediment–seawater interface and at depth, respectively. The boundary conditions for these species and the above model parameters were also used for the extended, multicomponent version of the diagenetic model (described in more detail below; Supplementary Tables 1 and 5). The porewater profiles of pH, [Mg²⁺], and [Na⁺] were, for this basic iteration of the diagenetic model, fixed at seawater levels. The Redfield ratio of 106:1 was used as the C/P ratio of the organic matter. There are only two reactions in this basic model, the decomposition of organic matter and the formation of CFA. A reactive continuum-type model³² was used to describe the decomposition of organic matter in this basic model iteration.

CFA is the largest P sink in the global marine P cycle³. In this study, we use the stoichiometry of a carbonate-rich francolite^11,33—Ca_9.54Na_0.33Mg_0.13(PO₄)_4.8(CO₃)_1.2F_2.48—although the stoichiometry of CFA can be variable. Previous studies have simplified the precipitation rate of CFA by presuming a linear relationship with excess dissolved phosphate (relative to an assumed equilibrium phosphate concentration, refs. ^5,8,10). However, experimental studies have shown that the rates of the crystal growth of calcium hydroxylapatite and calcium fluorapatite are functions of saturation state^15,16. Thus, in this study, we more completely parameterize the kinetics of CFA formation by assuming that the formation rate of CFA is characterized by a first-order relationship with its saturation state. When the saturation state of CFA (Ω_CFA) is higher than 1, the precipitation rate of CFA can be expressed as:

where k_CFA is a kinetic constant for CFA, obtained by fitting to FOAM porewater and sedimentary CFA profiles. Ω_CFA is expressed as:

where rCa, rNa, rMg, rPO₄, rCO₃, and rF represent the activity coefficients of each ion. K_spCFA is the equilibrium constant of CFA, which was found to be a function of carbonate activity¹¹:

In this parametrization, carbonate ion activity influences the saturation state of CFA through its effect on both IAP and Ksp. We have also used the same stoichiometry of CFA and a fixed K_spCFA (10^−99.7) without a relationship with carbonate activity in a series of runs (Supplementary Figures 5, 6, and 10) and found that this change only has a subtle influence on the model results (Figs. 2 and 5, and Supplementary Figures 5 and 10). It is also possible that stoichiometric abundance of CO₃²⁻ in CFA correlates with [CO₃²⁻] of porewater¹¹. However, the relationship between the solubility of CFA, the stoichiometric abundance of CO₃²⁻in CFA and [CO₃²⁻] of porewater are not fully understood¹¹. These uncertainties may influence the relationship between [CO₃²⁻] of porewater and CFA formation, which is currently a negative correlation (Fig. 3). Further, these will not influence our main conclusions as they do not significantly influence the relationship between ([{mathrm{Ca}}^{2 + }]) and the saturation of CFA (Ω_CFA). Following the approach of previous studies^8,10, the dissolution of CFA under Ω_CFA < 1 is not included, as CFA is highly insoluble in marine sediments.

Extended 1D multicomponent diagenetic model

Building from previous efforts^8,34,35,36, we also built an extended 1D multicomponent model incorporating the biogeochemical cycles of C, N, P, S, Fe, and Mn to more fully simulate the diagenetic P cycle in marine sediments (SEDCHEM). The model includes 15 solutes and 21 solids (Supplementary Table 1). A biodiffusion term was used to describe the mixing of sediment particles. A non-local function^37,38 was applied to describe the influence of bioirrigation on the exchange of solutes near the sediment–seawater interface. Combining the molecular diffusion, advection, bioturbation, and reaction terms, the mass balance of solutes and solids can be generalized to the following functions^30,31:

where D_B is the intensity of biodiffusion, a is the coefficient of bioirrigation and C_l0 is the solute concentration in open burrows, which is assumed to be equivalent to the solute concentration of the overlying water column. The attenuations in the intensities of biodiffusion and bioirrigation at depth are described by ({mathrm{D}}_{mathrm{B}}left( {mathrm{x}} right) = {mathrm{D}}_{{mathrm{B}}0} cdot exp left( { – left( {frac{x}{{xbt}}} right)^2} right)) and ({mathrm{a}}left( {mathrm{x}} right) = {mathrm{r}} cdot {mathrm{a}}_0 cdot exp left( { – frac{x}{{xbi}}} right)), respectively. D_B0 and a₀ are the biodiffusion and bioirrigation intensities at the sediment–seawater interface, xbt and xbi are the respective attenuation coefficients and r is a correction for the irrigation of Fe²⁺ and Mn²⁺ (refs. ^35,39). Model components and boundary conditions are shown in Supplementary Table 1. Model reactions, reaction rate laws and parameters can be found in Supplementary Tables 2–5. The following three sections delineate treatment of the diagenetic P cycle, pH, and adsorption in the extended 1D multicomponent diagenetic model. See Supplementary Text for model solution and model applications.

Diagenetic phosphorus cycle

Parameterizations of the diagenetic P cycle explored in this study share many features with previously published model exercises (e.g., refs. ^8,10,36). Three solid P species—organic phosphorus (orgP), iron-bound phosphorus (P_Fe), and CFA—as well as dissolved inorganic phosphate (ƩPO₄³⁻) are included in the extended version of our diagenetic model. In contrast to previous studies, however, our model also includes a more intrinsic reaction rate law for CFA formation, which involves all the major components of CFA and parameterizes the effect of pH on P speciation and burial. To simulate pH, our model also includes proper descriptions of adsorption, the diagenetic Fe cycle, and authigenic carbonate precipitation and dissolution.

Decomposition of organic matter drives diagenetic reactions. The chemical formula of organic matter can be simplified as CN_rNP_rP, where _rN and _rP are the molar ratios of organic nitrogen and organic phosphorus relative to organic carbon. The major pathways for the decomposition of organic matter, including aerobic respiration, nitrate reduction, Mn reduction, Fe reduction, sulfate reduction, and methanogenesis are included in the model. As the sequence (and thus spatial distribution) of these pathways is dictated by their respective energy yields⁴⁰, a Monod scheme was used to describe the operation of these reactions^8,10,36.

For the mineralization of orgC in the extended model, we used a multi-G model⁴¹ developed from a reactive continuum-type model³². The advantage of this model is that it not only can represent the reactive continuum of organic matter, but also is suitable to apply to the bioturbated zone of the sediment pile⁴¹. We divided organic matter into 12 G in this study. Although each G has its own first-order rate constant, their fractions in total organic matter are determined by a gamma distribution⁴¹. The rate constant for each G_i is

while the fraction of each G_i is

where a is the average lifetime of more reactive orgC and v is the shape of orgC distribution³².

The C/P ratio of particulate organic matter typically increases with depth in the sediment pile. Following previous methods^36,42, we assume that this variation is generated by different C/P ratios in different organic components, with less reactive (more refractory) organic matter having higher C/P ratios. Thus, we further divided the 12 G fractions into two pools (α and β) with different C/P ratios. With this procedure, it is possible to reproduce empirical sedimentary profiles of organic phosphorus and dissolved phosphate (Supplementary Figure 2). The α pool includes the first 2 G fractions (G₁ and G₂), whereas the β pool includes the remaining 10 G fractions (G₃–G₁₂). The C/P ratios of the two pools are shown in Supplementary Table 5.

Five iron phases are included in the modeled flux to the sediment–seawater interface: (1) highly reactive iron hydroxides (Fe(OH)₃^α); (2) less reactive iron hydroxides (Fe(OH)₃^β); (3) unreactive iron hydroxides (Fe(OH)₃^γ); (4) magnetite (Fe₃O₄); and (5) biotite (Biot). Demarcation of iron hydroxides by reactivity is similar to the treatment employed by previous models^8,36. The rate law for magnetite dissolution is reasonably well established⁴³. Iron may also occur as other silicate-bound phases, but we use biotite as a representation of silicate iron (see ref. ⁴⁴). At FOAM, biotite dissolution appears to be closely linked to porewater pH. As for previous models⁸, iron-bound phosphorus is assumed to be associated with iron hydroxides, and the P/Fe molar ratio is described using a constant γ and a variable θ (Supplementary Tables 1–5). Thus, the precipitation and dissolution of iron hydroxides are accompanied by, respectively, the scavenging and release of dissolved phosphate (Supplementary Table 2).

A major difference of this study from previous studies^5,8,10 is that we have parameterized the formation rate of CFA using its saturation state (see above), which is consistent with experimental results^15,16. A full description of this method can be found above, as well as in the Supplementary Text and Supplementary Table 3.

We have also included the iron phosphate mineral vivianite in the extended model. Although it has not been extensively documented in marine sediments, it is commonly found in restricted settings⁴⁵. Following ref. ⁴⁵, we have used the Michaelis–Menten kinetics for dissolved phosphate and iron in marine porewater to describe the formation of vivianite. Details of the formulation and parameters for vivianite can be found in Supplementary Tables 3 and 5.

Diagenetic pH simulation

We used a method similar to that of Hoffman et al.⁴⁶ to simulate pH variation during diagenetic processes. Instead of total alkalinity, we used total proton balance (TP) as an equilibrium-invariant and implicit differential variable (whose differentials (frac{{{mathrm{dy}}}}{{dt}}) are not treated as a main function and are only used to calculate the (frac{{{mathrm{dy}}}}{{dt}}) of the differential variable) to model pH. Although the results of the two methods are the same, the new method is more straightforward and less-demanding computationally. We define the total proton balance as the sum of those ions that can undergo proton exchange at the modeled pH range (~6–9), which can be written as:

As the first dissociation constant of phosphoric acid is high and the second dissociation constant of hydrogen sulfide is very low, we do not include H₃PO₄ and HS⁻ in the formulation of TP. Following Hoffman et al.⁴⁶, we treat TP as an implicit differential variable in the model, which removes the need to solve the algebraic equation numerically. The derivative of the proton concentration is:

where ([ {S_j}]) is the total concentration of the dissociation systems (big(mathop {sum }nolimits CO_2), (mathop {sum }nolimits NH_3), (mathop {sum }nolimits {mathrm{H}}_3{mathrm{PO}}_4), and (mathop {sum }nolimits H_2S big)). A derivation of (frac{{dH}}{{dt}}) and the formulation of (frac{{partial TP}}{{partial left[ {S_j} right]}}) and (frac{{partial TP}}{{partial H}}) can be found in the Supplementary Text. In the extended model, the proton mass balance is calculated using Eq. 11, and the (frac{{dTP}}{{dt}}) term is the sum of the mass balance functions of all its species calculated using Eq. 6.

Diagenetic adsorption simulation

Adsorption within marine sediments often involves exchange between protons and ions such as Fe²⁺ and Mn²⁺, which is important in modeling pH variation during diagenetic processes. Adsorption is treated here as a reversible linear equilibrium process. For instance, for the adsorption of Fe²⁺, K_Fe can be defined as the amount of iron in the adsorbed phase relative to Fe²⁺ in the solute, thus

where A_Fe is the concentration of adsorbed solid-phase iron, and (F = (1 – phi )/phi). In the model, Fe²⁺ is treated as a differential variable, whereas A_Fe is treated as an algebraic variable, with the assumption of instantaneous equilibrium between them. As Fe²⁺ is influenced by both the reaction and transport of Fe²⁺ and A_Fe, the derivative of Fe²⁺ is

where (RT_{Fe^{2 + }}) and (RT_{{mathrm{A}}_{{mathrm{Fe}}}})are the sum of the reaction and transport terms (excluding adsoption) of Fe²⁺ and A_Fe, respectively, which are calculated from Eqs. 6 and 7. The transfer rate ((RA_H)) of protons during the adsorption process is

Derivations of Eqs. 13 and 14 can be found in the Supplementary Text. The adsorption of NH₄⁺ was not included in the model as it does not influence pathways of P cycling. However, the C/N ratios of organic matter were tuned to fit the FOAM porewater profile (Supplementary Figure 2). As the model can, without parameterization of Mn²⁺ adsorption, accurately reproduce the FOAM dissolved Mn²⁺ profile (Supplementary Figure 2), the adsorption of Mn²⁺ was also not included.

Coupled carbon–phosphorus–oxygen cycle model

Our coupled carbon–phosphorus–oxygen cycle model integrates our newly developed diagenetic model with a global biogeochemical cycling model modified from the carbon–phosphorus cycle model of Van Cappellen and Ingall^4,47. The purpose of this coupled diagenetic-global biogeochemical cycling modeling exercise is to illustrate the relative importance of various environmental and biotic factors in controlling carbon and phosphorus cycling and atmospheric pO₂ levels. Flux equations, reservoir equations, and parameters are similar to those presented in Van Cappellen and Ingall^4,47. Detailed description of the reservoirs, fluxes, and parameters used in the model (and values for each of these) are listed in Supplementary Tables 9–12.

In contrast to the approach of Van Cappellen and Ingall^4,47, in our coupled model the marine sedimentary burial of calcium-bound phosphate (F58) is calculated from our diagenetic model (interpolated between time slices of 1 million years). Rates of organic carbon burial (which are, in turn, influenced by organic carbon remineralization) and organic carbon weathering are, in our model, influenced by atmospheric O₂ level, following the parameterization of the COPSE model⁴⁸. In particular, we have also, for these runs, parameterized the organic C/P ratio as a function of seawater oxygen level in our diagenetic model. This is achieved by parameterizing two factors rP1 and rP2 using the following functions:

To couple the diagenetic model and the carbon–phosphorus–oxygen cycle model, we built high-resolution look-up tables for CFA burial during the Phanerozoic, using the results of our diagenetic model. We carried out a series of model runs with different bottom-water oxygen concentrations and organic carbon fluxes to the sediment–seawater interface to build look-up tables for the deep sea and shallow oceans, respectively, which were then used to force P burial at each time step. More discussion of the coupled carbon–phosphorus–oxygen cycle model can be found in the Supplementary Text.

Source: Ecology - nature.com