1. Introduction

The global carbon (C) and phosphorus (P) cycles are directly coupled, as both elements are major components of organic matter and critical building blocks for life (Föllmi, 1996). Phosphorus is widely acknowledged to be the ultimate limiting nutrient on geologic time scales (e.g., Reinhard et al., 2017; Tyrrell, 1999). The global phosphorus cycle thus plays an important role in regulating the size of the biosphere, organic carbon cycling, ecology, climate and redox budgets. The co-evolution of climate change, P cycling, and primary productivity in the geological past (especially over late Cenozoic glacial-interglacial cycles and during carbon injections) has therefore been a topic of considerable interest (Boyle & Keigwin, 1987; Broecker, 1982; Filippelli et al., 2007; Ganeshram et al., 2002; Jiang et al., 2022; Komar & Zeebe, 2017; Winguth et al., 2012; Zhao, Zhao, et al., 2023). There has also been extensive recent interest in the effects of warming on P cycling (e.g., Griffiths et al., 2017).

Previous studies have suggested that temperature could be a strong lever on global phosphorus cycling, by means of its effect on the weathering of phosphorus-rich minerals such as apatite over geologic time scales (Guidry & Mackenzie, 2000, 2003). On the other hand, temperature may also have a significant influence on the burial fluxes of different P species (e.g., Ganeshram et al., 2002), which may partially offset the effect of increased P delivery from weathering during warmer times. The partitioning between P species has been proposed to substantially influence organic carbon cycling based on the assumption that oxic and anoxic burial environments are characterized by distinct organic C:P burial ratios (e.g., Lenton et al., 2018; Van Cappellen & Ingall, 1994, 1996). However, the current framework for global P cycling is currently hampered by limited understanding of the effect of temperature change on marine P burial and the formation of authigenic P phases in particular (cf. Van Cappellen & Berner, 1991).

Temperature shapes not only the reaction kinetics but also various boundary conditions for benthic diagenetic processes, which may in turn significantly impact marine P burial. Temperature generally correlates positively with chemical reaction rates, including the decomposition of organic matter (e.g., Aller, 1980; Arndt et al., 2013; Burdige, 2011; Klump & Martens, 1989; Knoblauch & Jorgensen, 1999; Middelburg et al., 1996; Robador et al., 2009; Sawicka et al., 2012; Weston & Joye, 2005; Westrich & Berner, 1988) and diagenetic precipitation of carbonate fluorapatite (CFA; Van Cappellen & Berner, 1991), which is the largest modern marine P burial sink (Ruttenberg & Berner, 1993). As documented extensively by previous work, temperature also has a significant influence on marine pH and dissolved oxygen levels (Breitburg et al., 2018; Caldeira & Wickett, 2003; Doney et al., 2009; Hönisch et al., 2012; Keeling et al., 2010). Temperature strongly regulates the physiological tolerances and thus the activities of benthic animals, including intensities of bioturbation (Gerino et al., 1998; Ouellette et al., 2004; Sun et al., 1991; Teal et al., 2008). It has also been suggested that climate change may have a significant influence on primary productivity and thus rates of organic matter delivery to marine sediments (e.g., Behrenfeld et al., 2006; Boyce et al., 2010). While many environmental parameters such as dissolved oxygen, sedimentary loading of organic matter and bioturbation have been previously demonstrated to have a significant influence on marine P burial (Dale et al., 2016; Föllmi, 1995, 1996; Ingall et al., 1993; Tarhan et al., 2021; Van Cappellen & Ingall, 1994, 1996; Zhao et al., 2020, 2021), the combined effects of temperature-related shifts in these boundary conditions have yet to be assessed. With this in mind we provide a new evaluation of the temperature sensitivity of both reaction kinetics and diagenetic boundary conditions in an attempt to move forward our understanding of the net effect of both past and future warming episodes on marine P burial.

In this study, we expand upon a comprehensive benthic diagenetic model (the SedChem model; Tarhan et al., 2021; Zhao et al., 2020) to provide a new look at the processes responsible for climatic shifts in phosphorus burial in continental margin sediments. Specifically, we incorporate current understanding of the temperature sensitivity of a variety of factors such as reaction rates, solubility product (Ksp), activity coefficients, diffusion coefficients, dissociation constants, seawater pH, seawater dissolved oxygen concentrations and bioturbation. With this modeling exercise, we aim to provide a more complete understanding of the compounded effects of temperature change on both phase-specific (organic P, CFA) and total reactive P burial.

2. Method

