Spatial continuous integration of Phanerozoic global biogeochemistry and climate

It is important to understand how Earth's surface conditions have changed over geological timescales and what has driven these changes. Much of this understanding comes from combining geological and geochemical data with global biogeochemical models


Why model Phanerozoic biogeochemistry?
By developing and running a long-term global biogeochemical model, two key questions about the operation of the ancient Earth can be addressed.Firstly, a model can use known metrics to make predictions for metrics which are unknown (a 'reconstruction' approach), and secondly a model can be used to interrogate the drivers of changes that are observed in global biogeochemistry (a 'hypothesis testing' approach).Through these approaches more can be learned about how the Earth's climate and surface chemistry has changed over time and what was responsible for these changes (e.g.Berner et al., 1983;Berner, 1991;Bergman et al., 2004;Arvidson et al., 2013;Goddéris et al., 2014;Lenton et al., 2018).These discoveries have important wider implications for untangling the relationship between surface conditions and biological evolution (e.g.Clapham and Karr, 2012), understanding how our planet will respond to present day climate change (e.g.Archer et al., 2009), and assessing the likelihood of habitable conditions being maintained on other planets (e.g.Cockell et al., 2016).
The field of global biogeochemical modelling is still relatively young, and much work remains to be done.When assessing the Phanerozoic as a whole, there is currently no global biogeochemical model that can satisfactorily replicate the major shifts in climate and surface chemistry (e.g.CO 2 levels) that have occurred, even at the broad ~100 Myr scale (Mills et al., 2019;Goddéris and Donnadieu, 2017).In this paper we discuss the reasons for this, and then draw on previous approaches to construct a new type of framework which has the potential to allow for a better understanding of Phanerozoic biogeochemistry and climate.

Current models and problems in timescales and computation
Early biogeochemical models computed fluxes in the global carbonate-silicate cycle to explore control of global surface temperature and atmospheric CO 2 concentration over the last 100 Myrs (Berner et al., 1983) and later the whole Phanerozoic (Berner, 1991).These computations were carried out in a 'nondimensional' or 'box model' framework, meaning that each process (e.g.degassing, weathering) and reservoir (e.g. the entire ocean) was represented by a single variable at the global scale.This approach is simple, but this simplicity meant that researchers could begin to test hypotheses about what controlled the CO 2 concentration of Earth's atmosphere over many millions of years without requiring powerful computing resources or sophisticated mathematics.
The current state of the art in Phanerozoic biogeochemical modelling is still largely dominated by box models.Specifically, by the GEOCARB (Berner, 2006;Royer et al., 2014;Krause et al., 2018) and COPSE (Bergman et al., 2004;Lenton et al., 2018;Mills et al., 2019) box models, which have been widely used to predict Phanerozoic climate and surface chemistry evolution, and to assess what drives major changes.Of course, the current generation of box models include more complexity than their predecessors: in recent versions both GEOCARB and COPSE are now solved using implicit numerical techniques designed for 'stiff' systems (where timescales vary substantially between processes), and both include a larger range of species and reactions.However, both models remain nondimensional.
An important limitation of nondimensional models is that they must represent surface properties such as temperature, precipitation and erosion as a single value at the global scale.Understandably therefore, the key process of continental weathering has been shown to be poorly predicted in box models when compared to spatially-explicit modelling (Donnadieu et al., 2006;Taylor et al., 2012;Goddéris et al., 2014).This is a major limitation because continental weathering is thought to be responsible for regulating atmospheric CO 2 , O 2 and marine sulfate concentrations over long timescales (Walker et al., 1981;Lenton and Watson, 2004;Bolton et al., 2006;Wortmann and Paytan, 2012), and also for driving complex effects on the carbon isotope record (Daines et al., 2017;Shields and Mills, 2017).
The reason that box models are still widely employed, despite this obvious weakness, is that spatially-resolved 3D Earth system models, even at low resolution (e.g.GENIE: Edwards and Marsh, 2005;Ridgwell et al., 2007) cannot be run over Phanerozoic timescales due to their complexity.For example, the 3D ocean model in the simplest GENIE configuration (e.g.Lord et al., 2018) uses a resolution of 36 × 36 grid cells, and 8 depth levels, or ~10,000 individual boxes, and can be run for 1 million model years (Lord et al., 2018;Colbourn et al., 2013), but would need to be speeded up by a further ~500-fold to run for Phanerozoic time.An accurate depiction of continental weathering processes also necessitates a more complex treatment of the atmosphere than available in GENIE, where this too is split into a 3D grid of cells, adding further complexity and further reducing integration times.When considering the fully-coupled atmosphere-ocean GCMs that are used for individual deep-time paleoclimate studies (e.g.EoMIP: Lunt et al., 2012), a standard simulation length is on the order of several thousand model years.Dynamically reconstructing the Phanerozoic carbon cycle with a General Circulation Model in any useful timeframe is currently impossible, and will be until a step-change in computing resources occurs.
Nevertheless, a useful compromise between the long run times of spatial models and the oversimplification of box models has been reached in the GEOCLIM model (Donnadieu et al., 2004(Donnadieu et al., , 2006;;Goddéris et al., 2014Goddéris et al., , 2017)).In GEOCLIM, a datastructure (or 'lookup table ') of previously-completed General Circulation Model simulations is used to inform surface process modelling in a biogeochemical box model.In the GEOCLIM approach, the current atmospheric CO 2 concentration of the box model is used to interpolate between the GCM runs to arrive at an approximate grid of continental runoff and surface temperature, which allows for a 2D approach to terrestrial weathering, returning bulk fluxes to the box model.GEOCLIM is a climate-informed biogeochemical model which can be integrated over multimillion-year timescales, and has led to several important breakthroughs in understanding the relationship between continental configuration, mountain uplift and long-term climate change (e.g.Donnadieu et al., 2004;Goddéris et al., 2017).But GEOCLIM cannot be run under a continually changing continental configuration, meaning that it is limited to producing steady state 'snapshots' of the Earth system, or relatively short (<10Myr) transient events (Donnadieu et al., 2011;Park et al., 2020).This means that COPSE and GEOCARB (or similar nondimensional systems, e.g.Arvidson et al., 2013) have remained the only modelling approaches that can produce predictions of atmospheric O 2 , marine sulfate concentration, marine nutrient inventories and geological isotope tracers over very long timescales, and they remain at the forefront of investigations into Phanerozoic and late Precambrian climate and biogeochemistry (e.g.Krause et al., 2018;Tostevin and Mills, 2020).

