What Models Tell Us About the Evolution of Carbon Sources and Sinks over the Phanerozoic

The current rapid increase in atmospheric CO 2 , linked to the massive use of fossil fuels, will have major consequences for our climate and for living organisms. To understand what is happening today, it is informative to look at the past. The evolution of the carbon cycle, coupled with that of the past climate system and the other coupled elemental cycles, is explored in the field, in the laboratory, and with the help of numerical modeling. The objective of numerical modeling is to be able to provide a quantification of the processes at work on our planet. Of course, we must remain aware that a numerical model, however complex, will never include all the relevant processes, impacts, and consequences because nature is complex and


INTRODUCTION
The ongoing rise in atmospheric CO 2 , sustained by the massive release of fossil carbon by human activities, is the cause of ongoing and future climate change (IPCC 2021).Understanding the link between the carbon release and the CO 2 concentration in the atmosphere requires an accurate description of the global carbon cycle, and of other coupled elemental cycles, such as those of oxygen, calcium, and nutrients such as phosphorus, nitrogen, and iron.One of the main challenges for the scientific community is the prediction of the future evolution of Earth's climate and of the carbon cycle.As there are no data for the future, past climates have been explored to validate methods of prediction of the future climate (Tierney et al. 2020).Most efforts have been concentrated on Quaternary climate (between the present and 2.58 million years ago, Ma), with the setup of international intercomparison programs such as the PMIP (Paleoclimate Modelling Intercomparison Project) initiative (Braconnot et al. 2007).
There is also value in going further back in time, exploring pre-Quaternary, or deep-time, climates.Models have been used to reconstruct the climate up to about 800 million years in the past (the Neoproterozoic) and even to explore the climate of the Archean.There are multiple reasons why a few pioneers started to run numerical models in the early 1980s to reconstruct deep-time climate (Barron 1983).One of the arguments often put forward is that the climate of the last million years was generally colder than the current and the future climate.Thus, even though it contained occasional warm interglacial events, the Quaternary might not be the most appropriate time to constrain climate models.
If we are looking for periods characterized by a warmer climate than the present one, we have to look to deep time (Tierney et al. 2020).Sedimentological studies reveal a plethora of deep-time examples of warmer climates, with flat equator-to-pole temperature gradients and without icecaps.The Cenomanian-Turonian period of the Cretaceous, which is considered to be the warmest climatic episode of the Phanerozoic, was among the first deep-time periods for which the climate was numerically simulated (Barron 1983).
When trying to understand the climatic consequences of an increase in CO 2 , we cannot limit ourselves to an approach that would consist of imposing a CO 2 level and simulating the corresponding climate because the climate changes caused by an increase in atmospheric CO 2 modify the processes involved in the global carbon cycle.These changes to carbon fluxes can amplify the CO 2 increase and thus act as positive feedbacks.Others can generate negative feedbacks within the carbon cycle, limiting the change in CO 2 concentration.The evolution of climate and that of atmospheric CO 2 must therefore be considered as intimately coupled.
The complex processes that link the climate to the carbon cycle require the implementation of numerical models in order to quantitatively estimate the importance of each of these processes.It is not currently possible to accurately reconstruct the world as it was tens or hundreds of millions of years ago, so the target for these models is to estimate the sensitivity of the Earth system to the various processes involved in the carbon cycle and in the other coupled elemental cycles.
In this review, we describe the main processes that control the evolution of the carbon cycle and climate on geological timescales (10 5 years or longer).We show how these processes are numerically modeled, emphasizing the contribution of numerical modeling to the understanding of the evolution of Earth's climate, as well as the limitations of current models.

TIMESCALES AND RESIDENCE TIMES
We define the surficial Earth system as the coupled atmospheric, oceanic, and biospheric surficial envelopes of Earth.These three envelopes contain about 40,000 gigatons of C (GtC), with the ocean being by far the biggest component (38,000 GtC).Carbon is rapidly cycled through these reservoirs.For instance, an atom of carbon typically stays in the atmosphere for less than 3 years before being transferred to the ocean or to the continental biomass.This mean residence time rises to a few thousands of years for the oceanic reservoir of carbon.In the following discussion, we arbitrarily fix the lower bound of the geological timescale to be 10 5 years.
The surficial Earth system is not a closed system; it is connected to the solid Earth.An atom enters the surficial Earth system following one of two distinct pathways: 1. Release from continental rocks, either through weathering of carbonate rocks or through oxidation of organic carbon-rich sediments outcropping at the surface of the continents.
Oxidative weathering of sedimentary organic carbon is estimated to be about 0.08 GtC/year (Lenton et al. 2018), and carbonate weathering at 0.14 GtC/year (Gaillardet et al. 1999).2. Degassing from the solid Earth, including arc volcanism, metamorphism, and magmatic processes along mid-oceanic ridges.The global degassing flux is estimated at 0.07 ± 0.04 GtC/year (Gerlach 2011).
Carbon is then removed by two processes: 1. Carbonate deposition on the seafloor (including sedimentary carbonates and production of carbonate minerals by the weathering of the oceanic crust).Global sedimentary carbonate deposition is estimated to be close to 0.24 GtC/year.The consumption by seafloor weathering is more poorly constrained, estimated at a maximum value of 0.03 GtC/year (Coogan & Gillis 2013).2. Sedimentary organic carbon accumulated in the continental environment and on the seafloor.This flux is of the order of 0.1 GtC/year (Lenton et al. 2018).
The residence time of carbon (the carbon content of a system divided by the sum of the output fluxes, at steady state) within the surficial Earth system is close to 150,000 years, a much longer timescale compared to the time required to mix the surficial Earth reservoirs.However, although 150,000 years might appear to be a long time period, it is short enough to allow huge fluctuations of the surficial Earth system at the geological timescale.Indeed, within a million years, the surficial Earth carbon reservoir can be renewed at least 6 times.In other words, a persistent difference of only 0.1 GtC/year between the global output flux and global input flux could cause Earth's surficial carbon stock to fall to zero within 400,000 years, with dramatic consequences for life on Earth.Compared to the hundreds of million years of the Phanerozoic (the past 539 million years of Earth history), 400,000 years is a short duration.It is also worth remembering that human activities release around 10 GtC/year into the atmosphere, which is huge compared to the geological carbon cycle disequilibria that is discussed in this contribution.
The consequences of this short residence time on the evolution of Earth's surficial environment have been identified by Garrels & Perry (1974) and nicely illustrated by Westbroeck (1991).At the geological timescale, the carbon stock of the surficial Earth system is a tiny reservoir.This reservoir is subjected to geological forces, which may disequilibrate it in such a way that it could either be emptied or grow massively within 1 million years.Both cases would have disastrous consequences for the biosphere: a global-scale glaciation in the first case or the development of excessively high temperatures in the second.The persistence of life on the surface of the planet requires that the carbon content of the surficial Earth system is always close to the equilibrium state.Because the outgassing of the solid Earth is an external forcing to the surficial Earth system, the only solution is the existence of a sensitivity of the outflows of carbon to the carbon content itself, allowing for a negative feedback loop.
Fortunately, there is one process operating at Earth's surface that is able to stabilize the atmospheric CO 2 level and hence the climate: the chemical weathering of exposed silicate minerals.Once exposed to the physical conditions of Earth's surface, silicate minerals are not thermodynamically stable, causing them to dissolve.The dissolution is promoted by the presence of carbonic acid (from dissolved CO 2 ) in rainwater and in soil solutions.The dissolution reaction can be written as follows (showing dissolution of a plagioclase by carbonic acid, release of cations in the rivers, and precipitation of secondary clay mineral): Although silicate rocks do not contain carbon, their dissolution generates alkalinity (e.g., HCO 3 − ), which is carried to the ocean by rivers, with the carbon coming from the atmosphere.When this alkalinity reaches the ocean, the saturation state of the ocean water with respect to carbonate minerals increases, generating an increase in calcium carbonate production by calcifying organisms to remove the alkalinity excess.In the above equation, two calcium ions are available for calcium carbonate precipitation and potential storage in the sediments.Each of those calcium ions takes with it 1 mole of carbon, originating from the atmosphere.Four are left behind in the ocean.It is therefore important to understand that the weathering of silicate rocks consumes CO 2 on geological timescales only if the process of calcium carbonate deposition in the oceans is functioning and can respond to the arrival of excess alkalinity.