We adapt a recently developed diagenetic model (SedChem; cf. Zhao et al., 2020) to explore the impact of temperature change on benthic biogeochemical cycling. This diagenetic model is a one-dimensional multicomponent model with a domain thickness of 300 cm. It includes full reaction-based representation of biogeochemical cycling for C, N, P, S, Fe and Mn. Modifications to SedChem for this study include parameterization of the impact of temperature on reaction rates, Ksp, activity coefficients, diffusion coefficients, dissociation constants, seawater pH, seawater dissolved oxygen concentrations and bioturbation. The range of warming explored in this study is 0–8 oC, as we have limited knowledge of the temperature responses, beyond this range, of several of these parameters, as discussed in greater detail below. This magnitude of warming is more than what has been forecasted for the next century by GCMs but ~50% of the Phanerozoic range (e.g., Guidry & Mackenzie, 2000; IPCC, 2021).

2.1. The effect of temperature on biogeochemical reaction rates

Chemical reaction rates generally increase at higher temperatures. In this study, the effect of temperature change on the rate constant (kT) of a biogeochemical reaction is parameterized using the Arrhenius equation (e.g., Arndt et al., 2013; Burdige, 2011):

kT=keEaRT

where Ea is the activation energy, R is the gas constant, T is temperature in Kelvin and k is the pre-exponential or frequency factor. The activation energy for each reaction is shown in table 1. Note that the collective apparent Ea for microbial processes includes the effect of temperature on both ecological factors and the degradability of organic matter (Arndt et al., 2013).

Table 1.Activation energies for the reactions in the model
Reactions Activation energy (kJ mol-1) References
R1 30-130 Burdige, 2011
R2 30-130 Burdige, 2011
R3 30-130 Burdige, 2011
R4 30-130 Burdige, 2011
R5 30-130 Burdige, 2011
R6 30-130 Burdige, 2011
R7 70 Strous et al., 1999
R8 95.81 Zhang et al., 2002
R9 29 Millero, Sotolongo, et al., 1987
R10 41.6 Chiriţă et al., 2008
R11 88 Nicholson et al., 1988
R12 46 Millero, Hubinger, et al., 1987
R13 52.7 King & Adamsen, 1992
R14 50 This study
R15 14 Yao & Millero, 1993
R16 4.5 Yao & Millero, 1996
R17 40 Harmandas & Koutsoukos, 1996
R18 50 This study
R19 104 Licht & Davis, 1997
R20 71.13 Rickard, 1975
R21 56 Yee et al., 2006
R22 50 This study
R23 71.2 Romanek et al., 2011
R23- 28 Gutjahr et al., 1996
R24 41.84 Sanjuan & Girard, 1996
R24- 33.47 Sanjuan & Girard, 1996
R25 70 This study
R25- 53.33 Wang et al., 2017
R26 47 Van Cappellen & Berner, 1991
R27 35 Rickard, 1997
R28 51 Hu & Jun, 2012; McMaster et al., 2008
R29 10 This study
R31 30 This study
R32 30 This study

Note: R23-, R24- and R25- represent carbonate dissolution.

Previous work, however, has indicated that changes in sulfate reduction rates with temperature do not follow the Arrhenius equation when the temperature is near to or higher than the optimum temperature of sulfate reducers (fig. 1; Finke & Jørgensen, 2008; Sawicka et al., 2012). In this instance, the temperature response of sulfate reduction can be described using a peaked function that includes the effect of temperature on enzyme activity (Alexandrov & Yamagata, 2007):

kT=kηf(T)(η1)+fη(T)

f(T)=eEa(TT0)TRT0

Where η is the ratio of levels of activation energy above and below an enzyme’s temperature (the optimum temperature), f(T) is the normalized Arrhenius equation and T0 is the temperature at which f(T) is equal to one.

Figure 1
Figure 1.The temperature responses of sulfate reduction rate (SRR) and bioturbation.

Temperature (°C) versus A. SRR (nmol cm-3 d-1); B. biodiffusion (Db; cm2 yr-1); and C. bioirrigation (Vb; cm yr-1). The in situ seasonal mean temperatures from the corresponding sites are also indicated by the arrows along the upper x-axis. The data for SRR are from Wadden Sea Station WW, Germany (Finke & Jørgensen, 2008); while Db and Vb are from the lower St. Lawrence Estuary, Québec, Canada (Ouellette et al., 2004). Empirical data are denoted by circles; the black and gray curves denote model outputs; simulations for the latter used equations (2) and (3). The shadings represent uncertainties in the estimate of maximum SRR (±10%), biodiffusion (±20%) and bioirrigation (±15%).