Overview
In this paper we link the spatial climate modelling procedure developed for GEOCLIM to the long-term biogeochemical processes in COPSE to produce a continuous-in-time global biogeochemical model that calculates Earth surface processes in 2D using realistic approximations of climate.We term this the Spatial Continuous Integration (SCION) model.A schematic representation of the model is shown in Fig. 1, where green arrows show dataflows and processing and black arrows show biogeochemical fluxes.As in GEOCLIM, a set of climate model (FOAM) outputs are arranged in a datastructure (Fig. 1A), allowing the model to access global 2D fields for surface temperature, topographic height and runoff at a number of pre-determined Phanerozoic timepoints and for a wide range of CO 2 concentrations.GEOCLIM is run for a chosen timepoint only, and the SCION approach differs from this by looking at two timepoints simultaneously and calculating all surface fields and processes for both of these at each model timestep.A weighted average of the resulting global fluxes is then used to compute surface processes at any time between the two chosen 'keyframe' timepoints, assuming a linear relationship in time between the two.This process logically extends to any timepoint during the Phanerozoic, allowing the model to produce continuous predictions for 540 million years.
Due to the linear interpolation required, and the relatively coarse resolution of the available FOAM timepoints (average ~20 Myrs between runs), this approach may miss certain features.For example, the movement of India northwards through the tropics during the late Mesozoic and early Cenozoic (Peirce, 1978) is represented by GCM runs where India is south of the equator (70 Ma) and then north of the equator (52 Ma), and the fluxes during the intervening period are calculated based only on these end-members, without considering the potential spike in precipitation and temperature around the equator itself.Running the model with a more finely spaced set of GCM simulations could reduce this problem in the future.

Interpolation routine
SCION uses the same FOAM GCM datastructure developed for GEOCLIM (Goddéris et al., 2014), and the biogeochemical box model is adapted from the latest version of COPSE (Tostevin and Mills, 2020).Thus SCION continues to use a single box to represent the atmosphere and ocean, and boxes to represent the sedimentary inventories of the    different chemical species (see Fig. 1).The model species are shown in Table 1, where each inventory evolves during the model run subject to the inputs and outputs described in Table 2.At each timestep the model requires information from the GCM datastructure to calculate surface processes.The model uses the current CO 2 concentration and geological age to generate a 'Gridstate' for the two keyframes that bracket the current model timepoint ('gridpast' and 'gridfuture'; see Fig. 1).These Gridstates consist of topographic height, runoff and air temperature, an example of which is shown in Fig. 2. At each timestep, 2D surface calculations are run for both gridpast and gridfuture, and a weighted average is taken for the final bulk flux.For example, for a generic GCM field F, we calculate: Where D U and D L are the distances in log space between the model current CO 2 value and the keyframe CO 2 values for the upper (U) and lower (L) GCM runs.Log space is used because the effect of CO 2 variation on global surface temperature and runoff is approximately logarithmic, and both are calculated based on the CO 2 values.We then run all surface 2D calculations using both gridpast and gridfuture.Finally, bulk fluxes to be passed to the box model are calculated by summing each of the gridpast and gridfuture fields over the land surface, and interpolating over time.
Where T F and T P are the fractional distances in time to the future and past GCM keyframes, and A grid is the gridbox area, which is dependent on latitude.