FROM THE WALKER NEGATIVE FEEDBACK TO THE REAL WORLD: 40 YEARS OF DEBATE
As early as 1845, about 40 years before the first description of the greenhouse effect by Arrhenius (Maslin 2004), a French engineer named Ebelmen foresaw the need for some sort of coupling between the consumption of carbon dioxide by the alteration of rocks and its return via magmatic degassing, in order to stabilize the carbon cycle (Ebelmen 1845, Berner & Maasch 1996).It was only in 1981 that Walker and collaborators proposed the existence of a negative feedback reaction able to modulate the carbon cycle in order to stabilize the climate of Earth.
The Walker negative feedback is a simple and elegant model to explain the relative stability of Earth's climate.Walker et al. (1981) proposed the existence of a negative feedback allowing a control on the atmospheric CO 2 evolution in response to the increase of the solar constant (brightening of the Sun) from the time of Earth's formation 4.5 billion years ago to the present day.They explicitly stated that the time evolution of atmospheric CO 2 has also been impacted by other important factors, such as the paleogeography, the planetary albedo, and possible influence of Overview of the feedback loop and causal links between the various component of the surficial Earth system.The green lines represent the silicate weathering negative feedback.A green arrow means a positive impact; a green dot represents a negative impact.The red arrows show the multiple connections between the various components of the geological evolution of Earth.
land plants on weathering rates.Walker et al. (1981) proposed the addition of a cog in a complex machine.They demonstrated that this additional cog is able to delay the evolution of Earth's surface temperature in response to the increase of the solar constant.They never claimed that this additional gear is the main controlling factor of the climate evolution of Earth.However, this feedback is nested in the core of the carbon cycle machinery and consequently plays a nonnegligible role (Figure 1).
We can now state the two laws of the Earth-thermostat controlling the geological evolution of the carbon cycle: 1.The sum of the carbon fluxes leaving the surficial Earth system must always be very close to the carbon flux entering the surficial Earth system, to avoid lethal fluctuations.2. To achieve this, negative feedbacks are required.The weathering of continental silicates can play this role.But there may be other solutions.It is worth mentioning that the existence of the silicate weathering feedback has been demonstrated from field data or data analysis at various timescales: from the decadal to the million years (Zeebe & Caldeira 2008, Gislason et al. 2009, Arnscheidt & Rothman 2022).
The first law cannot be violated because it is the mathematical consequence of the short residence time of carbon in the surficial Earth system, compared to the geological timescale.The second law is less constraining, and it allows us to consider new solutions rather than just the climate-dependent silicate weathering.Most numerical models use a kind of Arrhenius type of equation to calculate the CO 2 consumption by silicate weathering (Dessert et al. 2001, Oliva et al. 2003): 1.
This law has been calibrated on field measurements.E a is an apparent activation energy, generally around 40,000 J/mol.R is the ideal gas constant.Considering the data from which this equation has been calibrated, T is the mean annual temperature at the surface (in K), and runoff is the mean annual runoff (precipitation − evapotranspiration). k is a calibration constant that is fixed to fit the data.Area represents the surface of the considered continental zone: It could be the surface of a small watershed but also the surface of the continents.Given that runoff and temperature are functions of the atmospheric CO 2 , this equation can be also written as a function of the CO 2 partial pressure (Caves-Rugenstein et al. 2019): where R CO2 is the past partial pressure of CO 2 , normalized to the preindustrial value.