However, below this optimum temperature, the relationship between sulfate reduction rate and temperature can be approximately described using the Arrhenius equation (fig. 1, Finke & Jørgensen, 2008; Sawicka et al., 2012). Moreover, incubation experiments conducted under varying temperatures indicate that the optimum temperature for sulfate reducers is usually >8 oC higher than the in situ temperatures of surficial sediments at any water depth and latitude in the modern ocean (Finke & Jørgensen, 2008; Knoblauch & Jorgensen, 1999; Sawicka et al., 2012). On this basis, the Arrhenius equation is likely sufficient to describe the effect of temperature on microbial processes such as organic matter decomposition and sulfate reduction in surficial sediments, as we only explore a temperature change of 8 oC in this study. Therefore, for simplicity, we adopt an Arrhenius framework for the purpose of modeling the impact of temperature on benthic diagenetic processes under the range of boundary conditions considered for this study.

Previous work has indicated that rates of organic matter decomposition generally correspond to both age and temperature. We have used the same activation energy for all modeled pathways of organic matter decomposition, including aerobic respiration, nitrate reduction, manganese reduction, iron reduction, sulfate reduction and methanogenesis. Studies from coastal sediments suggest that the Ea of organic matter decomposition and sulfate reduction is in the range of 25 to 130 kJ mol-1 (fig. 1; Aller, 1980; Burdige, 2011; Klump & Martens, 1989; Knoblauch & Jorgensen, 1999; Middelburg et al., 1996; Robador et al., 2009; Sawicka et al., 2012; Weston & Joye, 2005; Westrich & Berner, 1988), which increases as organic matter reactivity decreases (Middelburg et al., 1996; Westrich & Berner, 1988).

For this study, we adopt this same range of Ea for organic matter decomposition. It is worth noting, however, that the effect of temperature on organic matter degradation may follow a complex pattern at different timescales, due to factors such as physiological adaptation and shifts in microbial community composition (e.g., Arndt et al., 2013; Arnosti et al., 1998; Hubert et al., 2009, 2010; Robador et al., 2009; Weston & Joye, 2005). Long-term temperature changes are likely to induce changes to the composition of sedimentary microbial communities (e.g., Hubert et al., 2009, 2010; Robador et al., 2009). While the effect of changes in microbial community makeup on rates of organic matter decomposition remains poorly understood, a microbial community characterized by a higher optimum temperature, selected under elevated ambient temperatures, may be expected to broaden the temperature range that should be included in the Arrhenius equation.

Notably, previous experimental studies show that the precipitation of apatite is also a function of temperature (e.g., Jahnke, 1981; Van Cappellen & Berner, 1991). The activation energy used in this study for CFA formation (47 kJ mol-1) is taken from Van Cappellen and Berner (1991). The temperature dependence of vivianite—a common P mineral in marine sediments with anomalously high reactive iron concentrations (e.g., the Baltic Sea; Reed et al., 2016)—formation rates is not well known. We have assumed a lower-end activation energy for vivianite formation (30 kJ mol-1).

2.2. The effect of temperature on bioturbation

Temperature has a strong effect on the metabolic and ecological activities of animals—including bioturbating animals—from respiration to locomotion and feeding (Gerino et al., 1998; Ouellette et al., 2004; Sun et al., 1991; Teal et al., 2008). Seasonal variability in temperature, among other factors, may strongly influence the ecological scope and activity of bioturbating animals. For instance, seasonal variability in the distribution of chlorophyll-a in sediments of Long Island Sound, along the eastern seaboard of the United States, suggests similar variability in average intensities of biodiffusion (DB), which are higher in summer than in winter and spring (Sun et al., 1991). Similarly, flume tank experiments using Long Island Sound sediments have indicated the seafloor erodibility—linked to more intensive bioturbation—is greater during the summer months (Rhoads et al., 1978; Yingst & Rhoads, 1978). Further in situ and laboratory measurements indicate that temperature has a strong impact on both biodiffusion and bioirrigation (Gerino et al., 1998). Experiments on the nereid polychaete Alitta virens from northern temperate latitudes at different in vitro temperatures also indicate that both biodiffusion and bioirrigation change with temperature (fig. 1; Ouellette et al., 2004). Results from 30-day experimental incubations with A. virens reveal that the intensity of biodiffusion increases nearly exponentially from lower temperatures characteristic of winter (~1 oC) to a peak during higher temperatures (~13 oC) characteristic of summer, followed by a notable decrease in biodiffusion intensity from 13 to 18 oC (Ouellette et al., 2004). The intensity of bioirrigation, conversely, appears to increase more gradually between 1 oC and 18 oC (Ouellette et al., 2004).