Continental weathering fluxes
Although this work uses the climate simulations from GEOCLIM, the weathering calculations do not follow other GEOCLIM publications.Basin-scale silicate weathering rates are calculated using the following parametric relationship, from West (2012), which combines dependencies on local runoff, temperature and erosion rate.
Where ω sil is the local silicate cation denudation flux and the kinetic term is defined by: Where dependencies on local runoff, temperature and erosion rates are calculated as: Erosion rate (ε) is calculated from topographic slope (s) and local runoff (Q) using the approach of Maffre et al. (2018).
The basaltic and granitic fraction of silicate weathering is calculated from the total silicate weathering rate based on the relative exposed areas of these lithologies (e.g.Berner, 2006b;Mills et al., 2014).This assumes a homogenous distribution.
Carbonate weathering is assumed to scale with runoff (e.g.Berner, 1994), where an additional parameter, k scale , is added to separate the present day rate from the spatial scaling effect.
Oxidative weathering ( f oxidw ), pyrite weathering ( f pyrw ) and gypsum weathering ( f gypw ) are all also assumed to be dependent on local runoff, given the general requirement for moisture to facilitate aqueous chemical reactions.But these processes also have other dependencies at the global scale, as in the COPSE model: Phosphorus weathering sums contributions from silicates, carbonates and organics, as in COPSE: Thus, all weathering fluxes are calculated in 2D from the spatial climate reconstructions.While the process is relatively detailed for silicate weathering, the other weathering fluxes currently only have a simple relationship to the climate, which may be improved with future work.In addition to the weathering fluxes, we also use information from the Digital Elevation Model (DEM) to assist in calculating the rate of gypsum deposition in the model.Records of evaporite gypsum burial show sporadic and massive events that usually cluster around the timing of basin restriction (Warren, 2010;Wortmann and Paytan, 2012), and as a rough approximation to this process we add a multiplier for the length of coastline (L C ) to the gypsum burial calculation: This acts to increase the rate of gypsum deposition when basins are closing, although it is a crude representation.

Nondimensional fluxes
Table 3 shows the nondimensional fluxes in the model, which generally relate to the degassing and burial of hydrosphere species, but also

Flux name Equation
Carbonate C degassing: Marine carbonate burial: Fe-phosphate burial: Ca-phosphate burial: Nitrogen fixation: P flux to land:

Table 4
Other model processes.

Process Equation
Carbon atmospheric fraction Seafloor weathering T effect: Temperature effect on vegetation: Oxygen effect on vegetation: Overall limitation of terrestrial NPP: Fire ignition probability scaling: Fire effect on terrestrial biomass: weathering effect: include some internal fluxes which are the distribution of weathered phosphorus between land and ocean, and the ocean-atmosphere fluxes of nitrogen species.All of these fluxes are taken directly from COPSE (Tostevin and Mills, 2020) and the reader is directed to the most recent COPSE review paper (Lenton et al., 2018), or original COPSE publication (Bergman et al., 2004) for further information on how they were constructed.These fluxes generally take the form of a present day rate multiplied by a series of scalings, which include the size of the parent reservoir, forcing factors, and non-flux calculations such as the degree of marine anoxia or the global average surface temperature (T surf ).
One alteration from the most recent version of COPSE is that sulfur degassing fluxes are not considered here, a reversion back to the original COPSE formulation which only considered sulfur weathering and burial (Bergman et al., 2004).As discussed in Lenton et al. (2018) this change makes little difference to the model biogeochemistry predictions, but the addition of sulfur degassing fluxes were not compatible with the COPSE calcium cycle.Although we do not consider a Ca cycle here, future work that does so is planned.

Non-flux calculations and fixed parameters
Table 4 shows the rest of the processes and calculations in the model, which also follow COPSE, and Table 5 shows all fixed parameters.The majority of these additional calculations relate to calculating the mass of the terrestrial biosphere and resulting biotic amplification of weathering (f biota : e.g.Schwartzman and Volk, 1989), which is the used to scale the silicate and carbonate weathering rates Where i denotes either basalt, granite or carbonate weathering.For this paper we continue to calculate the mass of the terrestrial biosphere and biotic weathering effects using globalised values for temperature, CO 2 concentration and O 2 concentration, as in COPSE.This is a major limitation of the current work, as the spatial climate and topography scheme would allow for a much more detailed representation of the terrestrial biosphere, although it represents a substantial amount of development work.A global vegetation model is included in the FOAM climate runs, but it embodies fixed assumptions about oxygen concentration, wildfires and nutrient uptake which are not compatible with the current approach.