NUMERICAL MODELS
We build numerical models when the system of interest is too complex to predict its evolution in response to external or internal changes.Numerical models are also a powerful method to evaluate the relative importance of a specific process within the carbon cycle.In this contribution, we focus only on models integrating the coupled evolution of the climate and of the carbon cycle (potentially coupled to other biogeochemical cycles) over geological timescales.All these models belong to a category called box models.A box model decomposes the Earth system into a variable number of boxes connected to each other by flows of chemical species.The boxes are assumed to be well mixed and homogeneous.The number and the size of the boxes depend on the spatial and temporal scale of the study.This type of model allows for very long simulations, with simulated times easily exceeding 1 million years.
In a box model, the time evolution of a given reservoir i containing C i moles of carbon (for instance, the carbon content of the atmosphere) is the solution of the mass balance between the input F i in and the output F i out fluxes.All these fluxes may be a function of time and a function of the size of one or several of the n simulated reservoirs: The system is said to be at steady state when all the time derivatives are equal to zero, meaning that the input of elements to a given reservoir equals the output.
For the present discussion, we consider a simplified system with two boxes: one box with a short residence time (7 years), similar to the residence of carbon in the atmosphere, and one with a longer residence time (500 years), close to that of the ocean (Figure 2).We denote the short residence time reservoir as atmosphere and the other as LTSR (for Long residence Time Surficial Reservoir).These two boxes exchange carbon.We assume that these exchange fluxes are proportional to the size of the reservoir from which the fluxes come, while the fluxes entering a given reservoir do not depend on the carbon content of that reservoir.At steady state, the inverse of the proportionality constant is equal to the reaction time of the reservoir relative to the considered flux.In addition, the reservoir with the shorter residence time is fed by a geological flux (if this reservoir is the atmosphere, the geological flux can be the solid Earth degassing, for instance).A schematic description of a simple 2-box model.C atm and C LTSR are the atmospheric content of carbon, and the amount of carbon stored in a Long residence Time Surficial Reservoir, respectively.All the exchange fluxes are assumed to be proportional to the reservoir where they come from.The k are the proportionality constants.F degassing is a forcing flux, such as the solid Earth degassing.F out is a flux leaving the surficial reservoirs and transferring carbon to very long-term storage, such as the oceanic sediments.This theoretical model is used here as an illustration of the key role played by the residence times.This should not be considered as a careful representation of the carbon cycle.
For this experiment, we set this geological flux at 0.1 GtC/year.To avoid accumulation of carbon inside the two connected reservoirs, we also add an output flux that may potentially be interpreted as the sedimentary sink of carbon.This flux is set proportional to the size of the LTSR.At steady state, this flux is equal to 0.1 GtC/year.Our simple model is not a realistic simulation of the real world, but the conclusions we draw from it are generally applicable.
The system is first assumed to be at steady state, meaning that the content of carbon in the two reservoirs is constant through time, with the input flux of each reservoir balanced by the output flux.Then, we impose on the system a progressive rise in the outgassing of carbon, from its background value of 0.1 GtC/year to 0.15 GtC/year.Such an increase is far from being negligible: From the available reconstructions of the past degassing rate of the solid Earth, an increase by 50% never occurred in less than at least a few million years (Goddéris & Donnadieu 2019).
When this rise in degassing occurs within a human timescale of 50 years, the maximum disequilibrium in the content of carbon of the atmosphere is calculated at 0.008 GtC/year, while the LTSR disequilibrium culminates at 0.049 GtC/year, which is 10% of the increase in the degassing rate.Those values rapidly decrease to 0.00006 GtC/year and 0.004 GtC/year when the change in the outgassing occurs at a geological timescale of 5 million years.The last value is less than 1% of the change itself.
Based on these calculations, a useful simplification is to set the derivatives of the surficial reservoir to zero when simulating the evolution of the carbon cycle at the geological timescale.This assumption allows us to change the differential equations describing the evolution of the carbon cycle at the geological timescale into a set of simple algebraic equations.This simplification is applied in several models, such as GEOCARB and its later versions (Berner 2004(Berner , 2006) )

(Table 1).
A full demonstration of the validity to set the derivatives to zero, as long as the forcings are geological, can be found in the appendix by François & Goddéris (1998), based on a realistic representation of the carbon and alkalinity cycles.

THE RAYMO TURNING POINT
Based on the Earth-thermostat principle, Berner et al. (1983) built the first numerical model able to calculate the evolution of atmospheric CO 2 levels over the past 100 million years.This model was later termed BLAG after the authors' surnames, and it incorporates the negative feedback of silicate weathering.The BLAG model was mainly driven by the CO 2 outgassing input of carbon into the surficial Earth system.The role of silicate weathering was limited to a stabilization process able to limit the changes caused by variations in the outgassing rate.Although the role of silicate weathering is crucial, it is rather passive within this source-driven framework.
But the hypothesis that the carbon cycle is mainly driven by the sources of carbon was challenged in the late 1980s and early 1990s.Maureen Raymo and coauthors published several articles that can be seen as a turning point in our understanding of the geological evolution of the carbon cycle (Raymo et al. 1988, Raymo 1991).The reasoning of Raymo and colleagues is based on the carbonate sediment record of the 87 Sr/ 86 Sr ratio in seawater (Veizer et al. 1999).This record was first published in 1982 by Burke and coauthors (1982).A more recent version is available by Veizer et al. (1999), and compilations continue to add detail (McArthur et al. 2020), yet the fundamentals of the Cenozoic curve are unchanged.The most salient pattern in this curve is a steep rise over the past 40 million years or so.Because strontium has a high isotopic ratio in continental rocks and a low ratio in hydrothermal sources, Raymo (1991) concluded that this isotopic signal is a record of rapidly increasing chemical weathering of the continents, in the face of hydrothermal processes that remain at a low level.This rapid increase is attributed to the uplift of the Himalayas and the Tibetan plateau.Increasing slopes and decreasing temperatures with altitude promote physical erosion.A more fragmented rock presents a larger contact surface to the circulating fluids, which accelerates the dissolution of minerals and thus the consumption of CO 2 .Tectonic uplift increases the global silicate weathering rate, consuming CO 2 , cooling the climate, and favoring the onset of glaciers in some regions.Those glaciers will, in turn, increase the overall physical erosion of the continents, further consuming CO 2 .But this hypothesis breaks the first law of the Earth-thermostat.Indeed, it introduces a positive feedback between climate cooling and silicate weathering.And when the first law is no longer applicable, with sinks greater than sources, CO 2 would disappear in a few million years, leading to extreme cooling and potentially a snowball glaciation.
When considering the intense debates of this period, it is clear now that there was a misunderstanding of the functioning of Walker's theory.The existence of a positive feedback loop within the carbon cycle does not exclude the presence of a negative loop.The two loops will interact, one destabilizing the climate, the other rebalancing it.The relative effectiveness of these two loops will allow for greater or lesser fluctuations in atmospheric CO 2 content and climate.Thus, nothing prevents erosion from increasing locally in response to tectonics, nor alteration from decreasing in other regions of the world as a result of global cooling (Goddéris & François 1995).The negative feedback loop, despite tectonic forcing working to destabilize it in some regions, has the ability to maintain the global flux of CO 2 consumption by silicate weathering at a level close to volcanic degassing.If degassing remained constant over the Cenozoic, as suggested by some reconstructions, then the global flux of weathering should also remain constant.