We parameterized the influence of temperature on biodiffusion and bioirrigation using the experimental data of Ouellette et al. (2004), (fig. 1), with the application of equations (2) and (3). We used this approach to calculate annual mean values for the intensity coefficients of both biodiffusion and bioirrigation (Db and Vb, respectively) through the application of a seasonal temperature cycle representing the northern temperate latitudes (fig. 2A; Ouellette et al., 2004). We calculated the potential change in annual mean values of Db and Vb at temperatures that are up to 8 oC higher than the temperature of modern bottom-water temperatures (<100 m water depth) at northern temperate latitudes (fig. 2). In this calculation, the temperature shift was applied to every day of a year. Additionally, we assumed that this relationship can be transferred to other continental margin sites. Our calculated annual mean Db values increase as temperatures increase up to 2 oC higher than modern annual average temperatures, but do not significantly change when temperatures are increased by 2 to 8 oC above this baseline (fig. 2). We interpret this to reflect infaunal biodiffusers approaching their thermal optimum (i.e., maximum intensities of biodiffusion). In contrast, calculated annual mean Vb values gradually increase in concert with an 8 oC increase in temperature (fig. 2), due to the higher optimum temperature of bioirrigators (~20 oC) relative to that of biodiffusers (~13 oC) (Ouellette et al., 2004).

Figure 2
Figure 2.Modeled temperature responses of annually averaged sulfate reduction rate (SRR), biodiffusion intensity (Db) and bioirrigation intensity (Vb).

A. Seasonal temperature variation for the two sites in fig. 1. B. Temperature responses of annually averaged SRR. C. Temperature responses of annually averaged Db. D. Temperature responses of annually averaged Vb. Red, Wadden Sea Station WW; Black, the lower St. Lawrence Estuary. The shadings in B, C and D represent uncertainties in the estimate of maximum SRR (±10%), biodiffusion (±20%) and bioirrigation (±15%). The red curve in B represents a temperature response following the Arrhenius equation. ∆T means a temperature increase at the corresponding site. The original data for each corresponding site can be found in figure 1.

We further applied the parameterized relative changes in Db and Vb in this predicted regional pattern to our diagenetic model. Note that the term Vb (in units of cm yr-1) is used in Ouellette et al. (2004) to represent the intensity of bioirrigation following a diagenetic model from Gerino et al. (1994); whereas the term a0 (in units of yr-1) is used in our model to represent the intensity of bioirrigation following the diagenetic model from Boudreau (1997). These two terms are roughly scaled to each other (a0~Vb/∆z, where ∆z is the distance between the modeled solute and the sediment-seawater interface), so we have directly applied the relative change of Vb with temperature to that of a0 in this study. It is worth noting that, due to limited empirical data on bioturbation-temperature feedbacks, the extent to which the scaling we employ can be readily applied widely across the global ocean remains undetermined and should be a key target for future studies.

2.4. The effect of temperature on the other parameters

We have also simulated the effect of temperature on the Ksp of CFA, activity coefficient of CFA components, diffusion coefficient and dissociation constant (fig. 3, and fig. S1). Previous studies have indicated that the Ksp of apatite decreases with increasing temperature (Harouiya et al., 2007; Van Cappellen & Berner, 1991). We used a linear relationship to describe the scaling between log(Ksp) and temperature (T):

log(Ksp)=log(Ksp0)0.04434(TT0)

where Ksp0 is the reference solubility product at T0, which is 12 oC in this study. Ksp0 is a function of carbonate activity (Jahnke, 1984; Zhao et al., 2020). The slope 0.044 was calculated from the data presented in Van Cappellen and Berner (1991). The relationship between temperature and the activity coefficients of [Ca2+], [Na+], [Mg2+], [CO32-], [PO43-] and [F-] were obtained from Millero and Pierrot (2002). The temperature responses of solubility products of calcite and aragonite, density of seawater, and the dissociation constants of H2CO3, HCO3-, H3PO4, H2PO4-, HPO42-, H2S, NH4+ and H2O were calculated using the R package “seacarb” (Gattuso et al., 2015), and the diffusion coefficients of dissolved components were calculated using the R package “marelac” (Soetaert, Petzoldt, & Meysman, 2010).

High atmospheric pCO2 linked with a warmer climate will induce a drop in seawater pH (e.g., Caldeira & Wickett, 2003). Although our previous work has suggested that bottom-water pH will, under near-baseline conditions, have a relatively limited influence on P burial (Zhao et al., 2020), we consider the effect of temperature on pH, pH = pH0 0.056*(TT0), where pH0 is the initial pH at T0 (fig. 3). This is a relationship that we derive from simulations using the model LOSCAR (Zeebe, 2012).

