A tectonically driven Ediacaran oxygenation event

The diversification of complex animal life during the Cambrian Period (541–485.4 Ma) is thought to have been contingent on an oxygenation event sometime during ~850 to 541 Ma in the Neoproterozoic Era. Whilst abundant geochemical evidence indicates repeated intervals of ocean oxygenation during this time, the timing and magnitude of any changes in atmospheric pO2 remain uncertain. Recent work indicates a large increase in the tectonic CO2 degassing rate between the Neoproterozoic and Paleozoic Eras. We use a biogeochemical model to show that this increase in the total carbon and sulphur throughput of the Earth system increased the rate of organic carbon and pyrite sulphur burial and hence atmospheric pO2. Modelled atmospheric pO2 increases by ~50% during the Ediacaran Period (635–541 Ma), reaching ~0.25 of the present atmospheric level (PAL), broadly consistent with the estimated pO2 > 0.1–0.25 PAL requirement of large, mobile and predatory animals during the Cambrian explosion.

T he oxygenation of the Earth system was a necessary condition for the rise of complex animal life [1][2][3][4][5] , which occurred in several steps [6][7][8] . The Great Oxidation Event (GOE)~2. 3 Ga saw a permanent rise in atmospheric oxygen from 9 <10 −5 PAL to 10-13~1 0 −4 -10 −1 PAL, but oxygen remained well below present levels throughout the Proterozoic 12,13 . A Neoproterozoic Oxygenation Event 14 has been proposed based on a range of indirect proxies (Fig. 1). The expansion of the oceanic inventories of Mo, V and Re, and a shift towards lower δ 82 / 76 Se values suggest a trend towards more oxidising ocean conditions 6,[15][16][17][18][19] across the Neoproterozoic-Phanerozoic transition, and cerium anomalies 20 point to at least a well-oxygenated shallow ocean environment after 551 Ma 21 (Fig. 1b). However, the redox state of the ocean clearly fluctuated, with a series of transient oxygenation events (Fig. 1c) getting somewhat more frequent through the Ediacaran and early-mid Cambrian 16 . These include partial and temporary oxygenation of deeper waters following the Sturtian 22~6 60 Ma, Marinoan 23~6 35 Ma and Gaskiers 24,25~5 80 Ma glaciations.
Several hypotheses have been put forward to explain a Neoproterozoic oxygenation, in which the observed ocean oxygenation trends are the result of a broader increase in atmospheric O 2 . Following the Snowball Earth glaciations, pulses of nutrients are suggested to have entered the oceans, enhancing primary production and burial of organic matter, releasing oxygen to the atmosphere 24,26,27 , but this would only have temporarily increased pO 2 28 (returning to the previous state after the pulse subsided). Alternatively, a sustained increase in terrestrial chemical weathering 29 , potentially amplified through selective weathering of phosphorus by early terrestrial ecosystems 30 and the presence of P-rich large igneous provinces 31 could have increased pO 2 . The expansion of an early land biosphere could also have restricted the oxidative weathering of reduced crustal rock 32 , reducing the major sink of O 2 . These mechanisms can increase atmospheric O 2 concentrations, but all also imply a rise in the average δ 13 C of carbonates, either by increasing organic carbon burial relative to carbonates or by restricting organic carbon weathering, and such a rise is not observed in the geological record at the time.
Here we explore an alternative mechanism of oxygenation, where changes in plate tectonics cause a rise in pO 2 during the late Neoproterozoic. An increased fraction of young zircon grains 33 at the Proterozoic-Phanerozoic transition points to increased continental arc volcanism, consistent with an increase in atmospheric CO 2 and warming of the planet from Cryogenian extreme glaciations to ≥20°C 34 during the Cambrian. Continental volcanic arc extent varies through time, broadly matching the detrital zircon age data 35 . Moreover, assuming that global arc magmatism is proportional to subduction, mantle depletion curves also point to a maximum in subduction during the early Phanerozoic 36 , supported by reconstructions of subduction-zone lengths from both plate-tectonic reconstructions 37 and kinematic modelling 38 . Whilst these methods all carry uncertainties, the geologic data all point to a general step-increase in subductiondriven degassing across the Proterozoic-Phanerozoic transition.
As there is no direct geochemical proxy for atmospheric oxygen levels, estimates largely rely on modelling approaches [39][40][41] , based on our understanding of the long-term carbon cycle (Fig. 2). For oxygen to build up in the ocean-atmosphere system, photosynthetically-produced organic carbon must be shielded from being re-oxidised through the burial of organic carbon in seafloor sediments. Inorganic carbonates are also buried in sediments following precipitation from seawater, and both buried carbon species are returned to the ocean and atmosphere through uplift and weathering, as well as via metamorphism and degassing. Any change to the burial and return fluxes of organic carbon will cause a change in atmospheric oxygen concentration.
The transformation of CO 2 to organic carbon via photosynthesis imparts a significant isotopic fractionation. Buried organic carbon contains a reduced fraction of 13 C atoms, leaving an excess in the ocean and atmosphere. Through this mechanism, changes in the rate of burial and weathering of organic carbon relative to the burial and weathering of carbonates alter the δ 13 C value of seawater, which is recorded in sedimentary carbonates. This effect is dampened at low pO 2 by the rapid adjustment of oxidative weathering rates to match carbon burial rates 12 (through changes in O 2 ), but not completely nullified unless oxidative weathering is the only O 2 sink. Whilst a step change in 82 Stronger negative cerium anomalies recorded in marine cements 20 (Ce anom ; pale blue dots), defined as Ce anom < 1, are indicative of more oxygenated regional to basin-scale conditions (the blue line is a fit through the Ce anom data). Both proxies suggest a trend towards more oxygenated ocean conditions during the Ediacaran period. c A summary of redox sensitive element (RSE; magenta) enrichment data (Mo, U, Re, V, Cr) from black shales, indicating intervals of widespread ocean oxygenation 16,17,22,23 . Below this is a summary of iron-speciation proxy data 24,25 , which records localised redox conditions, hence is divided by depth range. The intervals of the Sturtian, Marinoan and Gaskiers glaciations are also indicated. These redox proxies suggest a series of transient oxygenation events during the Cryogenian and Ediacaran periods carbonate δ 13 C likely occurred between the mid-Proterozoic and Phanerozoic baselines 42 , the later Neoproterozoic record does not show the sustained stepwise increase that would be expected if the long-term burial and weathering fluxes of organic carbon were increased relative to those of inorganic carbon. Leastsquares regression analysis of carbonate carbon isotope data shows a gradual negative trend over the whole Neoproterozoic (decline of~1-2‰), and no significant trend during the Ediacaran period where evidence for deep ocean oxygenation is found (see Supplementary Fig. 1).
Increasing the long-term tectonic input rate of CO 2 must result in an increase in the overall carbon burial rate to maintain steady state of the ocean-atmosphere system. If some of this carbon is buried organically then O 2 can rise, and if this additional burial follows the same distribution between organic carbon and carbonates as in the initial unperturbed system, then the oxygen concentration of the atmosphere can rise without any observable change in δ 13 C. This effect should be seen in carbon cycle models such as COPSE 40 and GEOCARBSULF 41 , although it has not been noted in any previous analyses using these models. In addition, a similar effect may have operated over the whole of Earth history, wherein cumulative mantle CO 2 input can drive long-term burial of organic carbon and planetary oxygenation 7 , but the model of this mechanism could not account for the timing or magnitude of any second oxygen rise after the GOE.
Here we investigate the impacts of the inferred increase in tectonic CO 2 input between the late Neoproterozoic and early Phanerozoic 33,35,37 by extending the COPSE Reloaded biogeochemical model 43 (see Methods) to the Ediacaran. Our model predicts a~50% increase in atmospheric O 2 during the Ediacaran. The predicted magnitude and temporal dynamics of changes in redox and 87 Sr/ 86 Sr compare well to available data, as does the lack of a secular trend in carbonate δ 13 C or δ 34 S.

