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

The marine phosphorus cycle plays a critical role in controlling the extent of global primary productivity and thus atmospheric pO2 on geologic time scales. However, previous attempts to model carbon-phosphorus-oxygen feedbacks have neglected key parameters that could shape the global P cycle. Here we present new diagenetic models to fully parameterize marine P burial. We have also coupled this diagenetic framework to a global carbon cycle model. We find that seawater calcium concentration, by strongly influencing carbonate fluorapatite (CFA) formation, is a key factor controlling global phosphorus cycling, and therefore plays a critical role in shaping the global oxygen cycle. A compilation of Cenozoic deep-sea sedimentary phosphorus speciation data provides empirical support for the idea that CFA formation is strongly influenced by marine Ca concentrations. Therefore, we propose a previously overlooked coupling between Phanerozoic tectonic cycles, the major-element composition of seawater, the marine phosphorus cycle, and atmospheric pO2.

Introduction

Although there are pan class="Chemical">fen>w constraints on the how the size of the biosphere has changed through Earth’s history, n>an class="Chemical">phosphorus (P) is commonly considered to be the major limiting nutrient for primary production on geologic time scales. It is also generally accepted that pan class="Chemical">P has a critical role in regulating the extent of organic carbon (C) burial and thus atmospheric oxygen levels (pO2)[1-7]. However, our understanding of the factors that govern the retention or regeneration of P during burial and subsequent chemical alteration of marine sediments (diagenesis) arguably remains limited. Traditionally, bottom-water oxygen concentrations have been assumed to exert the largest influence on P burial[4]. More recently, it has been suggested that the rise of bioturbation in the early Paleozoic significantly increased marine P burial[6,8]. However, although these factors may influence the burial of organic and iron-bound P phases, the dominant controls on the formation of carbonate fluorapatite (CFA)—the largest modern authigenic P burial flux[9]—are relatively poorly understood. As CFA accounts for >50% of marine P burial[9], this translates into considerable uncertainty regarding the chief factors responsible for global P cycling and, by extension, the processes regulating global organic carbon burial and atmospheric pO2 on geologic time scales.

The factors regulating pan class="Chemical">CFAn> formation and burial in marine sediments and the relative importance of envn>an class="Chemical">ironmental and biological boundary conditions in influencing this process have been underexplored. For example, the kinetics of pan class="Chemical">CFA formation are typically simplified for use in diagenetic models; commonly, CFA formation is parameterized as a linear relationship with dissolved phosphate (e.g., refs. [5,8,10]), and the roles played by other key CFA components such as Ca2+, CO32−, and F− are largely overlooked. In particular, calcium (Ca) is the most abundant ionic species in CFA (Ca9.54Na0.33Mg0.13(PO4)4.8(CO3)1.2F2.48, ref. [11]) and seawater dissolved Ca is the single most important Ca source for CFA formation in the upper marine sediment pile. Intuitively, the large swings in seawater Ca concentration during the Phanerozoic[12-14] could have significantly influenced the saturation state and thus precipitation rate of CFA. Here, using models and empirical data, we make a case that seawater Ca concentrations are a key factor shaping marine P burial. Our model results also suggest a significant role for seawater Ca on the evolution of Earth’s redox budget.

Results and discussion

Model evidence for a Ca control on marine P burial

To explore the idea that Ca may exercise a major control on n class="Chemical">pan class="Chemical">P cycling, we first built a simpn>le diagenetic model that includes only organic matter, organic pan class="Chemical">P, pan class="Chemical">Ca2+, CO32−, PO43−, F−, and CFA (see Methods). In contrast to the simplified rate law used in previous studies[5,8,10], this model includes a rate law that describes the formation rate of CFA as a function of its saturation state, which is consistent with experimental observations[15,16]. This exercise indicates that an increase in seawater Ca concentrations will drive a significant increase in CFA formation (Fig. 1 and Supplementary Figure 1), decreasing phosphate diffusion (i.e., P recycling) from the sediment pile back to seawater. In other words, this exercise confirms that Ca promotes CFA formation, which leads to greater P retention in marine sediments (Fig. 1 and Supplementary Figure 1).

Fig. 1

Diagenetic model of seawater Ca influence on CFA and P burial.

PCFA represents burial concentration of CFA-associated phosphorus in marine sediments. a–b The effects of seawater dissolved calcium concentration ([Ca]SW) upon PCFA and P burial efficiency. Labeled fluxes indicate rates of organic matter loading at the sediment–seawater interface.