Nondimensional forcings
Fig. 3 shows the dimensionless input forcings required to run the model.As we are using a spatial climate scheme, this forcing set is considerably smaller than other continuous Phanerozoic models like GEOCARB and COPSE.This is because it does not need to include dimensionless reconstructions of global erosion rates (e.g.Ronov, 1976;Hay et al., 2006), paleogeographic effects or land surface temperatures (Otto-Bliesner, 1995), which are difficult to apply as a single value at the global scale.Aside from the basalt and granite areas mentioned above, the SCION model requires a degassing rate 'D', carbonate burial depth 'B' (which scales the degassing rate for carbonates, e.g.Volk, 1989), and a set of biological forcings that represent the evolution of a terrestrial biosphere (plant evolution, 'E'), the effect of plant colonisation on global chemical weathering rates (biotic weathering, 'W') and the effect of sediment bioturbation on marine carbon and sulfur burial ('SB'), all of which follow the latest COPSE model (Tostevin and Mills, 2020).Two COPSE forcings which could be applied to SCION are omitted due to their uncertainty.These are the assumptions of a high C:P ratio and selective phosphorus weathering enhancement associated with early terrestrial nonvascular plants (Lenton et al., 2016).The intention with this model is to include only the most essential and most broadly supported model forcings, although these of course remain somewhat uncertain in their specific trajectories.The error window on the degassing rate represents the boundaries of two different estimatesthe lower from the length of subduction zones and rifts (Mills et al., 2017(Mills et al., , 2019;;Brune et al., 2017) and the upper from kinematic plate modelling of the total material destruction rate (Seton et al., 2012;Domeier and Torsvik, 2017).The error estimate on the basalt and granite areas is assumed to be ±20% (Mills et al., 2014).

Forwards model ensemble predictions compared to geochemical data
The SCION model is solved in MATLAB using the variable-step variable-order implicit routine for stiff systems (ODE15s: Shampine and Reichelt, 1997).Each plot shown is constructed from 1000 model runs in which several key uncertainties are tested: degassing rates and relative areas of basalts and granites are varied between the boundaries shown in Fig. 3; the weathering enhancement due to plant evolution is varied between no enhancement and a 7-fold enhancement (Lenton, 2001;Bergman et al., 2004;Quirk et al., 2015); the isotope fractionation during microbial sulfate reduction is varied between 20 and 40‰  (Canfield and Farquhar, 2009), and fractionations associated with photosynthesis are varied between 20 and 30‰ (land) and 25-35‰ (marine) (Bergman et al., 2004).
The evolution of the key model surface reservoirs CO 2 , O 2 and SO 4 are shown in in Fig. 4, with Fig. 5 showing the global average surface temperature and the model 'ice line' representing the rough positioning of the continental ice sheets.The ice line is estimated from the latitude at which the model continental surface temperature reaches that of the present day ice line latitude (~−10C e.g.Caldeira and Kasting, 1992).Each of these plotted metrics are readily comparable to independent geological and geochemical evidence, and can be used to validate the model and explore what processes may drive the variations in  Foster et al., 2017;Witkowski et al., 2018).C. Atmospheric O 2 concentration in the model (blue) plotted against the interpretation of the charcoal record (yellow; Glasspool and Scott, 2010).D. Marine sulfate concentration in the model (blue) plotted against fluid inclusion data (vertical bars; Horita et al., 2002;Brennan et al., 2004;Lowenstein et al., 2005), and based on sulfur isotope fractionation factors (yellow shading; Algeo et al., 2015).Phanerozoic surface chemistry and climate.Further direct comparisons can be made between the modelled stable isotope ratios δ 13 C, δ 34 S and 87 Sr/ 86 Sr, and their geochemical records, which are shown in Fig. 6.

Atmospheric CO 2
Panel 4B shows the model ensemble (and mean) atmospheric CO 2 concentration in parts per million (ppm), compared to recent proxy compilations (Foster et al., 2017;Witkowski et al., 2018).From Devonian to present day, the model predictions for variations in atmospheric CO 2 concentration follow the broad pattern shown in the data compilation reasonably well: CO 2 concentration declines through the Devonian and Carboniferous, reaches a nadir around the Carboniferous-Permian boundary, increases over the Permian and Triassic, decreases in the Jurassic, peaks again in the Cretaceous before declining over the Cenozoic.These dynamics are driven by a combination of the changing CO 2 degassing rate and the degree of amplification of silicate weathering at the global scale.From Jurassic to present the changes in CO 2 concentration quite clearly track the rate of CO 2 degassing, whereas for the Paleozoic and Triassic, a weathering driver is more apparent (as noted by Goddéris and Donnadieu, 2017).The CO 2 minimum at the Carboniferous-Permian is associated with the uplift of the Hercynian mountains in the tropics (Fig. 4A), which amplifies silicate weathering rates (Goddéris et al., 2017), whereas the shift to very high CO 2 in the Triassic is associated with the denudation of these mountains combined with the aridity of Pangea, which suppresses chemical weathering.
In general, modelled CO 2 is higher than the proxy data during the Paleozoic, and during the Cretaceous.However, model global average surface temperature is not in excess of proxy estimates during these times (Fig. 5; Scotese et al., 2021).This may be due to low climate sensitivity in the FOAM climate model, which means that an unrealistically large amount of CO 2 is required to raise the surface temperature to the level at which the carbon cycle is balanced.More detailed climate models (especially the recent CMIP6 models) tend to have higher climate sensitivities than the FOAM model used here (Tierney et el. 2020).
Another mismatch between the model and the data compilation is the failure of the model to capture low CO 2 as interpreted from phytane measurements from the Ordovician-Silurian.Previous hypotheses for this fall in CO 2 note the possibility for rapid weathering of volcanic arcs and/or decreases in degassing rate linked to arccontinent collision (Young et al., 2009;Macdonald et al., 2019), or large transient weathering rate amplifications associated with early  Cather et al., 2009).C. Global average surface temperature in the model (blue) plotted against global average surface temperature derived from paleo-Köppen belts and oxygen isotopes (Scotese et al., 2021).land plants (Lenton et al., 2012).The current SCION model does not provide an answer to how important each of these mechanisms was, but it is the ideal framework in which to test these ideas in the future, given that it can produce records of stable isotope ratios as well as weathering dynamics.Such tests are contingent on a clearer approximation of land plant dynamics and volcanic weathering in future versions of the model.