Results
Steady-state computations. Steady-state computations carried out at the Precambrian-Cambrian boundary (Fig. 3) show that the atmospheric oxygen reservoir size in COPSE is indeed strongly influenced by the relative rate of degassing, with higher degassing rates increasing the carbon throughput of the Earth system and atmospheric pO 2 . The burial fluxes of both organic and inorganic forms of carbon increase when degassing is increased, but the fraction of carbon burial that is organic (f org ) remains essentially constant, allowing the long term δ 13 C carb record to remain relatively stable (Fig. 3c). This demonstrates that degassing-driven oxygenation is robust to the multiple biogeochemical feedbacks that are present in the COPSE model, which includes the consideration that a fraction of subducted material is organic carbon, and thus will consume oxygen when recycled. Nevertheless, in order to test our hypothesis fully, it is important to view the whole spectrum of model outputs. The COPSE model produces estimates for δ 13 C carb , δ 34 S seawater and 87 Sr/ 86 Sr, and the current model baseline reproduces these records reasonably well over the Phanerozoic, providing 'ground-truthing' that the model processes are a sensible representation of global biogeochemistry 43 . If our hypothesis for degassing-driven oxygenation is reasonable, the model should be able to reproduce the key longterm trends in these isotope records over the Ediacaran period when it is subject to our proposed increase in degassing rates alongside other expected external forcings.
We first assess the model steady states with respect to changes to both degassing rate and the other major tectonic forcing in the model-the rate of uplift and erosion (Fig. 4). Uplift and erosion enhances continental weathering fluxes in the model, which increases phosphorus delivery to the ocean and stimulates organic carbon burial. Uplift also exposes fossil organic carbon to oxidative weathering, which is the major sink of O 2 over geological timescales. It is likely that the processes of uplift and erosion accelerated during the late Neoproterozoic in line with the collisions that formed the supercontinents Pannotia and Gondwana: reconstructed sediment abundances increase between the Ediacaran and Cambrian 44 , and rising 87 Sr/ 86 Sr ratios 45 in carbonates can be attributed to the weathering of more ancient crustal material.
In line with Fig. 3, the steady states shown in Fig. 4 show that increased degassing rates result in increased rates of burial of both organic and carbonate carbon, and increase atmospheric O 2 , whilst maintaining an almost-static δ 13 C ratio. Increased degassing drives a small reduction in seawater sulphate δ 34 S, due to a minor decrease in the rate of burial of pyrite sulphur relative to gypsum. In the model, rates of pyrite burial are determined based on the inferred rate of net microbial sulphate reduction, which is assumed to be positively influenced by the availability of organic matter and sulphate but is restricted by expanding deep water O 2 40,46,47 . This formulation produces a reasonable reconstruction of seawater δ 34 S over the Phanerozoic 43 , and in the case of increasing degassing explored here, the positive effects outweigh the negative and pyrite burial increases. However, as with the carbon cycle, burial of the oxidised form of sulphur (gypsum) increases by a similar amount, resulting in only a small change in δ 34 S isotopic fractionation. The 87 Sr/ 86 Sr ratio is decreased under an increase in degassing rates due to the increased input of less-radiogenic mantle-derived Sr.
An increase in uplift rate acts to increase the rates of burial of both organic and carbonate carbon at steady state, due to net increases in carbonate weathering-deposition and greater delivery of phosphorus. But an increase in uplift acts to slightly decrease O 2 overall, due to the enhanced weathering of organic carbon. Uplift and erosion act to decrease δ 13 C carb due to increased overall rates of burial of carbonates derived from carbonate weathering (as noted peviously 48 ). Increases to uplift and erosion act to increase seawater δ 34 S, because both the availability of organic C and the reduction in O 2 favour increased rates of microbial sulphate reduction and subsequent burial of pyrite. Finally, increased erosion rates drive a substantial increase in  The long-term carbon cycle. The cycle is composed of fluxes between atmospheric and oceanic carbon, organic carbon and carbonate. Carbon is moved from the atmosphere and ocean to the crust through burial (B), with weathering (W), degassing and metamorphism (D) offering mechanisms by which carbon is returned to the ocean/atmosphere system. Oxygen sources are displayed as blue, oxygen sinks as orange, and other processes are displayed as grey. Fractionation of buried organic carbon relative to atmosphere and ocean carbon is denoted ΔC. The δ 13 C value of the atmosphere and ocean reservoir can be altered by changing the proportion of organic and inorganic carbon burial, or by changing fractionation effects on either organic carbon or carbonate burial seawater 87 Sr/ 86 Sr due to the increased weathering contributions of granites and carbonates in the model, which are relatively radiogenic.
Monte-Carlo setup. We now run the COPSE model forwards in time for the Ediacaran period and compare to geochemical data for δ 13 C, δ 34 S and 87 Sr/ 86 Sr. The degassing rate follows Mills et al. 37 as in Fig. 1a, and the uplift rate is set to increase linearly from 0.5 to 2 times the present day rate. The increase in uplift is chosen in order to reproduce the magnitude of the observed rise in 87 Sr/ 86 Sr, but is also roughly consistent with the difference in reconstructed rates of sediment deposition between the Ediacaran and Cambrian 44 . In order to fully capture the uncertainty in model parameterisation, we run a Monte-Carlo analysis in which the COPSE model is run 10,000 times and the most important parameters are randomly sampled from an uncertainty window shown in Table 1.
We vary the assumed present-day fluxes in the model carbon cycle because these carry large uncertainties 43 , and will impact the relative effects of an increase in degassing. We vary the activation energy of seafloor weathering, as this process is a carbon sink that responds to increases in mid-ocean ridge production, thus should be increased under our assumption of increasing production and subduction rates. A highly-sensitive seafloor weathering sink would be expected to dampen the O 2 rise caused by an increase in degassing. Continental weathering activation energies are also varied, altering the relative strength of seafloor weathering. In the case of oxidative weathering, we test the range of possible responses to O 2 . We also consider different values for the longterm climate sensitivity. Finally, we experiment with adding a direct mantle flux of CO 2 and H 2 into the model: COPSE does not include direct mantle input at mid-ocean ridges, and the input of reduced H 2 gas is an additional oxygen sink that should be amplified by increasing degassing rates. Ridge CO 2 input is assumed to equal H 2 input (in terms of mol C and mol O 2 equivalent) to maintain balance of the linked carbon and oxygen cycles under an additional O 2 sink, this value is also consistent with measurements 49 . Finally, in addition to the parameter choices, the degassing forcing is resampled every 10 Myrs from the uncertainty window of Mills et al. 37 . This allows for quite rapid changes in degassing rate, but not above the rates of change that have been suggested for more recent time periods 50 .
The Monte-Carlo procedure follows that of Royer et al. 51 , who sampled a wider range of parameters in a model of roughly equal complexity (GEOCARBSULF 41 ), also using 10,000 runs. We choose to use a flat distribution between the parameter minimum and maximum estimates, with all values being equally likely, as we believe this is more representative of the real uncertainty than sampling from a normal distribution. The Monte-Carlo results for the carbon, oxygen and strontium cycles are shown in Fig. 5, with the dark and light shaded areas showing ± 0.5 and ± 1 standard deviation respectively over the whole experimental set.
Model predictions. Early Ediacaran pO 2 levels are predicted to be around~0.2 PAL. This is consistent with recent modelling of Proterozoic oxygen regulation and with constraints on pO 2 from the absence of detrital pyrite 12 . It is also below upper bounds of 0.4 PAL 47,52 or 0.5-0.7 PAL 28 inferred from ocean models, based upon the assumption of present-day nutrient levels and the presence of widespread deep ocean anoxia. Mean atmospheric oxygen is predicted to increase by around 50% during the Ediacaran period. This is primarily due to a rise in the organic carbon burial flux from~4 × 10 12 mol yr −1 to~8 × 10 12 mol yr −1 , which is driven by increasing phosphorus input and ocean phosphate concentration. These increases ultimately stem from the increased CO 2 input flux through degassing and associated increase in surface temperature and weathering, which then drive carbon burial in both organic and inorganic forms. It is possible to draw a line through the model uncertainty that represents no change or even a decrease in O 2 levels. However, this requires moving from one edge of the uncertainty window to the other and is therefore very unlikely. Figure 6 shows a histogram of the modelled change in pO 2, and shows that 97% of model runs produce an increase in atmospheric O 2 , and more than two thirds of runs produce an O 2 increase of between 25 and 75% of the initial value. Consequently, we infer with a high likelihood that there was a significant increase in pO 2 across the Ediacaran period.
Mean estimates in Fig. 5 show a slight reduction in f org and carbonate δ 13 C, whilst the plotted Ediacaran C isotope record shows a series of positive and negative excursions, the causes of which remain the subject of much active research e.g., 53,54 but the long-term average remains constant (see Methods and Supplementary Fig. 1), broadly consistent with our model. Increased weathering of carbonates, and uplift-driven weathering of older  Supplementary Figs. 2-4 and Supplementary Note 1). This results in a much poorer fit to the 87 Sr/ 86 Sr record, but the results for the carbon, sulphur and oxygen cycles are very similar, and as such we conclude that our proposal of degassing-driven O 2 rise is robust to any expected uplift rate changes. Figure 7 shows the sulphur cycle outputs from the Monte-Carlo experiment. With increased inputs of sulphur from continental weathering and degassing, oceanic sulphate concentration, alongside pyrite and gypsum burial rates, also increase  production. Again, as in the carbon cycle, the increase in both reduced and oxidised fluxes to the sediments leads to little change in the fraction of S that is buried as pyrite (f py ), and thus little change to δ 34 S. This is generally consistent with the geologic record, which shows little change in δ 34 S over the Ediacaran 55 , especially when considered against the variability window during the Phanerozoic (~30‰). The record does indicate a shift to slightly higher values towards the Cambrian, which we do not replicate. Nevertheless, a rise in δ 34 S would usually be interpreted as an increase in pyrite burial rates relative to gypsum (or reduction in pyrite weathering), thus it is unlikely that the process driving this change could deplete O 2 . In addition to this, we note that the modelled f py is lower than indicated in some reconstructions 55 . The low f py value is an inherent property of the COPSE model, which assumes relatively high rates of gypsum weathering and burial. The absolute value of f py is highly uncertain 56 , and the important prediction for this study is the change to f py that might be driven by our proposed mechanism, in which our model is in agreement with available data.