pan class="Chemical">PCFAn> represents burial concentration of n>an class="Chemical">CFA-associated phosphorus in marine sediments. a–b The effects of seapan class="Chemical">water dissolved calcium concentration ([Ca]SW) upon PCFA and P burial efficiency. Labeled fluxes indicate rates of organic matter loading at the sediment–seawater interface.

To further explore the role n class="Chemical">played by seawater Ca concentration and determine the relative influence of a range of envn>an class="Chemical">ironmental and biotic factors on P burial in marine systems, we also developed an extended, multicomponent reaction-transport marine sediment diagenesis model, which includes over 30 diagenetic reactions, a saturation state-dependent rate law for CFA formation, and complete, detailed parameterization of porewater pH, carbonate chemistry, and adsorption (see Methods and Supplementary Tables 1–5). We have also included the iron phosphate mineral vivianite (which is a significant P sink in some restricted marine settings today[8]). We calibrated our CFA rate law in our fully parameterized model using sedimentary porewater and solid-phase data collected from a well-characterized modern shallow marine site (ref. [3], Friends of Anoxic Mud site, FOAM) and the deep-sea ODP site 1226 (Supplementary Figures 2–3, see Supplementary Methods for further description). Using this extended model, we then explored the relative influence of a range of environmental and biotic factors on P burial.

Our model results show that the intensity of bioturbation, the magnitude of the pan class="Disease">organic matter fluxn> to the sediment surface, and the concentration of dissolved pan class="Chemical">oxygen in bottom pan class="Chemical">waters all have non-negligible and nonlinear effects on CFA formation and total P burial (Fig. 2 and Supplementary Figures 4–6). Our model, like other recently developed models[8], also indicates that increases in bioturbation intensity and depth can be associated with increases in CFA and total P burial, although enhancement of P burial by bioturbation is considerably more muted or even reversed at higher bioturbation intensities and depths (Fig. 2 and Supplementary Figures 4–6). Similarly, we find that CFA burial does not have a simple, linear relationship with organic matter flux or bottom-water oxygen levels (Fig. 2 and Supplementary Figures 4–6). Changes in bioturbation, organic matter flux, and bottom-water oxygen levels shape multiple factors controlling the saturation state of CFA, which can lead to nonlinear changes in P burial efficiency. For example, higher organic matter fluxes can promote CFA burial by providing more dissolved inorganic phosphate[17]. However, increased organic matter fluxes can also dampen CFA burial, when there is a sharp drop in the pH of pore waters in the upper portion of the sediment pile. Although our model provides new insights into P burial (Fig. 2), our results, overall, corroborate those of other modeling studies that have documented similar nonlinear behavior of CFA burial in response to variability in environmental factors[17].

Fig. 2

Nonlinear effects on CFA, total P burial, and Corg/Preac.

PCFA (circles) represents burial concentration of CFA-associated phosphorus in marine sediments. BE (diamonds) represents the burial efficiency of reactive phosphorus in marine sediments. Corg/Preac represents the burial ratio between organic carbon and reactive P. a–d The effects upon PCFA, BE, and Corg/Preac of bottom-water oxygen concentration a, seawater dissolved calcium concentration b, flux of organic carbon to the sediment–seawater interface (Joc) c, and bioturbation d. For bioturbation, the four parameters included in the model (DB0, a0, xbt, and xbi) were increased linearly from zero to modern average values, whereas a linear increase of 16% was applied to porosity. Bioturbation is parameterized here as a coupled biodiffusion and bioirrigation term. See Supplementary Tables 1–5 for full list of model parameters.

pan class="Chemical">PCFAn> (circles) represents burial concentration of n>an class="Chemical">CFA-associated phosphorus in marine sediments. BE (diamonds) represents the burial efficiency of reactive phosphorus in marine sediments. Corg/pan class="Chemical">Preac represents the burial ratio between organic carbon and reactive P. a–d The effects upon PCFA, BE, and Corg/Preac of bottom-water oxygen concentration a, seawater dissolved calcium concentration b, flux of organic carbon to the sediment–seawater interface (Joc) c, and bioturbation d. For bioturbation, the four parameters included in the model (DB0, a0, xbt, and xbi) were increased linearly from zero to modern average values, whereas a linear increase of 16% was applied to porosity. Bioturbation is parameterized here as a coupled biodiffusion and bioirrigation term. See Supplementary Tables 1–5 for full list of model parameters.