Atmospheric O 2
Atmospheric O 2 is relatively stable when compared to atmospheric CO 2 .This is because the model features strong feedbacks on both O 2 production and O 2 sinks, which are inherited from the COPSE model.Production-based feedbacks are the effect of wildfire (Kump, 1988) and photorespiration on limiting the terrestrial vegetation (Bergman Fig. 6.Model stable isotope records.A. Summary of topography change.Here 'ice sheets' are drawn where T < − 10 °C.B. Model carbonate δ 13 C (blue) compared to the geological record (Saltzman and Thomas, 2012).C. Model seawater δ 34 S (blue) compared to the geological record in evaporites, barites and CAS (Crockford et al., 2019) , 2004), and also the ability of marine oxygenation to increase the net burial rate of the nutrient phosphate and reduce the preservation potential of organic carbon (VanCappellen andIngall, 1994, 1996;Lenton and Watson, 2000;van de Velde et al., 2018).The pattern of atmospheric O 2 variation generally follows that of the COPSE model.A Cambrian drop in O 2 is driven by the initiation of marine bioturbation and associated decrease in organic carbon preservation (van de Velde et al., 2018), a rise through the Silurian and Devonian is driven by the evolution of first nonvascular and then vascular land plants, which increases global organic carbon burial rates (Bergman et al., 2004;Lenton et al., 2016), and the dynamics from Carboniferous to present tend to follow the changes in the rate of CO 2 input through degassinghere a larger input of carbon results in a larger output of both organic and inorganic carbon (e.g.Williams et al., 2019).The model is compared to the O 2 reconstruction of Glasspool and Scott (2010), which is based on sedimentary charcoal abundance, and generally compares well.But the model does not produce highly variable oxygen levels during the Mesozoic, due to the large number of negative feedbacks present.Models that split the ocean into multiple boxes (e.g.Slomp and Cappellen, 2007;Alcott et al., 2019;Wallmann et al., 2019) are better able to represent marginal settings where marine redox dynamics, and thus the burial rates of reduced species, can fluctuate much more rapidly than in the single-box ocean of COPSE which is currently used in this model.So, it is possible that oxygen concentrations are over-regulated in the current framework.

Marine sulfate
The SCION predictions for marine sulfate concentration are shown in Fig. 4D, and are compared to individual fluid inclusion measurements (Horita et al., 2002;Brennan et al., 2004;Lowenstein et al., 2005) as well as the broad area defined by the changes in sulfur isotope fractionations between buried sulfates and pyrite (Algeo et al., 2015).The model generally predicts a high concentration of marine sulfate for the Paleozoic, sitting above the combined proxy window, and clearly above some of the fluid inclusion data.The fit for the Permian to Jurassic is better, but later Mesozoic and Cenozoic model results are again higher than the fluid inclusion measurements.As with O 2 , marine sulfate is strongly regulated in the SCION model, as both the burial rates of pyrite and gypsum are assumed to vary linearly with sulfate concentration, which is inherited from the COPSE model.Therefore it is relatively difficult to force sulfate concentrations far from the present day.Large evaporite weathering and deposition events also remain poorly represented in the model, and these may substantially alter marine sulfate concentrations (Wortmann and Paytan, 2012;Shields et al., 2019;Shields and Mills, 2020).

Glacial record
The Paleoglaciation ice line latitude (e.g.Crowley, 1998) is a key piece of Phanerozoic climate information, and is one of the most robust datasets available on global temperature before the Cenozoic.The model reconstruction of this record (Fig. 5B) is based only on surface temperatures and is therefore highly uncertain, also during the Paleozoic the model uncertainty window is large (as for CO 2 ) due to the uncertainty around the effects of land plants on the carbon cycle.Overall, the Paleoglaciation record from the model is relatively consistent with data from the Cretaceous to the present day, supporting the idea that CO 2 was the key greenhouse gas over this period.The Paleozoic and earlier Mesozoic are not well reproduced, however.The model generally predicts ice sheet growth from the later Devonian until a Carboniferous-Permian peak, and a retreat during the Permian, which is qualitatively consistent with data, but it also predicts a period of ice sheet growth during the Triassic-Jurassic, which is not consistent with data.
As may be expected from the mismatch in CO 2 , the model ice line does not reproduce the Hirnantian glaciation in the late Ordovician.Here, the reasons for this mismatch are likely the same as for the mismatch in CO 2 as discussed above.For the Triassic-Jurassic, the model predicts a CO 2 concentration within the proxy window, but has a significantly lower surface temperature and much more ice than is evident in the geological record.This is an intriguing problem and may testify to the importance of a non-CO 2 greenhouse gas or other temperature forcing (e.g.albedo), or may simply be a limitation of the FOAM climate model (e.g.low climate sensitivity).