WEATHERING AND WEATHERABILITY
How can silicate weathering change the partial pressure of CO 2 if the sink of carbon through weathering must be close to the degassing rate?The answer is by changing the weatherability.The concept of weatherability was defined by Kump & Arthur (1997, p. 422).It is "the product of factors that affect rock chemical weathering other than climatic factors" (e.g., air temperature and runoff ).As reported by Caves et al. (2016), the dependence of the global silicate chemical weathering on atmospheric CO 2 can be expressed as a power function of the CO 2 level: n ranges from 0.2 to 0.6 in most models [e.g., LOSCAR (Long-term Ocean-atmosphere Sediment CArbon Reservoir) (Zeebe 2012), GEOCARB (Berner 2004), COPSE (Carbon, Oxygen, Phosphorus, Sulphur and Evolution) (Bergman et al. 2004, Lenton et al. 2018)].R CO2 is the ratio between the CO 2 level at any time in the past and its present-day value.The power function should account for all the direct climatic effects (impact of temperature and runoff changes on silicate weathering, which both depend on the atmospheric CO 2 level).k is the weatherability.
Given that the carbon cycle must be close to its steady state, we can write Combining Equations 4 and 5, we get the following simple relationship: If the solid Earth degassing rate F vol is known as a function of time, if k is known, and if n can be constrained, this simple equation can be solved to calculate the atmospheric CO 2 level.This equation is present in some form in models such as GEOCARB, LOSCAR, MAGic (Arvidson et al. 2006), and COPSE.Although this mathematical relationship is simple, it includes a large number of complex processes.n can fluctuate with age, depending on vegetation cover, for example.The weathering constant k depends on the presence or absence of mountain ranges.
Finally, the rate of degassing of the solid Earth remains difficult to constrain in the past (Goddéris & Donnadieu 2019).
In summary, with respect to the evolution of Earth's climate during the Cenozoic, the uplift of large mountain ranges (such as the Himalayas) can modify the ability of rocks to weather under the effect of physical erosion, which increases with the presence of slopes, glaciers, and freezethaw alternation.Many authors have proposed an increase in the weatherability of continental surfaces (k) as the cause of the decrease in CO 2 content over the same period.The central question is then by how much weatherability increased over the past 40 million years.The simplest way to evaluate this is to use reconstructions of sedimentary masses accumulated on the seabed from the erosion of continents (Wold & Hay 1990), but this is still a rough method.It is obvious that the weatherability did not change at the same time and with the same amplitude everywhere on the continents.Some regions were affected by uplift but not the global surface as a whole.Although a simple model's response to the global increase in physical erosion will be quite strong and unequivocal, this method is not entirely convincing.The seawater 87 Sr/ 86 Sr record has also been used to constrain the physical erosive impact on k in GEOCARB (Berner 2006).

WEATHERABILITY AND ISOTOPIC SIGNALS
The short residence time of carbon in the surficial Earth system implies that silicate weathering rates remain close to solid Earth CO 2 emissions.If this degassing rate stayed roughly constant, then so should the silicate weathering flux.But if these fluxes were indeed both roughly constant, how can we explain the rise in the 87 Sr/ 86 Sr of the seawater over the past 40 million years?The answer here came from field measurements and models (Edmond 1992, Goddéris & François 1995, Galy & France-Lanord 1999): What is changing over the Cenozoic is not the global Sr flux generated by continental weathering but its isotopic signature.Highly radiogenic rocks are exposed in the Himalayas, and their dissolution brings radiogenic strontium to the ocean, forcing the oceanic 87 Sr/ 86 Sr to rise, while the global silicate weathering flux may stay constant.
A similar ubiquity can be suspected in the interpretation of the 7 Li/ 6 Li of the ocean, although this is more complex because of the fractionation of lithium isotopes during secondary phase precipitation.This fractionation implies that the riverine signature is dependent on the clay precipitation efficiency and is not strictly a signature of the rocks being weathered, as is the case for Sr.But again, the 7 Li/ 6 Li of the ocean can be explained by a change in the riverine 7 Li/ 6 Li together and/or a change in the weathering flux itself (Misra & Froelich 2012).Once the Li cycle coupled to the carbon mass balance, several combinations are consistent with the rise in seawater 7 Li/ 6 Li: increasing riverine 7 Li/ 6 Li together with a decreasing Li flux or the reverse (Vigier & Goddéris 2015).
A constancy of global weathering has been suggested for the past 12 million years of the Cenozoic, based on the relatively constant 9 Be/ 10 Be ratio of marine sediments (Willenbring & Von Blanckenburg 2010).It was not correctly interpreted, as Willenbring & Von Blanckenburg (2010) concluded that silicate weathering rates were constant and therefore not the driver of Cenozoic climate cooling.But Caves Rugenstein et al. (2019) took up this discussion and correctly interpreted that the constancy of the Be signal was a signal of constant global alteration, in a context of increasing weatherability.