Despite significant increases in the comn class="Chemical">plexity of our extended diagenetic model, we find, in agreement with the results of our basic diagenetic model, that marine Ca concentrations markedly alter the efficiency of P burial. n>an class="Chemical">Calcium concentrations, by strongly influencing CFA saturation state, have a strong impact on overall P burial relative to other environmental factors (Figs. 2–3 and Supplementary Figures 4–7). Unlike other environmental forcings—such as bioturbation, organic matter flux, and bottom-water oxygen levels—higher seawater Ca concentrations always induced more CFA burial (Fig. 2 and Supplementary Figures 4–6), owing to its direct effect on the saturation state of CFA. Increases in CFA burial at high seawater Ca concentrations mediate decreases in the burial ratio between organic C and reactive P (Corg/Preac, Fig. 2, and Supplementary Figures 4–5). Given significant secular variability in seawater dissolved Ca concentrations over the last 80 million years as well as, more broadly, throughout the Phanerozoic (between ≤10 and 30 mm over the last ~500 million years)[12-14], shifts in dissolved calcium concentrations are likely to have been a key environmental factor controlling P burial efficiency over this interval.

Fig. 3

Influence of environmental factors on CFA burial.

a The effect of environmental forcings on CFA burial. The y axis denotes the concentration of CFA-associated P in sediments. b the effect of environmental forcings on the burial ratio between organic carbon and reactive P. The horizontal line represents the value for the shallow marine reference model run (see Supplementary Methods). All parameters were held constant, apart from the individual parameter varied for each sensitivity analysis. Bioturb is a coupled biodiffusion and bioirrigation term and Joc is the flux of organic matter to the sediment-water interface.

a The efpan class="Chemical">fen>ct of environmental forcings on pan class="Chemical">CFA burial. The y axis denotes the concentration of CFA-associated P in sediments. b the effect of environmental forcings on the burial ratio between organic carbon and reactive P. The horizontal line represents the value for the shallow marine reference model run (see Supplementary Methods). All parameters were held constant, apart from the individual parameter varied for each sensitivity analysis. Bioturb is a coupled biodiffusion and bioirrigation term and Joc is the flux of organic matter to the sediment-water interface.

Bottom-pan class="Chemical">watern> dissolved inpan class="Chemical">organic carbon (DIC), as well as dissolved pan class="Chemical">sulfate, magnesium, and phosphate concentrations appear to have relatively little impact on CFA burial (Fig. 3). There are likely two chief reasons why bottom-water DIC and phosphate concentrations do not significantly mediate CFA burial. The first is that the stoichiometric abundance of CO32−and PO43− in CFA is low, relative to that of Ca2+ (ref. [11]). Moreover, porewater DIC and dissolved inorganic phosphorus (DIP) are mainly generated by the decomposition of organic matter and/or the reduction of iron oxides[5,8], rather than diffusion from seawater. The Mg concentration of seawater also does not appear to exercise a strong control on CFA burial, likely because Mg is a minor component of CFA. Bottom-water pH also appears to not strongly influence CFA burial (Fig. 3). This is likely because sedimentary and porewater reactions, rather than seawater, are the most important source of the protons within the sediment pile. In sum, we find that bottom-water oxygen, the extent of organic matter loading, and bioturbation play a sizeable role in controlling marine P burial, consistent with previous studies[6,8]. However, our results also show that seawater Ca concentrations have a major role in controlling P burial.

Empirical evidence for a Ca control on marine P burial

A Ca control on marine pan class="Chemical">Pn> burial is also supported by n>an class="Chemical">P speciation data compiled from Cenozoic deep-sea sediment cores across the Pacific and Atlantic Ocean basins (Fig. 4; see Supplementary Methods; further characteristics of the compiled sites are also shown in Supplementary Table 6). After filtering data from the uppermost sediment pile that are still undergoing active diagenesis (see Supplementary Methods), these deep-sea data suggest relatively constant pan class="Chemical">CFA burial between 80 Ma and 40 Ma and a gradual decrease in CFA burial starting ~40 Ma (Fig. 4). This trend can be seen in both the Pacific and Atlantic sites (Fig. 4). Meanwhile, empirical Mg/Ca data from fluid inclusions, biogenic carbonates, and calcite veins in oceanic crust suggest a gradual decrease in seawater Ca concentrations also began ~40 million years ago (refs. [12-14], Fig. 4). Thus, secular trends in CFA burial correspond well with the timing of changes in seawater Ca concentration. Coupled with the prediction from our model that Ca concentrations will directly affect CFA saturation state and thus CFA burial, this covariation suggests that seawater Ca concentration is, through its influence of CFA formation, a major long-term forcing on global P burial.