Global average surface temperature
Fig. 5C compares the global average surface temperature in the model to the record derived from a combination of paleo-Köppen climatic belts and oxygen isotopes (Scotese et al., 2021).Overall, the long-term fit is acceptable, but there are clear mismatches in the late Ordovician (lack of cooling spike) and Triassic-Jurassic (too cold).The degree of cooling during the Late Paleozoic Ice Age is also not fully represented in the model.The potential reasons for the lack of Ordovician cooling and lack of Triassic warmth are discussed earlier.The model does not capture any of the shorter term (~1-5 Myr) temperature variations apparent in the proxy record.These variations are likely the result of processes not included in the model, such as LIP emplacements or bolide impacts (e.g.Scotese et al., 2021).
3.6.Carbonate δ 13 C Fig. 6B compares the SCION δ 13 C ensemble to the isotope record compiled by Saltzman and Thomas (2012).Many shorter-term (e.g.5-50 Ma) trends are not replicated by the model, but an overall trend of a rise in the mid-Paleozoic, and generally stable values through the Mesozoic, is reproduced.Although the decline in δ 13 C during the Cenozoic is not evident.This Paleozoic rise is due to land plant evolution and additional carbon burial, and the stability is partly due to the sensitivity of organic carbon weathering fluxes to atmospheric O 2 , which dampen any carbon isotope change associated with changes in organic carbon burial (Daines et al., 2017).It is also probable that the single-box ocean in the model, and lack of consideration of marginal environments, averages out any broader carbon isotope excursions which may otherwise occur.A key mismatch in the long-term fit is that the model fails to reproduce very high carbon isotope values during the middle to late Paleozoic.Previously the COPSE model did reproduce these values because it assumes a very high carbon-to-phosphorus ratio of early vegetation, which continues into the Carboniferous coal deposition period (Lenton et al., 2016(Lenton et al., , 2018)).This also results in estimation of higher O 2 concentrations.These high C:P ratios are not incorporated in this model as they would need to be part of a more detailed spatial terrestrial biosphere component, and the intention here is to concentrate on the effects of improving the representation of climate.It seems likely that a more complete treatment of the terrestrial biosphere may improve the model fit to the δ 13 C record.

Sulfate δ 34 S
As shown in Fig. 6C, the sulfate δ 34 S record as approximately the inverse of the carbonate δ 13 C record.Here the model is plotted against the compilation of Crockford et al. (2019).As with carbon isotopes, longterm changes are replicated reasonably well but short-term changes are not.Once again, an important factor here is the lack of spatial representation of the ocean, in addition to the heterogeneous nature of evaporite weathering (e.g.Wortmann and Paytan, 2012).In the model, long-term sulfur isotope changes are driven by the changing rate of pyrite burial, which responds to the prevalence of marine anoxia in the model, which is elevated in the Paleozoic, as well as the input of isotopically-heavy sulfate from continental weathering, which is restricted during the time of Pangaea.

Strontium 87 Sr/ 86 Sr
Seawater strontium isotope ratios in the model (Fig. 6D) are poorly predicted when compared to the data of McArthur et al. (2012), which is especially apparent given the high fidelity of the strontium geological record due to strontium's long marine residence time.Strontium isotopes are primarily controlled by the isotopic composition of weathered lithology (Brass, 1976), and previous attempts to predict the Phanerozoic strontium isotope record in full have been hampered by the lack of detailed surface lithological data available (e.g.Mills et al., 2014).The broad increase in 87 Sr/ 86 Sr values observed from Cretaceous to present is apparent in the model, and is driven by a sharp reduction in the area of exposed basalts assumed in the model, as well as a decline in hydrothermal inputs (Fig. 3).A clearer picture of the variations in 87 Sr/ 86 Sr over the Phanerozoic requires the mapping of the positions of different lithologies over time.While this is difficult, some progress can be made at least with respect to the positioning of volcanic terranes and LIPs (e.g.Lefebvre et al., 2013), and the SCION model will be capable of integrating these ideas in the future.