SPATIALIZATION
A striking feature of most numerical models simulating climate and CO 2 evolution at the geological timescale is that they do not account for spatialization-the difference in processes across different parts of the land surface.These are mostly global-scale models, and they consider continental weathering as a single global value, driven by a mean global temperature and a mean global runoff (e.g., in GEOCARB, LOSCAR, and COPSE), and by global weatherability.A posteriori, this methodology may appear surprising, given that many mathematical relationships included in these models are nonlinear functions of the climatic parameters.This is critical when considering the climatic parameters of importance for silicate weathering.In the global-scale models, if CO 2 rises, global average temperature will also rise, as will global runoff.The continents are simulated as a single grid element, where silicate weathering can increase with CO 2 , only if weatherability is held constant.But the real world response of runoff to a change in CO 2 is highly complex.Globally, runoff will probably increase with atmospheric CO 2 .But locally, it may decrease, and this will be dependent on the local relief.Another important factor is the paleogeography-the positioning of the continents on Earth's surface.The global runoff pattern will be strongly influenced by changes to paleogeography.
The first attempts to include spatialization of climate and silicate weathering in deep-time models of the carbon cycle date back to the late 1980s.Marshall et al. (1988) used an Energy Balance Model (EBM) to simulate the climate that divided Earth into latitudinal bands.The principle was to simulate climate variations as a function of latitude.The weathering rates were then calculated in each latitude band and summed to produce a global value.Using a climate model of the same type, François et al. (1993) built a numerical model coupling the global cycle of carbon, oxygen, phosphorus, and alkalinity to an EBM with 18 latitudinal bands.This work represented a significant advance in our understanding of the global evolution of Earth's surface environment.It is worth keeping in mind that, during the same period, GEOCARB II was still using very simple laws linking global air temperature and global runoff, respectively, to the atmospheric CO 2 partial pressure and the global average temperature (Berner 1994).However, despite this improvement of the calculation of the climatic fields using a more mechanistic model, EBMs do not resolve the water cycle.François et al. (1993) developed an equation allowing them to calculate the runoff for each latitudinal band as a function of the latitude and the calculated air temperature.
Taking a further step in this direction, Donnadieu et al. (2006) coupled a dynamic model of the main biogeochemical cycles (ten boxes, one of which describes the atmosphere and the other nine the ocean) to a General Circulation Model (GCM)-type climate model.An Atmospheric General Circulation Model (AGCM) is a complex numerical model describing the transport of matter and energy in the atmosphere in three dimensions.AGCMs can be coupled to a simple representation of the surface ocean that simulates the transport of heat as a diffusive process only-called a slab ocean.A more complete physical representation of the ocean can also be implemented, with diffusion and transport of heat and matter within the ocean.These models are called Atmosphere-Ocean General Circulation Models (AOGCMs).They include a physical representation of the water cycle and can be used to calculate continental runoff on a 2D-spatial grid covering the continents, together with the spatially resolved air temperature.
The model developed by Donnadieu et al. (2006) is called GEOCLIM, and it has been used for a large variety of studies.It was found that supercontinental configuration (for example, the Pangaea supercontinent that existed during the later Paleozoic and into the Mesozoic) would strongly limit the impact of changes to CO 2 on continental runoff because of the development of arid conditions over large continental surfaces.These periods are thus characterized by high steady-state CO 2 levels.Conversely, a wide dispersion of the continents should facilitate an increase of runoff with CO 2 (Goddéris et al. 2014).In addition to this, the climate sensitivity to atmospheric CO 2 (e.g., global warming per doubling of CO 2 concentration) is dependent on the continental configuration, ranging from about 4°for a CO 2 doubling during times of supercontinental configuration, or when large continental masses are located on the South Pole (such as in the early Paleozoic), to 2°for the present-day configuration.Finally, still using GEOCLIM, Goddéris & Donnadieu (2019) found that the paleogeographic control of the evolution of the carbon cycle and the climate is more pronounced in the Paleozoic.In the Mesozoic and Cenozoic, the continental dispersion appears to be large enough to force the Earth system into a glacial stage since the breakup phase of Pangaea in the Jurassic.As this was not the case, higher solid Earth degassing is required to keep CO 2 at a high level, especially during the Cretaceous.
Following the method employed in GEOCLIM, the SCION (Spatial Continuous Integration) model (Mills et al. 2021) extends the COPSE global biogeochemical model (which itself is a modification of GEOCARB) to include a spatial land surface informed by the FOAM GCM, which is interpolated over time as the continental configuration changes.The effects of continental distribution on the modeled atmospheric CO 2 evolution are similar to those observed in GEOCLIM, and the model predicts higher CO 2 levels during the time of the Pangaea supercontinent when compared to COPSE.An advantage of this new model is that the continuous whole-Phanerozoic simulations typical of COPSE and GEOCARB can now be run using a spatial representation of climate.
Adding a 2D representation of the continental surfaces allows for the simulation of the weathering fluxes accounting for the lithology.Richey et al. (2020) explored the sensitivity of the late Paleozoic climate to the spatially distributed weathering of mafic outcrops of variable size.Park et al. (2020) simulated the past 15 million years of global carbon cycle evolution, accounting for the presence of ultramafic rocks outcropping in a humid tropical Southeast Asia.Using an improved version of GEOCLIM, they found that the weathering of these ultramafic rocks is responsible for a global cooling of 3°C since 15 million years ago (Park et al. 2020).Reconstructing lithology further back into deep time becomes more challenging.At least large igneous provinces (LIPs) can be located on a map and assigned a lithological class (basalts).This has been done for the entire Cenozoic (Lefebvre et al. 2013) and for the late Proterozoic (Goddéris et al. 2017b).In both cases, the LIP weathering left its mark on the climatic evolution of Earth by reducing atmospheric CO 2