Fig. 4

Coupling of CFA burial and calcium over 80 million years.

a Changes in the burial ratio of CFA-associated P to total reactive P (a sum of organic P, CFA, and iron-bound P) in deep-sea sediments. Points represent mean values for each 2 myr bin, with the error bars represent one standard deviation (1σ). The black curve represents deep-sea model results from the coupled C-P-O model. See Supplementary Information for discussion of the data compilation. b Seawater Mg/Ca ratios recorded in CaCO3 veins of oceanic crust (circles)[14], fluid inclusions (triangles)[12] and echinoderm ossicles (squares)[13]. The curve denotes one of the estimates for seawater Ca concentrations through this interval[12].

a Changes in the burial ratio of pan class="Chemical">CFAn>-associated pan class="Chemical">P to total reactive pan class="Chemical">P (a sum of organic P, CFA, and iron-bound P) in deep-sea sediments. Points represent mean values for each 2 myr bin, with the error bars represent one standard deviation (1σ). The black curve represents deep-sea model results from the coupled C-P-O model. See Supplementary Information for discussion of the data compilation. b Seawater Mg/Ca ratios recorded in CaCO3 veins of oceanic crust (circles)[14], fluid inclusions (triangles)[12] and echinoderm ossicles (squares)[13]. The curve denotes one of the estimates for seawater Ca concentrations through this interval[12].

The influence of calcium on atmospheric oxygenation

In order to explore the influence of variation in marine Ca concentrations on global n class="Chemical">pan class="Chemical">organic carbon burial and atmospn>heric pan class="Chemical">pO2, we coupled our extended diagenetic model to a simple global pan class="Chemical">carbon cycle mass balance model. Specifically, we used the outputs of our extended diagenetic model to parameterize how P burial efficiencies vary as a function of bottom-water oxygen, extent of organic matter loading, marine Ca concentrations, and bioturbation in a global carbon cycle mass balance model (modified from ref. [4], Supplementary Tables 7–13). The primary goal of our global carbon cycle modeling approach was to test to what extent changes in seawater dissolved Ca concentration (as indicated by geologic archives) could have influenced contemporaneous atmospheric pO2 (Fig. 5 and Supplementary Figures 9–14). Full details of the carbon cycle mass balance model (and the full suite of parameters employed) are provided in the Methods and the Supplementary Information. We find that increases in seawater dissolved Ca concentration would, by inducing increased CFA precipitation and decreasing P recycling, lead to decreased seawater P concentration (Fig. 5). More importantly, an increase in CFA precipitation decreases the burial ratio between organic carbon and reactive phosphorus (Fig. 2)—this, at a constant P flux to the oceans, in turn will lead to a drop in atmospheric oxygen levels.

Fig. 5

Atmospheric oxygen levels respond to changes in calcium or bioturbation.

a Model input of seawater dissolved Ca concentration or bioturbation (0 represents no bioturbation, whereas 1 represents the reach of modern bioturbation). b–f Model outputs of variations in P burial associated with CFA formation, total marine P (a sum of dissolved inorganic phosphate, dissolved organic phosphorus and soluble particulate inorganic phosphorus), net primary productivity (NPP), organic P burial flux and atmospheric oxygen levels. pO2 shown in the figure is the actual percent by volume.

a Model input of sean class="Chemical">pan class="Chemical">water dissolved Ca concentration or bioturbation (0 repclass="Chemical">n>resents no bioturbation, whereas 1 represents the reach of modern bioturbation). b–f Model oun>an class="Chemical">tputs of variations in P burial associated with pan class="Chemical">CFA formation, total marine P (a sum of dissolved inorganic phosphate, dissolved organic phosphorus and soluble particulate inorganic phosphorus), net primary productivity (NPP), organic P burial flux and atmospheric oxygen levels. pO2 shown in the figure is the actual percent by volume.