Comparisons to the GEOCLIM and COPSE models for atmospheric CO 2
The SCION approach couples the spatial climate model dataset used in GEOCLIM to the nondimensional biogeochemical fluxes of the COPSE model over a continuous timeframe.It may be expected that the key model predictions for atmospheric CO 2 levels fall somewhere between the predictions of these two models, or may combine aspects of both systems.Atmospheric CO 2 predictions from SCION (ensemble mean), GEOCLIM and COPSE are plotted in Fig. 7 against the combined proxy record.GEOCLIM has been run under a range of different degassing rate estimates (Goddéris and Donnadieu, 2017) and we choose the version that is run under the relative degassing rate derived by Gaffin (1987) from sea-level inversion, because this curve is the most similar to the degassing rate used in this study.The degassing rate used in COPSE (Mills et al., 2019) is the same as the lower bound of the degassing rate used here in SCION.Overall, all models assume a degassing rate that gradually decreases over the Paleozoic, peaks in the Cretaceous then declines towards the present day, and all model CO 2 predictions show a gradual decline over much of the Paleozoic, and a Cretaceous peak.Another important note is that the GEOCLIM compilation shown here does not include the most recent work on the model to incorporate the effects of erosion and soil shielding, which alters the model results substantially but has not been carried out for all time periods (e.g.Goddéris et al., 2017).
When comparing the SCION model to COPSE, there are some clear differences in the CO 2 trajectories.Between the Devonian and the Jurassic, the models are almost in anti-phase.COPSE predicts a decline towards the lowest CO 2 values at the Devonian-Carboniferous boundary, followed by a rise over the Carboniferous and Permian, whereas SCION has the nadir in CO 2 around 50 million years later at the Carboniferous-Permian boundary.Similarly, COPSE predicts generally declining CO 2 between the mid-Permian and the mid-Triassic whereas SCION predicts a rise over the same period.When compared to the proxy dataset it can be seen that the SCION model has a trajectory that is in better agreement with the proxies, although the generally lower values for CO 2 in COPSE are sometimes closer to the proxy record.
This comparison demonstrates the degree to which CO 2 levels are likely to depend on the amplification of silicate weathering, and how this is very difficult to replicate in a nondimensional model.The significant CO 2 drop at the Carboniferous-Permian, corresponding to the Late Paleozoic Ice Age, is driven in SCION by the uplift of the Hercynian mountains in the topics (e.g.Goddéris et al., 2017), which combines high rainfall with high rates of erosion and reasonable surface temperatures to drive high rates of silicate weathering.In COPSE (and in GEOCARB and other nondimensional systems), rainfall and erosion rates must be prescribed at the global scale and this combination of conditions in the tropics cannot be representedthe global rate of erosion at the Carboniferous-Permian is not particularly high for example (Hay et al., 2006).From Jurassic to present, SCION and COPSE are in closer agreement over CO 2 trajectories, but SCION predicts much higher CO 2 during the Cretaceous and COPSE has low CO 2 over the Paleogene, neither of which are well supported by the proxy compilation.When comparing SCION to GEOCLIM, the CO 2 trajectories are similar from Cambrian to Cretaceous, but the SCION predictions are generally lower and therefore closer to the proxies, and feature a more prominent dip at the Carboniferous-Permian.Both models rely on the same climate scheme, and similar approximations of spatial silicate weathering, although SCION includes a representation of the effect of erosion on silicate weathering whereas the plotted version of GEOCLIM does not.The lower CO 2 levels in SCION also likely stem from the incorporation of functions from the COPSE model for global biogeochemistry, supported by the COPSE predictions also being generally lower than GEOCLIM for this period.SCION has inherited a number of strong negative feedbacks on CO 2 concentration from COPSE, which are not treated the same way in GEOCLIM.Firstly, SCION has a flux of seafloor weathering (hydrothermal carbonatization) which results in a net flux of carbon out of the ocean, and is temperature-sensitive (Sleep and Zahnle, 2001;Coogan and Gillis, 2013), whereas GEOCLIM does not.Secondly, SCION assumes that burial of organic carbon is temperature and CO 2 sensitive through both the direct effects on the terrestrial biosphere, and through changes to the delivery of the marine limiting nutrient phosphate.GEOCLIM does not replicate these functions exactly, although its marine biosphere will also bury more carbon when temperatures rise due to the spread of anoxia.Through the Cretaceous and Cenozoic, SCION predicts higher CO 2 concentrations than GEOCLIM because the assumed degassing rate is higher, and because it incorporates a dependency of weathering rates on local erosion, which has not yet been added to GEOCLIM for this period.

Controls on Phanerozoic biogeochemistry
This modelling exercise is the first to link a spatial climate scheme to global biogeochemical fluxes over the Phanerozoic in a continuous way.The result is a somewhat reasonable reconstruction of CO 2 , surface temperature and O 2 levels from Devonian to present day and a long-term (although dampened) agreement with isotope tracers δ 13 C and δ 34 S. Model results for marine sulfate concentration, the expansion of ice sheets and sedimentary 87 Sr/ 86 Sr ratios remain unsatisfactory, as does the general over-predictions of CO 2 levels.Our model supports the hypothesis that Phanerozoic CO 2 concentration has been largely controlled by changes to silicate weathering and degassing rates, and more specifically that the Late Paleozoic Ice Age was driven by mountain uplift, whereas cooling from the Cretaceous to present was driven largely by decreasing degassing rates, in combination with Himalayan uplift.
Phanerozoic oxygen levels are decoupled from CO 2 concentration.Although the organic carbon cycle produces O 2 from CO 2 (and vice versa), enough independent processes operate on the carbon cycle alone to remove any simple relationship between the two over long timescales.Sulfate concentrations are controlled by weathering and deposition events as well as changes to marine carbon cycling and redox and are thus difficult to reproduce in a model.Because the model predicts sulfate δ 34 S reasonably well, the times of mismatch between the SCION model and the δ 34 S record may be due to large evaporite deposition events, which do not alter seawater δ 34 S. As noted earlier in the paper, the strontium isotope record is very sensitive to the age and type of weathered silicates, which is currently not well represented in the model.Major features are missing form the model reconstructions for the Cambrian-Silurian.Here, predicted CO 2 concentration and ice sheet advance does not capture the late Ordovician glaciation, and sulfate levels are generally higher than proxies suggest.