VEGETATION
The role of continental vegetation has been included in numerical models of the deep-time carbon cycle mainly because of biotic impacts on silicate weathering rates.Studies of the weathering processes in modern environments have shown that the presence of plants may increase weathering rates due to the belowground accumulation of CO 2 , the release of organic acids, and the mechanical effect of roots (Drever 1994, Moulton et al. 2000, Porder 2019).This biological enhancement of global weathering reactions has been added in numerical models of the deeptime carbon cycle by modifying the weatherability factor.In nondimensional models, such as GEOCARB, the weatherability constant k is divided by 4 prior to the advent of land plants (representing the average effect seen in field studies).In the GEOCARB model (Berner 1994), such a drastic decrease in the weatherability constant generates very high CO 2 levels in the early Paleozoic.To avoid this, a direct dependence of global weathering rates on atmospheric CO 2 was added, which was suggested to simulate the effect of the lower pH of rainwater (Berner 2004).GEOCARB-type models do not account for the change in hydrology associated with plant colonization of the land because they do not resolve the water cycle explicitly.Ibarra et al. (2019) built a cascade of models: a 1D model of water vapor transport feeding a reactive transport model for weathering to simulate the impact of the rise of the vascular plants in the Devonian.They found that the colonization of the continents by vascular plants has increased the global weatherability by a factor of about 2, with half of this increase being related to the hydrological effect triggered by the enhanced transpiration on the continents, allowing water vapor to be transported further into dry continental interiors.
Le Hir et al. (2011) used GEOCLIM coupled to a dynamic vegetation model to explore the global changes during the progressive colonization by vascular plants during the Devonian.This study was performed without accounting for a direct weatherability increase.The objective was to quantify the impact of plant-driven changes in the continental albedo, rugosity, and evapotranspiration on the global water cycle.They found that the decrease in albedo warms the continents enough to significantly increase the continental runoff.This change in runoff promotes weathering reactions, reducing the partial pressure of atmospheric CO 2 .The calculated drop in CO 2 was from 7,000 ppm to 2,000 ppm between the Early and Late Devonian.Surprisingly, the air temperature at ground level remained about constant at around 22°C, the cooling generated by the CO 2 decrease being compensated for by the decrease in albedo.
There are currently no modeling studies that incorporate a dynamic, spatially resolved terrestrial biosphere evolving through the Paleozoic Era that also include the direct effect of a weatherability increase.Although the SCION model uses a spatially resolved climate (from GEOCLIM) and imposes a global weatherability increase due to plant evolution (like COPSE and GEOCARB), the weatherability change is simply applied to every continental grid cell regardless of habitability for plants.Future work will aim to link a dynamic global vegetation model (e.g., Gurung et al. 2022) in a similar manner to that of Le Hir et al. (2011).
Finally, Maffre et al. (2022) explored the impact of the plant colonization process on physical erosion and chemical weathering.They generated several random scenarios for the colonization phase.When a grid cell is colonized by vascular plants, chemical weathering increases and physical erosion is also impacted, either decreasing or increasing.They found that some scenarios generate a gentle CO 2 decrease, as continents are being progressively colonized.But other scenarios generate first an ample increase of atmospheric CO 2 .This is specifically the case for scenarios where the plants decrease physical erosion, promoting the growth of the regolith, which prevents the percolating water reaching the bedrock (see the next section).
There is currently a major effort to build process-based models of paleo-plants function to explore the interaction between ancient plants and their environments (Matthaeus et al. 2023).Once integrated into coupled ecological-climate models, these tools pave the way toward the reconstruction of the interaction between plants and climates in deep time.Those models are currently restricted to snapshot simulations.

EROSION
Since the early 2000s, the question of a possible control of the chemical weathering of rocks by physical erosion became increasingly widespread.From the field point of view, there is clearly a correlation between physical erosion and silicate rock weathering, as shown by Millot et al. (2002) and West et al. (2005).Several theoretical models have been developed (Gabet & Mudd 2009, West 2012, Maher & Chamberlain 2014).Those models focus on the silicate weathering sensitivity to climate, under various climatic conditions.For instance, Maher & Chamberlain (2014) showed that the sensitivity of silicate weathering to climate change is higher in areas where the residence time of water in soils is low (e.g., mountainous areas).Gabet & Mudd (2009) built a numerical model of weathering that includes the role of the thickness of the saprolite-a major contribution.Indeed, when thinking of physical erosion, we focus on places where this physical erosion is massive, such as in young and active mountain ranges (for instance, in Taiwan or New Zealand).As a result, a biased picture of continental surfaces seems to emerge from these studies in highly erosive areas.Outside of those regions, the continents look like vast flat surfaces where nothing happens.
However, the flat areas are indeed characterized by low weathering efficiency (low weatherability) because the saprolite is so thick that water cannot reach the fresh mineral at the saprolite-bedrock interface.Such systems are transport limited because the dissolution rates and the transport of solutes are limited by the capacity of the water to find its way into thick soil cover, typical of flat tropical environments.Using a coupled model describing physical erosion and weathering with a 2D spatial resolution, Carretier et al. (2014) showed that these flat areas weather 5 to 100 times slower than active orogens.Considering the present-day geography and the extent of flat areas, it is possible that flat areas consume as much CO 2 as the active orogens.But this needs to be explored.Also, for some large watersheds, the eroding mountains can provide fine-grained materials, without dissolving them, because the temperature is too low and because the residence time of water is too short.Those erosion products are then transported into the floodplains, where the temperature and runoff conditions promote weathering.This is the case in the Amazon, for example, where 50% of the CO 2 consumption occurs in the Andes, while the remaining 50% occurs in the floodplains (Moquet et al. 2011).
Up to now, the integration of physical erosion in large timescale box models such as GEOCARB or COPSE has been done by modifying the weatherability, with a corrective factor.As those models are global, lacking spatialization, this change in weatherability is applied to the entire continental surface.More importantly, the erosive factor is an external forcing parameter.This parameter is usually taken to be related to the total amount of sediments released in shallow seas over the Phanerozoic (Wold & Hay 1990).But physical erosion is not a forcing.Indeed, erosion is a function of the local slope but also of local climatic conditions (in particular runoff ) (Syvitski & Milliman 2007).Goddéris et al. (2017a), for the first time, accounted for physical erosion in the global-scale numerical model GEOCLIM.Here, physical erosion is calculated over the spatially resolved 484 Goddéris continental surfaces.The mass balance of the saprolites is calculated on each continental grid cell, depending on slope and runoff.The saprolites are removed at the surface by physical erosion and created at the bottom by the dissolution of the bedrock.The dissolution rate of bedrock is a function of runoff and temperature, as usual, but modulated by a function decreasing with the thickness of the saprolite.This function simulates the progressive shift from a weathering-limited to a transport-limited regime, when the thickness of the saprolites increases.The most important thing here is that the physical erosion is updated at each time step using the evolving climatic fields from a GCM.The main difficulty is the definition of the slopes for each grid cell, which is poorly constrained in the distant past but required for the calculation of physical erosion.According to Goddéris et al. (2017a), the slopes are fixed for each grid cell and calculated from a correlation between altitude and slope observed for the present day.The problem therefore becomes the knowledge of topographic height in the geological past, which remains to this day a real challenge.