Our model results suggest that Ca could exercise a strong external control on atmospheric n class="Chemical">pan class="Chemical">pO2 (Fig. 5). For exampn>le, an increase in seapan class="Chemical">water Ca concentration from 10 mm to 20 mm could lead to a >50% decrease in atmospheric pan class="Chemical">pO2, from modern levels 21% to ~10% (Fig. 5). There is a small rebound in the total marine P reservoir after a decrease induced by a shift in [Ca2+] (Fig. 5c), which is owing to the delayed decrease in atmospheric oxygen and thus iron-bound P burial.

Given the large swings in pan class="Chemical">Pn>hanerozoic seapan class="Chemical">water Ca concentration recorded by geologic archives[12-14,18,19], our results suggest that seapan class="Chemical">water Ca concentration could be one of the key factors shaping atmospheric pO2 evolution. For example, proxy records indicate low seawater dissolved Ca concentrations during the Carboniferous–Permian[12,20], which would have inhibited CFA formation and driven the development of the high atmospheric pO2 that, on the basis of both geochemical and paleontological archives[21-25], has been previously suggested to be characteristic of this time (interval 3 of Supplementary Figure 12). Fluid-inclusion data and carbonate mineralogical records are also commonly interpreted to indicate an increase in seawater dissolved Ca concentrations occurred during the early Cambrian[12,26-28], coincident with what has been previously suggested to have been an interval of ocean deoxygenation (e.g., refs. [6,23,29]).

Our modeling results indicate that seapan class="Chemical">watern> Ca concentration is likely a major driver of total pan class="Chemical">P burial, such that, over Earth’s history, high seapan class="Chemical">water dissolved Ca concentrations should correspond to significant increases in P burial efficiency and relatively low atmospheric pO2. A compilation of deep-sea P speciation data from the past 80 million years provides empirical support for the effect of seawater dissolved Ca concentrations on P burial and bolster the case that this Ca-driven feedback strongly impacts the global ocean-atmosphere system. These results provide a new view of the processes linking marine elemental cycles, tectonics and atmospheric pO2, and suggest that the major ion composition of seawater has been an important driver of biospheric change and atmospheric evolution throughout Earth’s history.

Methods

Basic 1D diagenetic model

The basic reactive-transport diagenetic model we develon class="Chemical">ped includes six components, including two solid-phase species: organic matter and CFA; and four solute compn>onents: [n>an class="Chemical">Ca2+], [DIC], [DIP], and [F−]. The mass balance functions for this basic model are[30,31]:where C is the concentration of solute, C is the concentration of solid, ϕ is porosity, D is molecular diffusion, ʋ is the solute advection rate, ω is solid advection rate, and R and R are, respectively, the solid and solute reaction rates. The effect of molecular diffusion is calibrated to tortuosity through , where D 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. [30]. The depth dependence of porosity was represented by ϕ(x) = ϕ∞ + (ϕ0 − ϕ∞) exp(−x/λ), where λ is the porosity attenuation length, and ϕ0 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, [Mg2+], 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[32] was used to describe the decomposition of organic matter in this basic model iteration.