Discussion
We can compare our model atmospheric pO 2 predictions to estimates of the oxygen requirements of early animal life forms. This assumes they lived in waters equilibrated with the atmosphere-e.g., in well-mixed (shallow) shelf seas or in benthic slope settings on down-welling margins or high-latitude regions of deep convection 47 . Other locations below the surface mixed layer of the ocean-e.g., seasonally-stratified shelf seas or the open ocean thermocline-tend to be depleted in oxygen due to net respiration of organic matter 47 . It has classically been argued that a minimum oxygen threshold exists for the evolution of animal life 1 , but oxygen requirement clearly depends on the type of animal, including their size, mobility, nervous system (information processing) and ecological habits. Sponges (Porifera)-which are the basal animals 57 -have low pO 2 requirements 58~0 .005-0.04 PAL. Hence their evolution was not limited by any of our predicted pO 2 levels, consistent with biomarker evidence that demosponges were present by~660-640 Ma in the Cryogenian 59,60 . The soft-bodied Ediacaran biota~575-540 Ma are currently interpreted as including a mix of stem-and crown-group animals 61 . The estimated pO 2 requirement of soft-bodied, thin, sheet-like animal forms (e.g., Dickinsonia), which are assumed limited by O 2 diffusion range from~0.01-0.03 PAL 62 to~0.06 PAL 63 . Modern benthic invertebrate analogues suggest a requirement of~0.1 PAL 64 , closer to the levels we predict during the early Ediacaran period. Whilst it is tempting to infer that our predicted steady oxygenation during the Ediacaran period enabled the evolution of progressively more complex Ediacaran animals, pO 2 may already have been sufficient beforehand. Furthermore, once bilaterian animals with a circulatory system evolved they could have tolerated lower pO 2 ; small bilaterian worms with a circulatory system are estimated to only need 0.0014-0.0036 PAL 65 . A better case can be made that our predicted rise in pO 2 enabled the increased size, mobility, nervous system, carnivory 3 Precise quantitative estimates of oxygen levels are difficult to make from a simple box model, in which many parameters are uncertain and one cannot be entirely sure that every relevant process is included. COPSE minimizes this uncertainty as much as possible by comparing its Phanerozoic predictions to multiple whole-Phanerozoic geochemical records 43 . Even so, whilst there remains high confidence in our qualitative result, the quantitative evolution of oxygen levels from~0.2 to~0.3 PAL that we show should be taken as an informed best guess.
A simple, unidirectional marine oxygen rise is overly simplistic, and shorter-term oscillations in the redox state of the ocean, at both local and global scales, clearly occurred during the Proterozoic-Phanerozoic transition 2,8,25,66 . Nevertheless, their increasing frequency 16,20 is consistent with an overall trend of oxygenation due to rising atmospheric pO 2 . To summarise, our predicted long-term Ediacaran oxygenation event driven by increased tectonic degassing is robust to model uncertainties, fits overall trends in geochemical data, and is consistent with existing inferences of pO 2 requirements of the Cambrian fauna. Whilst disentangling the many factors influencing faunal evolution is beyond the realms of this study, we provide the first quantitative prediction of Ediacaran oxygenation that is consistent with geochemical data and with estimated pO 2 requirements for the Cambrian explosion 1,2,5 .