THE NEW MOUNTAINS
Our understanding of the carbon balance of mountains has considerably evolved since the hypothesis of Raymo (1991).Two Copernican evolutions occurred, the first one in 1997 with the publication by France-Lanord & Derry (1997).From measurements of suspended solids in the Himalayan rivers and in sediments of the Bengal Fan, they showed that the very high erosion rates generate an extremely efficient burial of organic matter on the seafloor.This organic matter comes either from tropical forests that cover the southern foothills of the Himalayas or from marine primary production.This flux is estimated to store two to three times more carbon than the silicate weathering flux.Ten years later, Galy et al. (2007) came to the same conclusion.
The second revolution came in 2014.Torres et al. (2014) identified a new process particularly active in the Himalayas but that was initially identified by Spence & Telmer (2005) in Northwest Canada.In some mountainous areas, physical erosion exposes pyrite-rich outcrops.In contact with the atmosphere, the pyrites oxidize and release sulfuric acid.This acid causes the dissolution of carbonates, which releases carbon from the rock to the rivers, and ultimately the carbon from the rock degasses toward the atmosphere.From the point of view of this process, the physical erosion of mountain ranges could generate CO 2 .Torres et al. (2016) identified the same mechanism in the Andes.
But the modeling community largely ignored this new process.It was included in the global box model COPSE (in Mills et al. 2014), but the close balance of pyrite weathering and pyrite burial that model maintained completely nullified any long-term effect on atmospheric CO 2 .The absence of dynamic modeling of physical erosion prohibited more detailed calculation of the erosion fluxes that are a controlling factor for sulfuric acid release and for organic carbon burial.More recently, Maffre et al. (2021) simulated the release of sulfuric acid in mountain ranges, using GEOCLIM and its spatialized representation of weathering and physical erosion.In this theoretical study, they released reduced sulfur, which is oxidized and contributes to the dissolution of carbonate and silicate rocks.During the first stage of the numerical experiment, the dissolution of carbonates by carbonic acid releases CO 2 , and the atmospheric CO 2 content rises.But this is only a transient effect.Because oxygen is efficiently consumed during the reaction, the shallow epicontinental seas become progressively less oxygenated, promoting organic carbon preservation.After a short period of elevated CO 2 , its concentration decreases rapidly, after only a few million years, limiting the consequence of this process.

SOLID EARTH DEGASSING
The rate of CO 2 degassing from the solid Earth has been reconstructed in various ways over the Phanerozoic and tentatively into the Neoproterozoic.The question of the outgassing rate of the solid Earth is critical: It is clear that this flux is a first-order factor in the evolution of the carbon cycle.But constraining its amplitude and its time evolution is still a major challenge.One of the problems encountered in the reconstruction of outgassing rates comes from the fact that the global flux is the result of the addition of several components.These include the degassing rate of oceanic ridges (Domeier & Torsvik 2019) and different types of subduction zone volcanism and metamorphism (Lee et al. 2015, McKenzie et al. 2016), as well as continental rifting (Brune et al. 2017).
The first reconstruction was performed by Gaffin (1987), who used the evolution of sea level over the whole Phanerozoic to calculate a volume of the mid-oceanic ridge systems, which tends to be higher when hydrothermal activity is intense.The CO 2 degassing curve was then assumed to be proportional to this activity.The results show globally high degassing in the early Paleozoic and middle Mesozoic.Conversely, the end Paleozoic and the Cenozoic are characterized by low degassing rates.Other methods rely on the inventories of consumed lithosphere at subduction zones (Engebretson et al. 1992) or on the reconstruction of seafloor production rates (Cogné & Humler 2008).Recently, attention has been focused on continental arc volcanism (McKenzie et al. 2016).Visually, the reconstruction of continental arc volcanism over the Phanerozoic seems to be in agreement with the climatic fluctuations of Earth.However, the time resolution of the reconstructed arc volcanism is still quite low.For instance, McKenzie et al. (2016) showed that arc volcanism was limited at the time of the late Paleozoic ice age, but only around the end of the glacial period, and consequently can be considered only as a contributor to the climate cooling but not as the initiator of the ice age.Goddéris & Donnadieu (2019) tested a set of curves proposed in the literature by using the reconstructed histories in the GEOCLIM model.This exercise clearly showed the major impact of degassing on the carbon cycle and the poor constraint on the global outgassing flux.
Recently Müller et al. (2022) published a new CO 2 degassing curve for the Mesozoic and Cenozoic based on full plate tectonic modeling and a method for estimating the carbon content of the downgoing slab.A crucial point in this study is the decoupling of the CO 2 flux released during accretion at ocean ridges from degassing by arc volcanism over the Cenozoic.Müller et al. (2022) demonstrated that this decoupling is related to increased carbonate sedimentation on the pelagic seafloor.Indeed, the biological evolution of calcifiers led to a shift in sedimentary carbonate accumulation sites from the continental shelves in the Mesozoic to the deep ocean floor during the Cenozoic.As a result, ever-larger masses of carbonates were driven into subduction zones, leading to an increase in the rate of carbon outgassing by arc volcanism, which increasingly dominates the outgassing flow.This transfer of carbonate sedimentation had been described as early as 1989 (Volk 1989) and anticipated as a means of compensating for the increased chemical weathering of silicate rocks associated with the uplift of the Himalayas and the Tibetan plateau (Caldeira 1992).This result is particularly interesting and challenges interpretation of the 9 Be/ 10 Be record, which was assumed to imply constant global silicate weathering.Given a rigorous 4D modeling approach, the full plate method can also be extended to include other sources of CO 2 such as the heating of overlying sediments in continental arcs (Lee et al. 2015) and CO 2 sources at rifts (Brune et al. 2017).