pan class="Chemical">CFAn> is the largest pan class="Chemical">P sink in the global marine pan class="Chemical">P cycle[3]. In this study, we use the stoichiometry of a carbonate-rich francolite[11,33]—Ca9.54Na0.33Mg0.13(PO4)4.8(CO3)1.2F2.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 (Ω) is higher than 1, the precipitation rate of CFA can be expressed as:where k is a kinetic constant for CFA, obtained by fitting to FOAM porewater and sedimentary CFA profiles. Ω is expressed as:where rCa, rNa, rMg, rPO, rCO, and rF represent the activity coefficients of each ion. K is the equilibrium constant of CFA, which was found to be a function of carbonate activity[11]: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 (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 CO32− in CFA correlates with [CO32−] of porewater[11]. However, the relationship between the solubility of CFA, the stoichiometric abundance of CO32−in CFA and [CO32−] of porewater are not fully understood[11]. These uncertainties may influence the relationship between [CO32−] 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 and the saturation of CFA (Ω). Following the approach of previous studies[8,10], the dissolution of CFA under Ω < 1 is not included, as CFA is highly insoluble in marine sediments.

Extended 1D multicomponent diagenetic model

Building from previous efforts[8,34-36], we also built an extended 1D multicomn class="Chemical">ponent model incorporating the biogeochemical cycles of C, N, P, S, n>an class="Chemical">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 DB is the intensity of biodiffusion, a is the coefficient of bioirrigation and C 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 and , respectively. DB0 and a0 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 Fe2+ and Mn2+ (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

pan class="Chemical">Pn>arameterizations of the diagenetic pan class="Chemical">P cycle explored in this study share many pan class="Chemical">features with previously published model exercises (e.g., refs. [8,10,36]). Three solid P species—organic phosphorus (orgP), iron-bound phosphorus (PFe), and CFA—as well as dissolved inorganic phosphate (ƩPO43−) 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 simn class="Chemical">plified as CNP, where and are the molar ratios of organic n>an class="Chemical">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[40], 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[41] developed from a reactive continuum-tyn class="Chemical">pe model[32]. 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[41]. 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[41]. The rate constant for each Gi iswhile the fraction of each Gi iswhere a is the average lipan class="Chemical">fetime of more reactive orgC and v is the shapn>e of orgC distribution[32].

The C/pan class="Chemical">Pn> ratio of particulate organic matter typically increases with depn>th 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 pan class="Chemical">organic phosphorus and dissolved phosphate (Supplementary Figure 2). The α pool includes the first 2 G fractions (G1 and G2), whereas the β pool includes the remaining 10 G fractions (G3–G12). The C/P ratios of the two pools are shown in Supplementary Table 5.

Five pan class="Chemical">ironn> phases are included in the modeled flux to the sediment–sean>an class="Chemical">water interface: (1) highly reactive iron hydroxides (pan class="Chemical">Fe(OH)3α); (2) less reactive iron hydroxides (Fe(OH)3β); (3) unreactive iron hydroxides (Fe(OH)3γ); (4) magnetite (Fe3O4); 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[43]. Iron may also occur as other silicate-bound phases, but we use biotite as a representation of silicate iron (see ref. [44]). At FOAM, biotite dissolution appears to be closely linked to porewater pH. As for previous models[8], 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 difpan class="Chemical">ference of this study from previous studies[5,8,10] is that we have parameterized the formation rate of n>an class="Chemical">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 pan class="Chemical">iron phosphaten> mineral vivianite in the extended model. Although it has not been extensively documented in marine sediments, it is commonly found in restricted settings[45]. Following ref. [45], we have used the Michaelis–Menten kinetipan class="Chemical">cs for dissolved pan class="Chemical">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.[46] to simulate pH variation during diagenetic processes. Instead of total alkalinity, we used total proton balance (TP) as an equilibrium-invariant and impn>licit differential variable (whose differentials are not treated as a main function and are only used to calculate the 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 H3PO4 and HS− in the formulation of TP. Following Hoffman et al.[46], 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 is the total concentration of the dissociation systems , , , and . A derivation of and the formulation of and can be found in the Supplementary Text. In the extended model, the proton mass balance is calculated using Eq. 11, and the 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 n class="Chemical">protons and ions such as Fe2+ and n>an class="Chemical">Mn2+, 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 Fe2+, KFe can be defined as the amount of iron in the adsorbed phase relative to Fe2+ in the solute, thuswhere AFe is the concentration of adsorbed solid-phase iron, and . In the model, Fe2+ is treated as a differential variable, whereas AFe is treated as an algebraic variable, with the assumption of instantaneous equilibrium between them. As Fe2+ is influenced by both the reaction and transport of Fe2+ and AFe, the derivative of Fe2+ iswhere and are the sum of the reaction and transport terms (excluding adsoption) of Fe2+ and AFe, respectively, which are calculated from Eqs. 6 and 7. The transfer rate () of protons during the adsorption process isDerivations of Eqs. 13 and 14 can be found in the Supplementary Text. The adsorption of NH4+ 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 Mn2+ adsorption, accurately reproduce the FOAM dissolved Mn2+ profile (Supplementary Figure 2), the adsorption of Mn2+ was also not included.

Coupled carbon–phosphorus–oxygen cycle model

Our coupled n class="Chemical">pan class="Chemical">carbon–pn>an class="Chemical">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 pO2 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 apn class="Chemical">proach 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 (interpn>olated between time slices of 1 million years). Rates of n>an class="Chemical">organic carbon burial (which are, in turn, influenced by organic carbon remineralization) and organic carbon weathering are, in our model, influenced by atmospheric O2 level, following the parameterization of the COPSE model[48]. 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.