Future directions
A relatively small number of model improvements should help better understand some of the mismatches between the SCION predictions and the geological record.Firstly, a spatial treatment of mafic versus felsic silicate weathering is essential.This is easily implemented in the model by mapping volcanic terranes onto the digital elevation model, allowing their contributions to weathering to be more accurately assessed (e.g.Lefebvre et al., 2013;Goddéris et al., 2017).It is expected that this would lead to a better fit to the strontium isotope record, and will also change the model CO 2 and temperature predictions because mafic lithologies are more reactive.Specifically, arc weathering has been proposed as a driver of Ordovician glaciation (Young et al., 2009;Macdonald et al., 2019) and Neogene cooling (Park et al., 2020) and these hypotheses can be tested in the model.
Secondly, a number of model inconsistencies are likely tied to the poor representation of the terrestrial biosphere.Currently the model has a globally-averaged terrestrial biomass which impacts weathering rates and carbon burial, but there is clear potential to add a spatially-resolved biosphere which would likely behave in a different way by integrating local effects of photosynthetically active radiation, water availability and temperature.A more detailed model terrestrial biosphere may alter the predictions for δ 13 C during the late Paleozoic, and may alter CO 2 and O 2 levels throughout the simulations.Land plant evolution has also been proposed as a driver of late Ordovician glaciation (Lenton et al., 2012).In addition to altering the model chemical cycles, a more sophisticated land biosphere and terrestrial surface scheme could aid in predicting properties of paleosols which could be compared to the geological record (e.g.Retallack, 1997Retallack, , 2009)).
Another achievable modification is improved representation of the ocean.SCION currently uses a single-box ocean model inherited from COPSE, and ultimately from the very first long-term biogeochemical models (Walker et al., 1981;Berner, 1991).Recent work has shown that using multiple ocean boxes to represent shelf environments reduces the amount of negative feedback in the system and allows for more rapid changes in global redox and nutrient cycling (e.g.Slomp and Cappellen, 2007;Alcott et al., 2019;Wallmann et al., 2019).It is possible that this may increase the variability of atmospheric O 2 in the model, and will also impact the sulfur cycle through redox changes.
Finally, this framework may later be expanded to consider geochemical tracers that were not easily integrated into current box models.Two examples are the lithium cycle and δ 7 Li, and the marine δ 18 O ratio.Both can be explicitly modelled in SCION because it can approximate ocean surface temperature, and can also approximate the balance between weathering and clay formation.

Fig. 1 .
Fig. 1.SCION Earth System Model schematic.A. Visualization of the General Circulation Model datastructure used to produce the evolving model 'gridstate'.B. The underlying biogeochemical box model.Green arrows show dataflows and processing, black arrows show biogeochemical fluxes.See text for further details.

Fig. 2 .
Fig. 2. Example GCM fields.This example is for 300 Ma and 350 ppm CO 2 .The model uses a lookup table of climate fields, derived from FOAM runs of 21 continental reconstructions at varying CO 2 levels.A. Topography used for climate model run.B. Output air temperature at surface level.C. Output continental runoff.

Fig. 3 .
Fig. 3. Nondimensional forcings.A. Tectonic forcing controlling input rate of CO 2 .B. Land area forcings affecting the strontium isotope composition of silicate weathering fluxes.C. Biological forcings affecting carbon burial and weathering.All are based on the COPSE model.

Fig. 5 .
Fig. 5. Phanerozoic model results for ice cap latitude and global temperature.A. Summary of topography change.Here 'ice sheets' are drawn where T < − 10 °C.B. Ice line latitude: the extent ice caps in the model (blue) plotted against the geological record (yellow; as compiled byCather et al., 2009).C. Global average surface temperature in the model (blue) plotted against global average surface temperature derived from paleo-Köppen belts and oxygen isotopes(Scotese et al., 2021).

Fig. 7 .
Fig. 7. Comparison of atmospheric CO 2 evolution in different types of biogeochemical model.The central SCION predictions (continuous, spatial surface model) are compared to GEOCLIM (snapshots, spatial surface model), and COPSE (continuous, nondimensional surface model), as well as to the combined proxy CO 2 record (Yellow;Foster et al., 2017;Witkowski et al., 2018).
Benjamin J. W. Mills is Associate Professor at the University of Leeds, UK.He received a PhD in biogeochemical modelling in 2013 from the University of East Anglia.His work focuses on the evolution of Earth's surface environment over geological time, and specifically on the processes controlling planetary surface temperature and marine redox state.His main contribution is on developing novel theoretical frameworks to test hypotheses in geochemistry and evolution.Yannick Donnadieu is CNRS research director at the Centre Européen de Recherche et d'Enseignement des Géosciences de l'Environnement (CEREGE), Aix-en-Provence, France.He is a paleoclimatologist with strong skills in Earth System modelling who works at bridging theoretical models with observations and measurements.His work aims to quantify primary and secondary feedbacks occurring within the Earth System.His research has been devoted to improving our understanding of processes controlling atmospheric CO 2 concentration at geological timescales and to identify ocean atmosphere interactions occurring on Earth as a result of plate tectonics.Yves Goddéris is CNRS research director at the laboratory Géosciences Environnement Toulouse, Toulouse, France.He received his PhD from the University of Liège in Belgium in 1997.He works on the geological evolution of global biogeochemical cycles, focusing mainly on the behaviour of the continental weathering system in the context of plate tectonics.He is also involved in the modelling of weathering reactions in response to anthropogenic forcing of the carbon cycle and climate.

Table 1
Model chemical reservoirs.

Table 2
Model reservoir mass balance.

Table 3
Model flux definitions.

Table 5
Model parameters.