SEAFLOOR WEATHERING
The CO 2 consumption by the alteration of seafloor basalts in contact with seawater was first described in the early 1980s.This additional carbon flux is potentially an important process in the carbon cycle evolution (Isson et al. 2019).Like the weathering of continental rocks, this process releases alkalinity into the waters circulating inside the seafloor.This alkalinity is then removed from the seawater via the precipitation of calcium carbonates in the seafloor (Staudigel & Hart 1983, Spivack & Staudigel 1994).Seafloor weathering was recognized early on as a significant carbon sink (Brady & Gislasson 1997).The magnitude of the sink has been estimated at 0.02 GtC/year (Coogan& Gillis 2013), able to play a significant role in the geological evolution of the carbon cycle.François et al. (1993) developed a numerical model of the carbon cycle integrating the role of seafloor weathering in order to understand the evolution of climate over the Cenozoic, mainly based on the seawater 87 Sr/ 86 Sr.The seawater 87 Sr/ 86 Sr has rapidly risen since about 40 million years ago, and this rise was first interpreted as a rapid increase in global silicate weathering (Raymo 1991).But as mentioned before, the sinks and sources of carbon must always be close to steady state.Through the introduction of seafloor weathering, the mass balance of the inorganic carbon cycle moves from solid Earth degassing = continental silicate weathering, which is the first law of the Earth-paleothermostat, into the following expression: solid Earth degassing = continental silicate weathering + seafloor weathering.
If the degassing rate stays constant, a way to increase the continental silicate weathering rate over the Cenozoic is to assume that seafloor weathering has decreased at the same time (Goddéris & François 1995).Because continental weathering provides radiogenic strontium to the ocean, while seafloor weathering and water-rock interaction at mid-oceanic ridges release less radiogenic strontium, the addition of seafloor weathering to the conceptual framework allows an increase in the strontium isotopic ratio of the ocean, without invoking a change in the isotopic composition of the source rocks.
But why would seafloor weathering rates be higher at the end of the Mesozoic compared to today?One possibility is that a higher ridge generation rate in the late Mesozoic allowed for high rates of seafloor weathering (e.g., Gillis & Coogan 2011, Mills et al. 2014).Brady & Gislason (1997) were the first to suggest a strong temperature dependence.Given that the temperature of the deep sea was 10 to 20°C above the present-day value, ocean floor basalts should dissolve faster.But a higher dissolution accompanied by a greater storage of carbonates inside the ocean crust might also be due to a different chemical composition of the oceans (especially the content of K) in the late Mesozoic.Le Hir et al. (2011) tried to account for this flux during the Snowball Earth glacial episode in the Neoproterozoic by simulating the dependence on pH and on the affinity of the seawater with respect to the minerals composing the oceanic crust.They found that seafloor weathering can strongly increase the duration of the total glaciation, given its efficiency under low pH conditions, and due to the accumulation of CO 2 in the surficial Earth system.

CONCLUSIONS
Many questions have been considered using numerical models of the long-term carbon cycle coupled to climatic simulation.But in the end, few of them have been answered definitively.A great suite of models exist (e.g., Table 1) that include different processes and operate over different timescales, usually with models that include more detailed processes also having a more limited timeframe of operation.This means that while meaningful conclusions have been drawn on carbon sources and sinks over Phanerozoic and later Neoproterozoic time, it has not been possible to include all of the relevant mechanisms (at appropriate resolution) into a single model.We strongly believe that numerical models should be used to identify and quantify processes impacting the evolution of the carbon cycle and climate, and not as tools able to reconstruct past environments precisely, for the simple reason that the complexity of nature cannot be fully integrated into a model.In addition, when going back into deep time, we have to acknowledge that most of the data have been lost forever.Short-term fluctuations progressively disappear from the record when going back in time because the time resolution decreases with geological age.Because of this decrease, it is almost certain that we are losing sight of the connections between short-term events and multimillion-year climate and carbon cycle evolution.Deep time should be considered as a timescale continuum, while the models generally consider only the long-term components alone.
Over the past 30 years, models have greatly improved our understanding of carbon cycle evolution.One of the most exciting research topics here has been the simulation of the role of land vegetation on the long-term evolution of the global carbon cycle and climate.Moving from basic corrections applied to the global silicate weathering rate, newer models now include the impact of land plants not only on weathering but also on surface rugosity, albedo, evapotranspiration, and physical erosion.The most recent models also now include ecological processes adapted to deep time, simulating the behavior of extinct species.Another important breakthrough is the ability to take into account the spatial resolution of processes, allowing researchers to spatially simulate, for example, the weatherability of the continental surfaces.
From the numerical models listed in Table 1, it seems to us that the historical trend has been to refine models that were initially quite coarse.Then, depending on the object of study, the oceanic part of the model may also be refined by increasing the number of oceanic reservoirs in the box model.One of the problems of this approach is the total absence of 3D ocean dynamics, which would react directly to CO 2 variations during the simulations.This leads us to a model that we have not discussed because it follows the opposite path: the GENIE (Grid ENabled Integrated Earth system) model (Ridgwell et al. 2007).Unlike the other models discussed in this review, GENIE is a dynamic ocean model that calculates the chemical composition of the water and the dynamics of ocean transport at the same time.But running GENIE over millions of years requires a very long computer run time.Using the 36 × 36 × 16 ocean resolution, a 1 Myr model run took around two months to complete (Colbourn et al. 2015, Lord et al. 2016).This timeframe limitation is the reason why GENIE includes a rather coarse representation of the climate and a simplified description of the weathering processes.Regarding weathering, the GEOCLIM model includes the most mechanistic description of this process.But when models are refined, we rapidly lack data about the boundary conditions such as the topography of the continents in the past and the bathymetry of the seafloor.
So we have to keep in mind that some required information to constrain mechanistic models has disappeared forever.This may push modelers to experiment with a more stochastic method, generating uncertainties in the input data and boundary conditions, and performing ensembles of runs, instead of only one solution.
It is a tradition to end a review with some perspectives for future research.We anticipate exciting developments might depend on the following: 1.A better quantified CO 2 degassing flux, taking into account all of the different carbon releasing processes.2. Clearer understanding of the terrestrial vegetation.To what degree does it enhance weathering directly versus indirectly, and what were the paleo-distributions of vegetation, and their interactions with the climatic parameters, including the water cycle?3. Calculating the elemental budget of mountain ranges, accounting for all the processes that generate or sequester carbon.4. Including the behavior of marine ecosystems into geological Earth system models.

Table 1 Brief overview of the models coupling the biogeochemical cycles to climate at the geological timescale Model Reference (nonexhaustive) a Timescale (years) Steady-state carbon cycle Climate model Geochemical model type Modeled cycle(s)
a We choose to cite only one reference per model.478Goddéris• Donnadieu • Mills (Longman et al. 2022)nnualreviews.orgAccessprovidedby109.16.57.26 on 10/31/23.See copyright for approved use.concentration during times when LIPs were exposed.Smaller and high-latitude LIPs such as the ∼15 Ma Columbia River Basalt Group, however, were found to have a negligible effect on global silicate weathering [e.g., as explored in the SCION model(Longman et al. 2022)].