Figure 3
Figure 3.Temperature response of individual environmental parameters for the temperature-sensitive model runs.

A. pH of marine bottom waters, adopted from LOSCAR model results (see Zeebe, 2012); B. coefficient of biodiffusion intensity (Db); C. KspCFA0, note that KspCFA=KspCFA0102.3307log([CO23]rCO3); D. dissolved oxygen concentration of marine bottom waters ([O2]BW); The black line represents the temperature effect on oxygen solubility, whereas the blue line represents the inclusion of a temperature effect on oceanic anoxia through enhanced P weathering; e. coefficient of bioirrigation intensity (a0). See Supplementary figure 1 for the temperature response of activity coefficients of [Ca2+], [PO43-], [CO32-], [F-], [Na+] and [Mg2+] in pore waters and the influxes of iron and manganese phases for the temperature-sensitive model runs. ∆T means a global temperature increase relative to modern values. The unit μM represents μmol/L.

The effect of temperature on seawater dissolved oxygen concentrations was calculated using the R package “marelac” (fig. 3; Soetaert, Petzoldt, & Meysman, 2010). Note that this only includes a temperature effect on oxygen solubility. To include the effect of temperature on P weathering and its feedback on ocean anoxia (e.g., Komar & Zeebe, 2017; Van Cappellen & Ingall, 1994, 1996), we also explore a scenario in which oxygen levels decline more quickly with an increase in temperature (i.e., a linear decrease in dissolved oxygen levels from 150 μM to 25 μM driven by a temperature increase from 12 oC to 20 oC).

The link between climate, primary productivity and carbon loading to sediments is a controversial topic. Further, we would expect this relationship to vary between different warming events—and to be highly spatially variable. Moreover, the relative impact of this feedback would likely be timescale dependent. The direct effect of temperature on primary productivity may occur on seasonal to annual timescales. However, the impact of temperature on primary productivity through P weathering is likely to only become important over long (10–100 kyr) timescales. Nonetheless, there is some basis for exploring this relationship. Observations across the modern ocean indicate a lack of clear directionality in the relationship between marine net primary productivity (NPP) and sea surface temperature (Behrenfeld et al., 2006; Boyce et al., 2010; Vancoppenolle et al., 2013). Model results from CMIP5 suggest that marine NPP will decrease by 3–10% by 2100, whereas CMIP6 results show substantially less reduction of marine NPP by 0–3% (Kwiatkowski et al., 2020; Moore et al., 2018). During climate warming events in the geological past, such as the PETM, indirect effects of temperature on primary productivity, such as enhanced P weathering and P regeneration, may have boosted primary productivity (Komar & Zeebe, 2017; Zhao, Mills, et al., 2023). However, the combined effects of these process on primary productivity are still within the range of ±30% (Komar & Zeebe, 2017). As such, we applied a range of values in organic carbon loading to sediments (0-30% variation over 8 ˚C warming, with organic matter loading linearly scaled with the extent of warming).

2.5. Other modifications and model solution

The SedChem model was modified to include variations in the C/P ratio of organic matter as driven by redox conditions. This is achieved through the addition of a scaling factor that enhances P regeneration (fanP= 1.6) via anoxic pathways for organic matter decomposition such as nitrate reduction, Mn reduction, Fe reduction, sulfate reduction and methanogenesis (Dale et al., 2016). The C/P ratio of all organic matter delivered to sediments is set at 106/1. We have also added a scaling factor to represent the enhanced rate of organic matter decomposition under oxic conditions (foxC= foxP= 2) (Katsev & Crowe, 2015). The rate of CFA formation is represented by the following expression:

RCFA=kCFAH0.35eEaCFART(ΩCFA1/τ1)

where kCFA is the rate constant for CFA formation (0.012 mol L-1 yr-1), H is the concentration of H+, EaCFA is the activation energy for CFA precipitation, ΩCFA is the saturation state of CFA and τ (τ=18.87) is the number of ions per formula unit of CFA (Van Cappellen & Berner, 1991). Vivianite is also included, following the previous versions of our model (Tarhan et al., 2021; Zhao et al., 2020).

As described in greater detail in Zhao et al. (2020) and Tarhan et al. (2021), the model was written in R (R development core team, 2006). The R package ReacTran was used for the transport terms (Soetaert & Meysman, 2012). The model was solved via the Method of Lines using the Variable Coefficient Ordinary Differential Equations (VODE) solver (Soetaert, Petzoldt, & Setzer, 2010). All model simulations were run to steady state. Other than the temperature effects of parameters explored in this study, all the other parameters are the same as in the baseline model for shallow seawater in Zhao et al. (2020).