Methods
Regression analysis of carbonate carbon isotopes. Supplementary Fig. 1 shows plots of carbonate δ 13 C for the Neoproterozoic (left, yellow) and the Ediacaran subset (right, red) 67 . For both, least-squares regression is performed on the whole dataset and on data that is binned into 5 Myr and 10 Myr groups. The Neoproterozoic data show an average of 2-3‰ and an overall trend towards lower values over time. The Ediacaran data show an average of around 0‰ and no clear trend. Thus we argue that there is no step-change towards more positive values over either the Neoproterozoic Era or the Ediacaran period. This absence of a stepincrease is problematic for many proposed mechanisms that infer an increased burial rate of organic carbon relative to carbonates over these intervals (see main text).

COPSE model alterations.
For this work we use the COPSE (Carbon Oxygen Phosphorus Sulphur Evolution) Reloaded model 43 , which is an Earth system box model that computes changes in the C, O, P, S and N cycles over geological timescales to reconstruct long-term climate and sediment geochemistry. It was designed for the Phanerozoic Eon but has been extended previously in simplified forms to test hypotheses for Proterozoic oxygen controls 12,29 , reasoning that the same key processes control the major geochemical cycles. Organic carbon burial in the model is dependent on the concentrations of the limiting nutrients nitrogen and phosphorus, with P being the ultimate limiting nutrient. P input is controlled by weathering fluxes, which are dependent on global temperature and also supply alkalinity and cations, resulting in carbonate precipitation. COPSE is a fullydynamic predictive model in which all fluxes are controlled by the internal processes, rather than being prescribed. Only external forcings (e.g., degassing rate, uplift rate) are prescribed and the resulting predicted changes in global biogeochemical cycling can be tested against isotopic records of carbon, sulphur and strontium 43 .
Tectonic changes and biosphere evolution are imposed as external forcings. See Lenton et al. 43 for the Phanerozoic COPSE model runs and comparisons to δ 13 C, δ 34 S and 87 Sr/ 86 Sr data. In this work, we run the model for the Ediacaran period to test the global response to an increase in CO 2 degassing rates. The degassing rate (D) forcing is prescribed following a quantitative reconstruction for the late Neoproterozoic 37 , and the uplift/erosion (U) forcing is assumed to linearly increase based upon the model fit to the strontium isotope system. Land area forcings for weathering of carbonates, granites and basalts are held at their present-day value given the lack of data, and terrestrial biosphere forcings are turned off (as in the Cambrian part of the Phanerozoic run). We run the model 10,000 times using the Monte-Carlo approach of Royer et al. 51 . Full details of the Monte-Carlo parameter changes and procedure are in the main text. Aside from the standard run denoted above, we also run the model with constant uplift, the results of which are shown in the Supplementary Information (Supplementary Note 1, Supplementary Figs. 2-4).
The model equations used in this version of the model are documented below. Whilst a number of external forcings are set for this work, all model processes remain the same as the published model 43 , with the exception of the addition of direct mantle input of carbon and reducing power (modelled as H 2 gas). We do not include equations relating to the impact of vegetation on the long-term Carbon cycle (i.e., plant-assisted weathering) as these processes were not occurring in the Precambrian and the forcing set used here reduces them to zero. Full model equations including those omitted are documented in Lenton et al. 43 , and we direct the reader to this publication for more detail on the model itself. RO 2 and RCO 2 denote concentrations of oxygen and carbon dioxide relative to presentday. Subscript zeros represent the present-day size of fluxes and reservoirs. Model runs were performed using the MATLAB ODE suite variable timestep solvers 68 in parallel on 6 cores. The COPSE code has been made freely available at https://github.com/sjdaines/COPSE/releases List of fluxes Phosphorus weathering: P delivery to land surface: Land organic carbon burial: P delivery to oceans: Marine new production: newp ¼ r C:P Á min 30:9ðN=N 0 Þ=r N:P ; 2:2ðP=P 0 Þ ð Þ ð 5Þ Marine organic carbon burial: Optional oxygen dependence for mocb: Marine organic phosphorus burial: Calcium-bound phosphorus burial: Iron-sorbed phosphorus burial: Nitrogen fixation: nfix ¼ k 3 P À N=r N:P P 0 À N 0 =r N:P 2 for N r N:P <P; else 0 ð11Þ Marine organic nitrogen burial: Denitrification: Granite weathering: Basalt weathering: Silicate weathering: Carbonate weathering: Oxidative weathering Marine carbonate carbon burial: Seafloor weathering: Pyrite sulphur weathering: Gypsum sulphur weathering: Pyrite sulphur burial: Gypsum sulphur burial Organic carbon degassing: Carbonate carbon degassing: Pyrite sulphur degassing: Gypsum sulphur degassing: Other calculations Atmospheric CO 2 : COPSE Reloaded uses the global temperature function from Berner and Kothavala 69 rather than that of Caldeira and Kasting 70 . Global Temperature: where k c = 4.328°C and k l = 7.4°C Ocean anoxic fraction: Reduced Gas Flux: A flux of reduced gas, modelled as H 2 , is added to the model to explore the possibility that increased subduction and degassing rates may act to lower O 2 levels by delivering more reductant from the mantle. The rate of input is defined by an assumed present-day rate (k rgf ) and scales with the relative degassing rate. To maintain balance of the oxygen cycle at present day, the present day oxidative weathering flux is reduced by k rgf . This modification, in turn, requires an additional source of carbon to the surface system to maintain balance, which is represented by mid-ocean ridge degassing of CO 2 with magnitude equal to k rgf . We also assume that the large sedimentary reservoirs in the model (C, G, PYR, GYP) do not change over time. This avoids further model extension to link these reservoirs to the mantle 71 and is justified because the changes in the relative sizes of these reservoirs over the timescales we wish to consider are negligible (an average of 2.5% change over 100 Myrs in the model of Hayes and Waldbauer 71 ).
Strontium isotope system. This study follows Lenton et al. 43 whereby the strontium cycle and its isotopes are implemented following Francois and Walker 72 and Vollstaedt et al. 73 with some improvements to the formulation described in Mills et al. 74 .
Strontium fluxes Mantle Sr Input: Basalt weathering input: Granite weathering input: Inputs from carbonate sediments Burial in carbonate sediments: Removal in seafloor weathering: The relative proportions of the burial and seafloor weathering removal fluxes of strontium are assumed to follow the same proportions as the corresponding fluxes in the carbon system, with the total flux dictated by assuming present-day steady state for oceanic Sr concentration.
Output from metamorphism: Although there is no fractionation of Sr isotopes associated with the input and output fluxes to the ocean, decay of 87 Rb to 87 Sr influences the 87 Sr/ 86 Sr ratio over long timescales (and is responsible for the differing 87 Sr/ 86  The isotopic composition of the ocean and the sediments are calculated by first creating reservoirs consisting of Sr concentrations multiplied by their isotopic ratios, where δSr X denotes the 87 Sr/ 86 Sr ratio of reservoir X: dðOSr Á δSr ocean Þ dt ¼ Sr granw Á δSr granite þ Sr basw Á δSr basalt þSr sedw Á δSr sediment þ Sr mantle Á δSr mantle À Sr sedb Á δSr ocean ÀSr sfw Á δSr ocean ð44Þ dðSSrÁδSr sediment Þ dt ¼ Sr sedb Á δSr ocean À Sr sedw Á δSr sediment À Sr metam Á δSr sediment ð45Þ The ocean 87 Sr/ 86 Sr ratio is calculated by dividing the new reservoir by the known concentration: The carbonate sediment 87 Sr/ 86 Sr ratio includes an additional term to account for rubidium decay within the sedimentary reservoir: Where here Δt is time elapsed since the start of the model run. The rubidiumstrontium ratio of sediments is calculated to achieve the average crustal 87 Sr/ 86 Sr of 0.73 75 .

Data availability
The model data that support the findings of this study are available from the corresponding author upon reasonable request.