3. Model results and discussion

3.1. Full temperature-sensitive results and sensitivity tests

We performed a series of runs to simulate the influence of warming temperatures associated with temperature change on marine phosphorus burial in continental margin sediments (figs. 4 and 5). As described above, we have included the effect of temperature on a variety of factors such as reaction rates, Ksp, activity coefficients, diffusion coefficients, dissociation constants, seawater pH, seawater dissolved oxygen concentrations and bioturbation. Following the approach of previous studies (Tarhan et al., 2021; Zhao et al., 2020), we have calibrated the SedChem model using empirical pore water and sediment data from the particularly well-characterized shallow marine FOAM (Friends of Anoxic Mud) site in Long Island Sound, along the eastern seaboard of the USA, and which is characterized by a heterolithic lithology (sand, silt, mud), as is typical of shallow-water continental-margin settings (e.g., Aller, 1977; Berner et al., 1978; Canfield, 1988; Ruttenberg, 1990). A revised version of SedChem, calibrated to include temperature responsivity for a range of variables and parameters, was used to investigate the influence of temperature on marine P burial along continental margins. In these runs (hereafter referred to collectively as full temperature-sensitive runs), temperature is allowed to vary from modern temperatures observed at the FOAM site (12 oC) to up to 8 oC higher (i.e., 20 oC), and drives temperature-sensitive responses of biogeochemical reactions, diffusion coefficients, dissociation constants and the other forcings (table 1, fig. 3 and fig. S1). An upper limit of 8 oC warming is utilized here due to the limited nature of empirical data constraining the temperature response of bioturbation beyond this range. The maximum temperature in the bioturbation dataset we draw upon (that of Ouellette et al., 2004) does not exceed 8 oC warmer than the modern summer temperature (fig. 1). The results of our temperature-sensitive modeling exercise demonstrate that both CFA content and P burial increase with increasing temperature when a temperature effect on oxygen solubility is included (fig. 4). Specifically, sedimentary CFA concentrations and the P burial flux are elevated by ~40% and ~20%, respectively, following exposure to 8 oC of warming above baseline conditions. This is mainly due to the positive relationship of temperature and CFA formation rate. However, in the scenario with a significant oxygen drop due to temperature-enhanced P weathering and its feedback on oceanic anoxia, the positive effect of temperature on CFA formation is strongly dampened by oceanic anoxia (fig. 4). Given that our simulations were run with a relatively small reactive iron load (as is characteristic of most modern normal marine settings), vivianite does not, in our model outputs, play an important role in controlling P burial rates.

Figure 4
Figure 4.Modeling results exploring the temperature sensitivity of marine P burial, with the inclusion of temperature effects on all environmental parameters explored in this study.

A. CFA concentration in the final burial sediments; B. organic phosphorus concentration in the final burial sediments; C. total phosphorus burial flux. We examine a range of activation energies for organic matter decomposition, including 80 kJ mol-1 (solid curves), and 30 kJ mol-1 to 130 kJ mol-1 (shaded regions). The black curves represent the temperature effect on oxygen solubility, whereas the blue curves represent the inclusion of a temperature effect on oceanic anoxia through enhanced P weathering.

Figure 5
Figure 5.Modeling results exploring the effect of temperature on marine P burial, with the inclusion of temperature effects on each environmental parameter.

A., D., G. CFA concentration in the final burial sediments; B., E., H. organic phosphorus concentration in the final burial sediments; C., F., I. total phosphorus burial flux. We examine the effect of temperature sensitivity of CFA precipitation rate (A–C), bioturbation (D–F) and dissolved oxygen levels (G–I) on P burial. The black curves in (G–I) represent the temperature effect on oxygen solubility, whereas the blue curves represent the inclusion of a temperature effect on oceanic anoxia through enhanced P weathering.

We also present sensitivity tests in which there is a temperature effect on only CFA precipitation rate, bioturbation or dissolved oxygen levels, individually, in order to examine the strength of each of these forcings (fig. 5). We elected to explore bioturbation and dissolved oxygen levels in these sensitivity analyses because these factors have been previously demonstrated to be sensitive to temperature variation (e.g., Ouellette et al., 2004; Soetaert, Petzoldt, & Meysman, 2010) and previous modeling work has indicated they are among the most important factors regulating the magnitude of CFA burial (Zhao et al., 2020).

In the sensitivity analyses of CFA precipitation rate, the Ksp and the rate constant of CFA, and the activity coefficients of [Ca2+], [Na+], [Mg2+], [CO32-], [PO43-] and [F-] were allowed to vary jointly with temperature. In this simulation, we observed an increase in temperature-enhanced CFA formation and thus total marine P burial flux within the sediment (fig. 5A-C). The magnitude of CFA formation under these conditions (a ~25% increase under an 8 oC temperature increase) and total marine P burial flux (a ~16% increase associated with a temperature increase of 8 oC) in this simulation is similar to the full temperature-sensitive run with the temperature effect on oxygen solubility discussed above. This indicates that the temperature responsivity of CFA precipitation rate is the most important factor regulating the positive correlation between temperature and CFA and total P burial, as observed in the full temperature-sensitive run that includes representation of the effect of temperature on oxygen solubility.

In contrast, simulations in which only the temperature sensitivity of bioturbation was considered suggest that temperature-mediated shifts in the intensity of bioturbation (with biodiffusion and bioirrigation considered jointly) have a limited effect on CFA and total P burial (fig. 5D-F). This is likely generated by competing effects of biodiffusion and bioirrigation on P burial (see Tarhan et al., 2021).

Depending on the resulting magnitude of change in dissolved oxygen levels, this will have either a small or substantial influence on sedimentary CFA concentrations (fig. 5G-I). When sensitivity of dissolved oxygen levels to temperature is included in model simulations, the amount of CFA burial decreases (fig. 5G). This is due to the inverse relationship between temperature and dissolved oxygen and the positive effect of bottom-water oxygen on CFA burial (Tarhan et al., 2021; Zhao et al., 2020). When only this temperature effect on oxygen solubility is included (and other temperature sensitivities are excluded), the shift in CFA formation with temperature is small.

However, when oxygen levels drop significantly at higher temperatures (i.e., a linear decrease in dissolved oxygen levels from 150 μM to 25 μM driven by a temperature increase from 12 oC to 20 oC) due to the increase in the relative importance of temperature-enhanced P weathering (which can, in turn, drive an increase in ocean anoxia on long (10–100 kyr) timescales), sedimentary CFA abundances decrease dramatically in response to a rise in temperature (8 oC warming drives a ~25% decrease in CFA concentration). In this simulation, we have applied an endmember relationship representing a significant shift in ocean oxygen level with temperature. However, the actual shifts in oceanic oxygen levels during climate warming events are likely characterized by pronounced spatial heterogeneity, and may be influenced by additional factors such as the background climate state, the weatherability of P on land and the configuration of continents and oceanic basins (e.g., Pohl et al., 2022).

Our previous work has shown that seawater Ca concentration substantially influences CFA burial (Zhao et al., 2020). Thus, in order to investigate the temperature response of P burial predicted for geological periods characterized by evidence for high seawater Ca concentrations (cf. Horita et al., 2002), we also performed a set of runs with boundary conditions characterized by high bottom-water Ca concentrations (fig. S2). In this set of runs, seawater Ca concentration is not a function of temperature. Under high seawater Ca concentration, the magnitude of shifts in CFA concentration and total P burial flux do not change substantially relative to our temperature-sensitive model runs employing boundary conditions characteristic of the modern ocean (fig. S2). However, the relative shifts in CFA concentration and total P burial flux are substantially smaller under a scenario of high seawater Ca concentration (fig. S2). This suggests that the temperature response of P burial could be significantly muted during geological intervals characterized by relatively high seawater Ca concentrations.

3.2. Monte Carlo simulations

Given the uncertainties tied to the temperature sensitivities of bioturbation and sedimentary organic matter loading, we additionally performed Monte Carlo simulations (fig. 6). The range of parameters explored is presented in Figure S3, and a uniform distribution is used to sample the parameters. The results show that the directionality of change in CFA concentration is variable. As discussed above, this likely reflects the competing effects of temperature and development of anoxic conditions on CFA precipitation.

Figure 6
Figure 6.Monte Carlo results exploring the impact of temperature on marine P burial.

A. CFA concentration in the final burial sediments; B. organic phosphorus concentration in the final burial sediments; C. total phosphorus burial flux. The relationships between temperature and environmental parameters for the Monte Carlo simulation that are different from the temperature-sensitive model run are shown in fig. S2. The number of runs is 10,000 for the Monte Carlo simulation. The vertical colored bar indicates the density of results.

3.3. Implications for the global P cycle

The effect of temperature on weathering of P from the continents has been relatively well-established (Guidry & Mackenzie, 2000, 2003) and has been applied to simulate P cycling during the Phanerozoic, including past extreme warming events such as the PETM (e.g., Guidry & Mackenzie, 2000; Komar & Zeebe, 2017). It has been widely accepted that increases in P weathering should boost primary productivity, resulting in increased ocean anoxia as well as, via anoxia-driven positive feedback mechanisms, enhanced P regeneration from sediments on long (10–100 kyr) timescales (e.g., Komar & Zeebe, 2017; Van Cappellen & Ingall, 1994, 1996). Our results, however, suggest that warming may also stimulate CFA formation by enhancing CFA precipitation rates, which may in turn mute expected increases in benthic recycling of P from organic matter during climate warming events, as a greater proportion of organic P will instead be transformed to CFA through “sink switching” in seafloor sediments. For instance, in the simulation that includes representation of the effect of temperature on oxygen solubility (fig. 4), the concentration of CFA increases by 2 μmol g-1, although the concentration of organic P decreases by 0.8 μmol g-1. However, where increases in productivity driven by increased P weathering fluxes in a warmer climate state drive increased ocean anoxia, this may ultimately limit the extent to which warming leads to increased CFA formation.

Temperature change in the deep oceans will likely keep pace with temperature changes on land over timescales exceeding oceanic mixing times (~2 kyr). On these thousand-year and longer timescales, the effect of temperature on P weathering fluxes (Guidry & Mackenzie, 2000, 2003) could, in theory, be similar to its effect on CFA burial in marine sediments if there is no significant development of oceanic anoxia, as the activation energies for apatite dissolution and formation are likely to be roughly similar. However, if significant oceanic anoxia develops, the direct temperature effect on CFA precipitation rate would likely be strongly buffered by the effect of temperature on CFA burial by means of enhanced P weathering and oceanic anoxia. Further, the shift in the extent of P solubilization during weathering will be dependent on the changes in vegetation and hydrologic cycling with warming.

While we have not attempted to reproduce changes in P burial across multiple regions and settings—given that each site would be characterized by distinct boundary conditions (including temperature shifts) and sedimentary P compositions—our results do provide a general first-order estimate of the effect of climatic shifts (and temperature change in particular) on P burial. More spatially resolved modeling of P burial during warming—where boundary conditions vary—would be an obvious follow-up to this work. In many ways, this work highlights the need for additional empirical and modeling work on temperature-driven effects on benthic P cycling. As noted above, we have not explored reactive iron-rich sediments, where vivianite formation is likely to become more important than apatite formation (e.g., Reed et al., 2016). Nonetheless, our results highlight that there are multiple positive and negative feedbacks influencing phosphorus burial in continental margin sediments during warming events and suggest that increased apatite formation during warming may be an important but underappreciated aspect of the response of the P cycle to marine warming.

4. Conclusion

Using a benthic biogeochemical (diagenetic) model with a comprehensive reaction network, we have provided a first-order estimate of the effect of temperature change on marine P burial. We show that warming can drive increased marine phosphorus burial in continental margin settings, by means of the strong positive relationship between temperature and the biogeochemical reaction rates associated with authigenic P formation. The temperature effect of P burial observed in our model results is consistent with the observed increases in phosphorite burial during the last interglacial period both in the Arabian Sea and along the Peru margin (Burnett, 1980; Ganeshram et al., 2002; Garrison & Kastner, 1990; Schenau et al., 2000; Shimmield & Mowbray, 1991). Moreover, when we consider the temperature sensitivity of a range of biogeochemical and environmental variables highlighted by previous studies as important modulators of marine P cycling, we observe that, under a scenario of strong oceanic anoxia under higher temperatures (i.e., a >80% decrease in oxygen levels), the impact of warming as a driver of enhanced marine phosphorus burial in continental margin settings will be strongly muted by oceanic anoxia generated by enhanced P weathering. Increases in authigenic P burial during climate warming episodes in Earth’s past (and future) may, under certain boundary conditions, have dampened P regeneration from organic matter, impacting extents of oceanic anoxia and organic matter burial. We suggest that the temperature response of the authigenic P flux should be included in future biogeochemical models simulating marine nutrient cycling across climate change episodes.


ACKNOWLEDGMENTS

M.Z. is funded by the programme of the Chinese Academy of Sciences (E251520401) and the IGGCAS Key programme (no. IGGCAS-202201). We thank Jonathan Payne and two anonymous reviewers for critical feedback that improved this manuscript.

SUPPLEMENTARY INFORMATION

https://earth.eps.yale.edu/~ajs/SupplementaryData/2023/Zhao

https://data.mendeley.com/datasets/bcz4hgd9j4/1

Editor: C. Page Chamberlain, Associate Editor: Jonathan Payne