# Phanerozoic icehouse climates as the result of multiple solid-Earth cooling mechanisms

Andrew Merdith, Thomas Gernon, Pierre Maffre, Yannick Donnadieu, Yves Goddéris, Jack Longman, R. Dietmar Müller, Benjamin Mills

► To cite this version:

Andrew Merdith, Thomas Gernon, Pierre Maffre, Yannick Donnadieu, Yves Goddéris, et al., Phanerozoic icehouse climates as the result of multiple solid-Earth cooling mechanisms. Science Advances, 2025, 11 (7), pp.eadm9798. 10.1126/sciadv.adm9798. hal-05463327

HAL Id: hal-05463327

https://hal.science/hal-05463327v1

Submitted on 17 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

CLIMATOLOGY

# Phanerozoic icehouse climates as the result of multiple solid-Earth cooling mechanisms

Andrew S. Merdith<fnref n="1" />,²*, Thomas M. Gernon<fnref n="3" />, Pierre Maffre⁴,<fnref n="5" />, Yannick Donnadieu⁴, Yves Goddéris<fnref n="5" />, Jack Longman⁶, R. Dietmar Müller⁷, Benjamin J. W. Mills²

The Phanerozoic climate has been interrupted by two long “icehouse” intervals, including the current icehouse of the last ~34 million years. While these cool intervals correspond to lower atmospheric CO₂, it is unclear why CO₂ levels fell, with hypotheses suggesting changes in CO₂ degassing rates or modification of silicate weathering through changing continental lithology or paleogeography. Here, we construct an Earth System Model that integrates these proposed cooling mechanisms in detail. The model can reproduce the broad geologic record of ice cap expansion, allowing us to infer the primary drivers of long-term climate change. Our results indicate that recent icehouse climates required a combination of different cooling mechanisms acting simultaneously and were not driven by a single known process, potentially explaining why icehouses have been rarer than greenhouses over Earth history.

Copyright © 2025 The
Authors, some rights
reserved; exclusive
licensee American
Association for the
Advancement of
Science. No claim to
original U.S.
Government Works.
Distributed under a
Creative Commons
Attribution License 4.0
(CC BY).

INTRODUCTION

Long-term drivers of climate

Atmospheric CO₂ levels control global temperature variations over geologic time, with the concentration of CO₂ over these timescales determined by the long-term carbon cycle (1–3). Several fundamental tectonic drivers of the carbon cycle over 10– to 100-million year (Ma) timescales have been proposed over the past 40 years, which tend to be partly dependent on one another. Degassing of Earth's lithosphere adds CO₂ to the atmosphere (4–6). Weathering of silicate minerals— especially those exposed in continental arcs (5, 7, 8) or ophiolite-bearing suture zones (9)—removes CO₂ from the atmosphere and ultimately buries it as carbonate minerals (10). Weathering also liberates nutrients needed for photosynthesis, which draws down atmospheric CO₂ and sequesters carbon as organic matter (11–13).

Changes to Earth system through the Phanerozoic have funda-
mentally modified how these burial, degassing, and weathering pro-
cesses operate. The most notable example is the evolution and
colonization of Earth's surface by land plants during the Silurian-
Devonian, which irrevocably changed both the carbon cycle and
weathering processes (14–16). Since this enhancement of weather-
ing intensity and organic carbon burial by the emergence of land
plants, subsequent variations in the carbon cycle are thought to have
been driven primarily by the paleogeographic and paleotectonic
evolution of Earth's continents. The slow movement and shifting of
the continents have caused variations in the CO₂ degassing flux
(17), changed patterns of global precipitation and runoff (18), and
exposed different lithologies to enhanced weathering in the tropics
(7–9, 19, 20). The intersection of these different mechanisms is hy-
pothesized to control the overall degassing and weathering fluxes
that set surface temperatures on Earth. Many studies have sought to

reconstruct these processes over the Phanerozoic, and a succession
of “Earth System” or biogeochemical models has attempted to inte-
grate some combination of these mechanisms to reproduce the Pha-
nerozoic CO₂ or temperature record and identify the drivers of key
climatic excursions (fig. S1) (2, 21–24).

Empirical studies of the long-term carbon cycle have linked the
temporal occurrence of continental arc degassing to climate warmth
(5), their weathering to cooling (8), and the occurrence of low-
latitude ophiolite bearing sutures have also been linked to cooling
events (9). Modeling studies have suggested that mountain uplift
events could also have driven climatic cooling (25, 26). However,
there remains little agreement on which processes are most impor-
tant for causing icehouse intervals as we lack a continuous Earth
System Model (ESM) that can disentangle how all these processes
combine to drive observed changes. Continuous, global biogeo-
chemical models for the Phanerozoic, such as GEOCARB(SULF)
(2, 22) and COPSE (Carbon-Oxygen-Phosphorous-Sulfur Evolu-
tion) (21, 27) have attempted to integrate these weathering and de-
gassing processes but are nondimensional and not spatially explicit
and therefore cannot track the spatial distribution of weatherable
lithologies or hydrology through time. Furthermore, they are un-
able to spatially resolve global temperature. As a result of these
limitations, previous models cannot adequately reproduce the cy-
cles between icehouses and greenhouses evident in the geologic
record (fig. S1). Here, we use plate tectonic modeling and global
lithological maps to reconstruct the spatial and temporal variation
of solid-Earth degassing rates and weatherable lithologies from
subduction-related environments. We then incorporate these vari-
ables into a long-term 3D ESM [pySCION—a version of the SCION
(Spatial, Continuous Integration) model (23) written for the python
programming language] to assess the relative importance of the
cooling mechanisms they represent. Because of the considerable un-
certainty in the biogeography and biotic weathering ability of early
plants, alongside the unusually short duration and substantial extent
of the Hirnantian glaciation (28), we do not currently expand this
analysis into pre-Devonian times. Although a dynamic terrestrial
and marine biosphere and their organic carbon burial fluxes are in-
cluded in our analysis, we do not assess here how plant colonization
affected Earth.

<footnote n="1">¹School of Physics, Chemistry and Earth Sciences, University of Adelaide, Adelaide, SA, Australia. ²School of Earth and Environment, University of Leeds, Leeds LS3 9JT, UK.</footnote>

<footnote n="3">³School of Ocean and Earth Sciences, University of Southampton, Southampton, UK. ⁴Aix-Marseille Univ, CNRS, IRD, INRA, Coll. France, CEREGE, Aix-en-Provence, France.</footnote>

<footnote n="5">⁵Géosciences-Environnement Toulouse, CNRS-Université Paul Sabatier, Toulouse, France. ⁶School of Geography and Environmental Sciences, Northumbria University, Newcastle-upon-Tyne, UK. ⁷EarthByte Group, School of Geosciences, University of Sydney, Sydney, NSW, Australia.</footnote>

<footnote>*Corresponding author. Email: andrew.merdith@adelaide.edu.au</footnote>
# RESULTS AND DISCUSSION

## Assessing carbon sources and sinks over geologic time

Our ESM integrates a detailed reconstruction of the global CO₂ degassing rate with the global paleogeography and distribution of weatherable lithologies on the continents, overcoming the limitations of previous approaches. Degassing from the Earth's mantle and lithosphere is considered the main source of atmospheric CO₂ over million-year timescales (6). Degassing estimates have traditionally been taken from a variety of proxies (1), but recently full-plate reconstructions made using GPlates have allowed for quantitative four-dimensional (4D) approach through the mapping of dynamic processes at subduction zones over Phanerozoic time. Using this framework, a recent study coupled thermodynamic modeling to subducted carbon inventories to quantitatively assess the mobility of carbon in downgoing slabs and degassing through the Mesozoic and Cenozoic (6). Crucially, this study was able to quantify the contribution of pelagic sediments to arc-degassing over time, which increased markedly with the evolution of pelagic calcareous nannoplankton during the Cretaceous Period. Thus, our implementation of it as part of our degassing curve allows us to remove from our model the “burial depth factor,” an arbitrary scaling of degassing for pre-Cenozoic times used in all carbon-cycle models (29) that represents the increasing contribution to degassing budgets from pelagic calcareous nanoplankton since the Cretaceous.

Our analysis spans the Phanerozoic and includes degassing estimates from mid-ocean ridges, continental arcs, and continental rifts (Fig. 1A). We compile estimates of present-day fluxes from the literature (see Materials and Methods, data S1, and Fig. 1, B and C) and scale them to relevant tectonic processes to estimate how the carbon flux in different tectonic environments have changed through time (Fig. 1B). We use estimates of seafloor production and consumption taken from full-plate reconstructions to scale contributions from mid-ocean ridges and subduction zones and use the lengths of peri-continental subduction zones extracted from a full-plate model (30) to scale the amount of carbon that could be released at continental-arc environments.

We deliberately omit degassing estimates from large igneous provinces (LIPs), including siliceous LIPs (sLIPs) from our compiled degassing curve, as their individual volcanic activity (and inferred degassing contribution) typically spans only a few million years (31, 32). In contrast, our focus is the 10– to 100–milion year (Ma) variability recorded in proxy records. This does not discount the contribution to warming or perturbations to the Earth system by either LIPs or sLIPs on short (<5 Ma) timescales, where they likely cause substantial change (33, 34). Example model runs using a basic LIP degassing database (35) are shown in the fig. S14, which underscore this point. They also demonstrate how longer-lived LIPs should have lower average degassing rates, so while some may degas over intervals of >10 Ma, the modeled climate on these timescales does not show appreciable long-term changes. The degassing contribution of the Siberian Traps, which is the only LIP we include in our “default” model runs as it supplied an order of magnitude more CO₂ than any other Phanerozoic LIP (36–38), still only led to elevated surface temperatures for <5 Ma, which are shown in the temperature proxy record (Fig. 2, B and C).

Figure 1B depicts our total carbon degassing flux (shaded blue area, black outline) relative to present-day levels. It shows (i) high CO₂ degassing during the Cretaceous; (ii) modern levels during the late Paleozoic and Early Mesozoic, and (iii) below present-day levels

in the Devonian. Our uncertainty range is calculated by dividing the minimum and maximum limits from the composite curve by the 10th and 90th percentiles (respectively) at present-day (see Materials and Methods for details). In our merged curve, diffuse degassing around continental arcs and degassing from ridges and rifts are the most significant contributors to CO₂ over our analysis (Fig. 1 and fig. S2), although at present-day their contributions broadly overlap with degassing derived from subducted pelagic sediments (Fig. 1C). Tracking these data back through time underscores the importance of quantifying continental deformation (that is, mainly arcs and rifts) through time, as these regions plausibly contribute between 40 and 60% of the total solid-Earth degassing budget at any time (Fig. 1, B and C).

In the original SCION model (and in GEOCARBSULF and COPSE), Earth's surface lithology is globally uniform, and a fraction of the continental area is assumed to be volcanic (35, 39). For the current model, we incorporate a digital database of ophiolite-bearing sutures (9) and estimates of peri-continental subduction zones (inferred to represent the possible extent of continental arcs) extracted from a plate model (30). These are two distinct but overlapping lithotectonic features, with continental arcs being the surface expression of subduction-driven volcanism and ophiolitic-bearing sutures being obducted oceanic crust when two pieces of continental crust collide. For each time step, we produce a grid overlaying the paleogeography that records the fractional content of each cell that contains either arcs (primarily intermediate rocks) or ophiolite-bearing sutures (primarily ultramafic-mafic rocks) (see the Materials and Methods; Fig. 1, D to G, and fig. S4). To include important but tectonically inactive lithological and topographic features, we propagate each raster grid of active pericontinental subduction zones by one grid time step (~20 Ma), thus capturing the effects of “relict arcs” (7). Once the arc, suture, and relict-arc grids were constructed, the contributions of each to silicate weathering through time were calibrated using an arc and suture “enhancement factor” (AEF and SEF, respectively), a multiplier applied to the fractional content of the respective grids (the AEF is also applied to the relict arc grids). Weathering enhancement factors were calculated using a least-squares optimization on the weathering equations (40) at present-day and determining the enhancement factors necessary to best match the riverine silicate flux across major world rivers (41). This analysis found an AEF of 7 and SEF of 20 best reproduced the present-day data (table S1 and fig. S4).

## Phanerozoic proxy data

To evaluate our model, we consider four different types of proxy data (Fig. 2) that can be used to infer transitions between icehouse and greenhouse conditions: (i) ice cap latitude, (ii) global average surface temperature (GAST), (iii) low-latitude sea surface temperature, and (iv) pCO₂ levels. We consider ice cap extent to be the most robust indicator of icehouse-greenhouse transitions, but because no proxy data are able to provide a perfect indication of global temperature, we present our results against all four datasets. For each of these datasets, we add a line of best fit to compare our model quantitatively, although because of data coverage, some datasets are more amenable to this than others.

Our ice cap latitude estimates combine a database of glacial deposits (42) and Phanerozoic glendonites (43), reconstructed to their paleocoordinates using a full-plate model (30). Uncertainties in this dataset include the reconstructed latitudes (44) and the potential
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Sources and sinks of CO₂ implemented within pySCION. (A) Summary of different tectonic environments from which carbon is degassed. (B) Tectonic forcing curves that we use to scale each individual flux back in time divided into continental rifts, continental arcs, mid-ocean ridges, and subducting marine sediments. Our total summed curve used to drive the model is also depicted. (C) Present-day measured or modeled flux from different sources [corresponding to (A)], including from (i) continental rifts (97, 104); (ii) from subduction of marine sediments (6, 110); (iii) assimilated from lower crustal sources on the upper plate (106); (iv) diffuse degassing from the flanks of continental arcs (104); and (v) from mid-ocean ridges (6). (D to G) Spatial boundary conditions (land-sea/paleogeographic maps) (119, 120) of pySCION, depicting our lithological maps that are used to enhance the silicate weathering process.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Overview of proxy data used to constrain paleosurface conditions. Lines of best fit are shown in each panel (in blue). (A) Ice line extent estimated from preserved lithological indicators of cold climates (tillites, dropstones, and glendonites). The heatmap shows the temporal distribution in 10° latitudinal bins, considering the age uncertainty provided in the underlying dataset, and is colored by abundance. The maximum ice extent is a smoothed line (moving window of 11 Ma) requiring at least two or more deposits to be in a latitudinal bracket (explaining why in the Carboniferous this does not extend to 90°). For (C) and (D), a line of best fit and 95% confidence interval (solid blue line and interval) was calculated using a nonparametric locally weighted scatterplot smoothing (LOWESS) function (www.statsmodels.org). (B) GAST with uncertainty (blue line and envelope) in (45, 46) and GAST with uncertainty (2 sigma, indigo line and envelope) in (47). (C) Equatorial average surface temperature derived from oxygen isotopes (51). We have depicted the data as a heatmap to show the concentration. (D) pCO₂ proxy data (57, 58). Ma, million years.</figcaption></figure>
that glendonite formation might reflect the upwelling of cold water rather than cold surface temperatures (43). Figure S4 depicts both datasets separately rather than combined, as in Fig. 2A.

GAST estimates are taken from two recent compilations: first, from Scotese *et al.* (45), with uncertainty from Van der Meer *et al.* (46), who use maps of climate-sensitive lithologies through time to define paleo-Köppen belts, with shorter-term variations based on oxygen isotopes. The continuous nature of this record, and its normalization against long-term changes in seawater oxygen isotope ratios, makes it an attractive benchmark against which to compare model results. In addition, uncertainty in reconstructed latitudes, this method relies on the assumption that Köppen belts, and their indicators are broadly consistent over time (Fig. 2B). We also compare our results to the most recent “PhanDA” curve (47), which uses a data assimilation approach to combine temperature estimates from δ¹⁸O and palaeoclimate models to estimate GAST. Although this record does not correct for variations in seawater δ¹⁸O (48), there is reasonable consensus between both approaches with regard to the timings and durations of icehouses and greenhouses despite some differences in global average temperature. The most notable deviations pertain to the trajectory of late Palaeozoic cooling between the late Devonian and late Carboniferous and the magnitude of warmth during the Cretaceous and early Cenozoic, where estimates from (47) are ~5° to 10°C warmer than those of (45) (Fig. 2B). We also plot the low-latitude ±30° sea surface temperature proxy from δ¹⁸O of shelly marine fossils (49, 50). The data we use here predominantly come from a recent compilation (51), supplemented by studies focusing on more discrete time periods (Fig. 2C) (52–55). As with PhanDA, a key uncertainty with this proxy is whether predicted warm ocean temperatures in the Early Paleozoic are realistic or whether they reflect changing oxygen isotope composition of seawater (3, 50, 56).

A range of proxies are used to constrain past CO₂ levels based on either photosynthesizing organisms (e.g., stomatal density in leaves, δ¹³C in liverworts, and alkenones) or on seawater pH through boron isotopes in foraminifera and pedogenic carbonates. We use three recent compilations (57–59) that provide proxy data back to the Ordovician. Uncertainty ranges for these proxies are based on different methodological steps, particularly the estimation of soil CO₂ levels for the paleosol record (60) and the boron isotopic composition of seawater for reconstructing paleo-pH (61).

## Modeled Phanerozoic climate and driving mechanisms

We now run the SCION model through the Phanerozoic. CO₂ degasing follows our compiled curve in Fig. 1B, and weatherable lithologies are mapped onto the model land surface (Fig. 1, D to G). We evaluate our model by comparing its prediction of the geographic extent of ice sheets against the proxies from Fig. 2A, with a secondary focus on surface temperatures and CO₂ levels in Fig. 1 (B to D), which are less certain. SCION assumes that an ice sheet is present if the grid cell continental temperature of the model output has a mean annual temperature below -10°C (e.g., 62); as a sensitivity test, we also plot the predicted ice curve for -12°C and -8°C. Figure 3 shows our “default” model ensemble, which includes all processes and considers the uncertainty bounds in the CO₂ degassing rate and uncertainty in the background weatherability of the non-arc and nonsuture continental grid cells through time. We also show the calculated Wasserstein Distance—the mean distance between two distributions—between the uncertainty range of our results and that of the proxies in each panel.

We find that the model (Fig. 3) can replicate several key features of the Phanerozoic ice cap paleolatitude record (42), showing cooling from a warm early Devonian to a Carboniferous-Permian icehouse, followed by warming into a Triassic greenhouse, cooling through the Jurassic cool interval, into a Cretaceous hothouse and lastly into a Cenozoic icehouse (63–65). The model produces more mixed fits to the CO₂ proxy record, compiled global average temperature record (45, 47), and low-latitude oxygen isotope record, although the uncertainty in these records is much higher. Key misfits in the combined records, identified by the model-proxy distance metrics shown at the bottom of each panel, occur during (i) the Devonian—where ice cap expansion and CO₂ drop are roughly matched, but temperature proxies are not (max ~10 °C distance); (ii) the Permian, where temperature proxies are well-matched, CO₂ proxies are highly variable, but the model ice cap extent is located at higher latitudes than the associated proxy record (max ~50° distance); (iii) the Early Cretaceous, where the model has temperature within the wide uncertainty in different reconstructions but has higher CO₂ [and max ~3000 parts per million (ppm) CO₂ distance] and has more restricted ice caps than the geologic record suggests; and (iv) the Eocene, where the model is colder and has lower CO₂ than the proxy records suggest. The estimation of ice sheet extent in our model is likely complicated by the underlying paleogeographic maps. For example, the mismatch in ice expanse during the Permian is possibly because the underlying paleogeographic-topographic model is not conducive to producing low temperatures at low latitudes, where, at this time, localized and regional topography may be more important. We note that our model predicts the iceline extending to ~60° latitude, which is reasonably consistent with the bulk of the Gondwanan ice sheet (66). Our results over the Early Triassic also differ from the GEOCLIM model (67) because we add a weathering dependency on erosion and also have abundant low- to mid-latitude ophiolitic sutures at this time that appear to help drive some component of Early Triassic cooling (fig. S6).

Overall, despite some apparent misfits, our model makes a generally reasonable prediction of long-term greenhouse-icehouse trends at the 10- to 100-Ma scale within the uncertainty offered by multiple different proxy datasets. It performs better [with a mean Wasserstein distance to modeled GAST (45) of ~3.3] than the original SCION model (Wasserstein distance of 4.07) and previous simpler and more idealized model frameworks (e.g., GEOCARBSULF and COPSE, Wasserstein distances of 3.41 and 4.57 respectively; fig. S1). We now seek to better understand what is driving these model predictions by examining how individual processes act to warm or cool the model climate. To do this, we systematically rerun the model with individual drivers turned off to quantify the influence of each individual driver on changing the modeled GAST curve relative to the default GAST curve of (Fig. 3D) (45)—we choose this proxy curve as a reference for this section because (i) it is continuous in high resolution over the Phanerozoic, (ii) it does not reach a maximum beyond which it can no longer vary (e.g., 90° for the ice line), and (iii) it is robust to changing seawater oxygen isotope compositions. We assess the impact of the climate drivers by calculating the ratio between the Wasserstein Distance of the default model distribution and each individual run with a driver switched off, with the temperature distribution at each time step. This allows us to test the influence of several previously suggested drivers of Phanerozoic Icehouses (Fig. 3E) (5, 14, 25).

In our model, the Devonian-Carboniferous cooling is predominantly driven by changes in paleogeography and a decrease in
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** **pySCION model results for the Phanerozoic.** Results from our default run are plotted in (A) to (D), where gray lines show 10% of the ensemble (1000 runs per ensemble), and the solid amber lines are the mean ± 1 SD. In some cases where ensemble runs overlap, the lines are black [e.g., (A)]. For comparison, proxy data are in red-brown, and their lines of best fit and confidence are shown in blue (see Fig. 2 caption for details). Solid black lines show the distance between proxies and model. (A) Ice line results. The default model run assumes a critical temperature of -10 °C, and the mean values of alternative runs assume critical temperatures of -8° and -12°C that are represented by the coarse and fine dashed amber lines, respectively. (B) Atmospheric CO₂ results. (C) Modeled equatorial (±30°) average surface temperature against equatorial (±30°) sea surface temperature from δ¹⁸O (51). (D) GAST and uncertainty of Scotese *et al.* (45) in blue and Judd *et al.* (47) in indigo. (E) Relative influence of each climate driver when it is turned “off” (no change relative to present-day)—a measure of how influential each driver is through time. Temperature, ice line, and CO₂ for each series of runs are provided in figs. S9 to S14. Ma, million years.</figcaption></figure>
degassing, although arc weathering also plays a prominent role (Fig. 3E). Both changes in palaeogeography and degassing have been individually proposed to initiate glaciation (5, 14, 25). The former through a reduction in sea level (68, 69) and emergence of continental shelves and islands (70–72) and the latter by a reduction in degassing that has been linked to a slowdown in subduction rates as inferred by a lull in the magmatic zircon record (5) triggering cooling. However, our model suggests that it is only through a combination of these processes acting together that an extensive icehouse could be initiated (figs. S5 to S8). Our results suggest that once initiated, the prolonged icehouse was maintained predominantly by the paleogeographic configuration of Earth, wherein the low-latitude Variscan-Hercynian Orogeny (which formed in the Permian from the collisions between Gondwana and Laurussia to form Pangea) increased physical weathering, maintaining cool temperatures (25).

The warm Mesozoic climate state in our model is also most sensitive to changes in degassing and paleogeography, specifically the increase in volcanic degassing driven by the initial rifting stages of Pangea breakup between 245 and 200 Ma. Here, we recognize an interplay between both mechanisms, whereby degassing of continental rifts increases the build-up of $pCO_2$. However, this degassing occurs without forming new ocean basins, thus, in a palaeogeographic sense, preventing an increase in supply of nutrients derived from chemical weathering (fig. S6) reaching young ocean basins, resulting in warming (although our model is not yet able to capture the full extent of this as weathering products are always assumed to reach the ocean). Arc systems and ophiolitic sutures still contribute to enhanced silicate weathering (8); however, their contribution does not change substantially over the Mesozoic as their spatial footprint remains mostly constant, and their exposure is often balanced by hydrological changes (i.e., less rainfall and runoff). Our modeled warmer Mesozoic climate is punctuated by ephemeral glaciation in the Jurassic, in reasonable, but not perfect, agreement with the geologic record of glacial deposits at this time (42). Our analysis suggests that this short-term decrease in global temperature (evident in all three temperature proxy curves) and increase in high-latitude glacial deposits is mostly driven by a reduction in degassing from both rifts and mid-ocean ridges (Fig. 1) and an increase in erosion and weathering following the breakup of Pangea (figs. S6 and S7).

Last, cooling during the Late Cretaceous and Cenozoic is most sensitive to our modeled decrease in solid-Earth degassing (Fig. 2), alongside changes to paleogeography and arc weathering. In our compiled degassing curve, the decrease is principally driven by the reduction in degassing around continental arcs, as the decrease in mid-ocean ridge degassing is offset by the increase in degassing driven by the subduction of pelagic sediments over this time (6). Tighter quantifications of the contributions of diffuse degassing around arcs and rifts would potentially change these sensitivities. Changes in paleogeography, such as those resulting from the closure of the low-latitude Tethyan Ocean and formation of the Himalayas as well as exposure of arcs and ophiolitic sutures around Southeast Asia, also played a role in increasing the efficacy of the silicate weathering cycle (fig. S8). While previous studies have proposed each of these drivers acting in isolation to explain this climate shift (5, 19, 26), our model results suggest that all three played important roles in contributing to global cooling.

## No single driver for Phanerozoic icehouses

To further explore whether any of the key proposed processes could exert a first-order control on Phanerozoic climate, we rerun the

model subject to a single climate driver (Fig. 4 and figs. S9 to S13). As a baseline reference, we ran the model with all drivers turned off. As expected, this produces a predominantly icehouse Earth similar to that of the late Cenozoic (black dashed lines in Fig. 3 and fig. S7). We find that when we run the model systematically with only each discrete driver, the model produces an overall poor fit for the Phanerozoic icehouse-greenhouse record. While the weathering of volcanic arcs is important, it does not change substantially over the Phanerozoic and, therefore, is unable to independently force major climate change (Fig. 4 and fig. S15). Weathering of ophiolite-bearing sutures changes substantially but constitutes only a minor component of the total weathering budget (Fig. 4 and fig. S15).

Degassing and paleogeography show a greater ability to drive climatic change when considered alone. While some of these model runs yield reasonable fits for certain time periods, we observe numerous discrepancies with the climate record when each driver is considered in isolation (Fig. 4). Notably, these discrepancies are far greater than those observed in the “combined drivers” model (Fig. 3). In particular, the “degassing only” run can reproduce a Cretaceous hothouse, but it is unable to reproduce a warmer Triassic (Fig. 4). No single mechanism can reproduce a warm Devonian, with all four single-driver runs producing an ice-dominated Devonian world (Fig. 4A). This strongly suggests that over the Phanerozoic, the combination and intersection of multiple different drivers were required to drive substantial variations in climate on 10- to 100-Ma timescales. While these results do not negate the correlations between individual events and climatic excursions, for example, tropical exposure of sutures with global cooling (9) or the correspondence of continental arc degassing to the Mesozoic greenhouse (5), it does suggest that single processes in isolation are incapable of explaining Earth’s long-term climate state.

## Model limitations

Despite some of the advances that our model can provide, there remain several key limitations that should be addressed in future work. The model uses a single paleogeographic reconstruction and a single series of paleoclimate model runs to inform the climate emulator, meaning the uncertainty in these variables cannot be assessed. Moreover, the model grid times (i.e., times at which paleogeographic maps are used to drive the paleoclimate models) are quite coarse, on average, spanning ~25 Ma, but at times the gap is 40 to 50 Ma. Two consequences of this approach are, first, that many key paleogeographic transitions (such as smaller terranes and blocks entering/leaving tropical zones, emergence of continental shelves, and formation of islands arcs) are missed in our modeling strategy. Machine learning has been used to try to fill these gaps but is not capable of reproducing hydrology to a high standard (73). Second, we are unable to easily test alternative paleogeographic configurations, such as Pangea-A versus Pangea B [although see discussion in (74) about a resolution to this controversy] or differing locations of terranes in the Permian Meso-Tethyan ocean. A strength of the (py)SCION framework moving forward is that as additional paleoclimate models are generated, they can be incorporated within the climate emulator in the model—irrespective of the spatial resolution—to better resolve these fundamental questions.

It is well established that the FOAM (fast ocean-atmosphere model) general circulation model has a low climate sensitivity (75) and, therefore, requires more $CO_2$ to increase temperature than more modern climate models. In addition, FOAM has only been
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Results of the model where only one driver is turned on. Proxy data in each panel is the same as those shown in Fig. 3. The default run (all drivers on) is plotted as black lines with gray-shaded area. The figure shows (A) ice line, (B) pCO₂, (C) equatorial average surface temperature, and (D) GAST of Scotese *et al.* (pale red) (45) and Judd *et al.* (47) (indigo). “Arc” and “suture” only refer to the modeled contribution to silicate weathering of continental arcs and ophiolitic bearing sutures, respectively. Note that organic carbon burial is still active in these runs. Ma, million years.</figcaption></figure>

run for changes in CO₂, ignoring other greenhouse gases. This may be why, despite our model predicting temperatures often within the expected ranges (i.e., 15° to 25°C), our pCO₂ estimates are sometimes far higher than the proxy data suggest. The FOAM model also has quite a coarse spatial resolution (7.5° × 4.5°, corresponding to 48 longitudinal cells by 40 latitudinal cells), which is not high

enough to properly differentiate finer paleogeographic features such as the opening or closing of oceanic passages. It also struggles to properly capture small-scale drainage systems, such as those that might occur on the seaward side of a continental arc. Consequently, while we consider our model capable of capturing 10- to 100-Ma trends, it is unable to resolve trends on a finer spatial or
temporal scale, even if they might be driven by the processes that we model.

A notable biogeochemical limitation of the model is the simple treatment of the organic carbon cycle. Despite implementing a spatially resolved weathering scheme for lithotectonic features and the silicate weathering process, the model implementation of oxidative weathering of sedimentary organic carbon is only related to local runoff, whereas erosion rate and lithology are likely also important (11, 76). Marine organic carbon burial is not calculated spatially despite the geologic record suggesting it is heterogeneous (77), and the model also does not feature a spatially resolved dynamic vegetation model for calculating land-derived organic carbon burial. Such a model has been developed for the Mesozoic and Cenozoic and has been shown to alter SCION model predictions (78), but it has not yet been developed for the Paleozoic floras. Overall, we consider that a more sophisticated paleogeography, climate emulator, and organic carbon cycle would improve the model and could help further reconcile model results with proxy data during the misfit periods that we note. However, we believe that these advancements are unlikely to reverse the major cooling or warming trends in our results, as the net effect of amplified continental weathering is still likely to be global cooling. Although organic carbon weathering is a CO₂ source, it is accompanied by nutrient delivery and a reduction in atmospheric O₂, which both encourage additional organic carbon burial.

We conclude that the prolonged Phanerozoic icehouses during the Carboniferous-Permian and Cenozoic were initiated and maintained through the interaction of multiple “climate drivers” rather than being the result of a single dominant process. This helps to explain why icehouse periods have been relatively rare during Earth history, why so many different drivers have been proposed, and why icehouse-greenhouse transitions have previously proven difficult to reproduce in more simplified carbon cycle models (3). Our results suggest that changes in paleogeography, which result in spatial-temporal variations in the hydrological cycle, alongside changes in degassing, are the two most important mechanisms. Continental arc weathering intensity is substantial during warm periods, in agreement with (8), but is not able to cause an icehouse transition by itself. Our findings suggest that glacial climates on Earth depend on the interplay between multiple global processes that occur simultaneously, carrying important implications for understanding Earth’s long-term habitability. The long-term temperature of our planet appears to be mechanistically regulated toward an ice-free state rather than a glacial climate.

# MATERIALS AND METHODS

## SCION and pySCION

The model runs were performed using pySCION, a pythonic version of the SCION model (23). The underlying SCION model (written in MATLAB) is available at https://github.com/bjwmills/SCION. Then, porting the model to python, the framework is extended to incorporate the impact of different lithologies on the weathering equations. The pySCION model functions as an amalgamation of the COPSE biogeochemical forward model (27) and the spatially resolved GEOCLIM model (24). Here, the pySCION model uses the data structure of GEOCLIM to resolve key hydrological and climate fields (e.g., runoff and temperature) at each grid time step to calculate erosion and weathering. However, unlike the GEOCLIM model, which is non-continuous, in between grid time steps, the COPSE biogeochemical

model is used to integrate fluxes between time steps. Thus, pySCION can leverage spatially resolved data through time, allowing us to interrogate the effects of long-lived geologic features on the carbon cycle. The model is available as part of the submission in data S2. Full model equations for SCION and pySCION are provided in (23). They are also available at the original SCION GitHub repository.

## Arc and suture weathering grids

The same approach was taken to include both continental arcs and ophiolitic-bearing sutures in the model, whereby we calculate the fractional grid area of each cell that was composed of each lithotectonic feature that represent estimates of the spatial extent of intermediate and mafic-ultramafic rocks, respectively. For continental arcs, we used the plate model of Merdith et al. (30) [the Phanerozoic portion of this model was adopted from models in (79–81)] to extract all subduction zones within a specific distance of a polygon (representing continental crust) at the pySCION grid times. The distances used were 250 and 436 km, which represent the mean and mean plus one SD distance of arc volcanoes from trenches at the present-day (82). We refer to these arcs as “pericontinental” arcs rather than strictly continental arcs because it is possible that our analysis includes some arc systems that may not be geologically classified as continental arcs. For sutures, we used the database of Macdonald et al. (9), who assembled a record of ophiolite-bearing sutures arising from arc-continent collisions during the Phanerozoic. At each time step, the maps of arcs and sutures were rasterized onto a regularly spaced rectangular grid using the Phanerozoic reconstruction model provided in Macdonald et al. (9)—a modified version of (72, 83). A particular shortcoming with this approach is translating linear features into area representative domains (as the digitization provided is in the form of polylines). The resolution of the rasterized grid for sutures was 18.75 • 18.75 km, reflecting the scale of the maps they were taken from (1:3,750,000) (84), and is slightly finer than estimates of ophiolite thickness from other compilations, where ophiolite sheet width typically ranges between 40 and 60 km (85, 86). For the arc data, we used a coarser grid cell size of 100 km, equivalent to the order of magnitude of width of continental arcs from global compilations and analyses (7, 87). In addition to rasterizing the data to these fine(r) grids, we also rasterize it to the resolution of the pySCION model (48 • 40) (e.g., fig. S3). Thus, at this point of our analysis, we have two grids for each of arcs and sutures (four in total), one representing the resolution of the data and the second at the resolution of the pySCION model.

To determine the fraction of each pySCION grid cell, we use the latitude-longitude bounds of the coarser pySCION resolution map (fig. S3C) to cluster the finer resolution map. For each grid cell at the coarse resolution, we sum the area of corresponding arc or sutures from the finer resolution map to calculate a fraction of the total grid cell area of the coarser grid that is either arc or suture. Because of differences between the underlying paleogeography and the reconstruction used, coupled with the absence of a rotation file of the underlying SCION paleogeography, some manual manipulation was required to better align the reconstructed sutures and arcs with their position on the land-sea masks.

In cases where sutures overlie arcs or relict arcs, we do not subtract the area contribution attributed to sutures from both arcs and relict arcs, on the basis that sutures are features that are emplaced tectonically onto existing arc systems. Conversely, in situations where arcs overlie relict arcs, we subtract the area of active arcs from
relict arcs to prevent double counting these regions. As the resolution of our digitization of ophiolite-bearing sutures is 18.75 km, we multiply their abundance by 3 to better reflect their preserved width [e.g., see compilation by Porkoláb et al. (86)]. The digitized width of arcs is roughly equivalent to their actual width, so we do not use a multiplier for their abundance.

## Calibrating weathering

As the standard weathering in SCION assumes an average lithology and reactivity across the entire globe, we introduce a weathering factor to represent the enhanced weathering that is expected from silicate rich rocks. We do this by optimizing the parameters in the SCION weathering equations (40, 88) in the modeled 0-Ma timeslice against the present-day measured weathering flux from a series of drainage basins (40) using a least-squares optimization from the python library SciPy. We include in our modeled weathering flux the distribution of pericontinental arcs and ophiolitic sutures and introduce an AEF and SEF to the weathering equations to isolate what enhancement is necessary to match present-day weathering rates, as measured from dissolved loads in world rivers. We use the same AEF value for our enhancement of relict arcs. The values we let vary, their range, and the results of the optimization are summarized in table S1. Bounds for K, K<sub>w</sub>, and σ are taken from Maffre et al. (40).

For the AEF and SEF, we enforce a lower bound of 7 and an upper bound of 20. These enhancements were selected to reflect present-day enhanced weatherability of silicate rocks relative to an average weatherability of continental crust. The lower limit of 7 was established as minimum increased reactivity (35). We established an upper limit based on the area of peri-continental arcs and ophiolitic-bearing sutures in our model at present-day, and their total contribution to the silicate weathering cycle. It has been suggested that volcanic silicate weathering accounts for 30% of total weathering at present-day (22) on the basis of 87Sr/86Sr and 187Os/188Os curves. In our analysis, the combined present-day area of pericontinental arcs and ophiolitic bearing sutures is 1.5% of Earth's surface area. Thus, we infer that a 20-fold weathering enhancement is a reasonable upper limit.

To assess our optimized values for arc and sutures factors (in addition to the other parameters), we plot them against measured data from present-day (40) (fig. S3, E to G). We find that our optimized enhancement factors produce a better fit to present-day riverine flux (fig. S3G) than the unoptimized run with no enhanced weathering (fig. S3F). In particular, there are no outliers that overestimate weathering, indicating that our results (and weathering enhancement factors) are likely to be conservative.

To implement these grids within the pySCION framework, we multiply the fractional content of each land cell that is “ARC,” “RELICT ARC,” or “SUTURE” according to the appropriate weathering multiplier (that is, AEF and SEF) using the following equations

$$ \omega_{arc} = \omega_{silw} \times [1 + K_{arcfrac} \times (AEF - 1)] \quad (1) $$

$$ \omega_{relictarc} = \omega_{silw} \times [1 + K_{relictarcfrac} \times (AEF - 1)] \quad (2) $$

$$ \omega_{suture} = \omega_{silw} \times [1 + K_{suturefrac} \times (SEF - 1)] \quad (3) $$

These equations are then combined as

$$ \omega_{\text{totalsilw}} = \omega_{\text{silw}} + \omega_{\text{arc}} + \omega_{\text{relictarc}} + \omega_{\text{suture}} \quad (4) $$

## Degassing

Previous models, including COPSE (27) and GEOCARB(SULF) (89), invert Phanerozoic sea level to approximate crustal production at mid-ocean ridges and then use that measure as a proxy for degassing (89, 90) (i.e., CO₂ output from the lithosphere into the atmosphere). To quantify degassing in our model, we divide it into five different tectonically driven sources (Fig. 1A): (i) continental rifts; (ii to iv) convergent margins, including carbon sourced from the (ii) subduction of carbonaceous ocean sediments, (iii) degassed from the magma chambers (either actively or passively), or (iv) degassed diffusively on the flanks of volcanic arcs; and (v) mid-ocean ridges. Our analysis assumes that the partitioning of carbon from solid earth processes within these environments is linear through time—the same proportions of carbon are being emitted now as in the past, and the only factor changing this amount is variation in tectonic fluxes. Our approach is therefore to treat each source of degassing separately, determining first an appropriate present-day flux, then deriving a scaling factor able to scale the present-day flux back in time, and, lastly, combining the scaled measurements into a single “total” degassing curve. We use the 10th and 90th percentiles of the present-day flux to divide the minimum and maximum estimates to express the degassing curve as a flux relative to present-day. When we fully preserved the present-day uncertainty (i.e., the present-day ratio of the compiled curve is the same as the present-day relative curve), our uncertainty envelope encompassed all temperatures between 15° and 25°C over the Phanerozoic.

Our final flux is assembled from the raw data before being smoothed using a Savitzky-Golay filter with a 21-Ma moving window and polynomial order 1. Data S1 provides the data we used. We selected a 21-Ma moving window as this provided a good balance to smoothing short term perturbations that could be artifacts in the underlying plate- tectonic model while still preserving variations on a 10-Ma scale.

## Scaling mid-ocean ridge processes

To scale the processes occurring at mid-ocean ridges, we first derive a tectonic forcing curve constrained by three (I to III) separate curves estimating a gross tectonic flux at ridges and/or subduction zones. We use the range covered by all three curves to accommodate uncertainty in gross estimates of pre-Pangea subduction and spreading flux. (I) A recent study (17) calculated a tectonic degassing curve back to 550 Ma. They used a hybrid approach of crustal consumption curves calculated from full-plate reconstructions (91) for the Mesozoic-Cenozoic and zircon data for the Paleozoic (80, 92) (where the zircon flux is proposed as a proxy for subduction flux). For our mid-ocean ridge, we use their hybrid curve (derived from zircons and full-plate models). Our second and third estimates of tectonic flux use crustal consumption (II) and production (III) curves from a plate model (30). We assume that over 1- to 10-Ma timescales crustal production (i.e., area of new seafloor created at mid-ocean ridges) is equivalent to crustal consumption (i.e., area of seafloor consumed at subduction zones), which can be demonstrated quite for the past 150 Ma (6), but before this time, there are differences between the two curves. We extract both crustal production and consumption curves from 550 to 0 Ma. To create our final tectonic forcing curve that encompasses all three individual curves, we take the group minimum and maximum across all three curves at each time step.
## Scaling continental arc processes

To scale the processes occurring at continental arcs, we adopt a similar approach to our method for mid-ocean ridges. We constrain continental arc flux by three separate curves (I to III), one derived from zircon spectra (17) and the second and third from the subduction flux extracted only from pericontinental subduction zones using two different distance thresholds. We deliberately refer to these arcs as pericontinental, rather than continental arcs, as limitations within plate models (e.g., rigid polygons and uncertainty about extent of continental crust back in time) make it an imperfect estimate of “true” continental arc length. As with the mid-ocean ridge flux, we use the range covered by all three curves to accommodate uncertainty in gross estimates of continental arc flux. (I) The study of Marcilly et al. (17) also produced a degassing estimate derived completely from zircon. As their curve is explicitly defined from arc-derived zircon, we take this as a direct proxy for continental arc flux through time. Our second and third estimates of tectonic flux use crustal consumption curves from the plate model (30). However, this time we filter the subduction zones in the model by their distance to the nearest polygon in the plate model to extract those that could plausibly be a continental arc (rather than an oceanic arc), in the same manner we applied to produce the continental arc silicate weathering grids. We used two distance thresholds (82), representing the (II) mean distance (250 km) and (III) mean distance plus one SD (436 km) from a modeled subduction zone to volcanic arc at present day. We then multiplied these lengths by the time-sensitive, average convergence rate extracted from the plate model (30) to determine an estimate of minimum and maximum pericontinental arc flux. To create our best estimate of continental arc flux curve that encompasses all three individual curves, we take the group minimum and maximum across three at for each time step.

## Rifts

CO₂ can be degassed diffusely from active continental rifts by mobilization of subcrustal reservoirs and migration via extensive fault and fracture networks (93, 94). There is a large range in degassing measurements across rift environments (95), with fluxes from the Madagi-Natron Basin in East African Rift (EAR) system thought to possibly contribute up to ~1 Mt C/a (93), while the back-arc Taupo Volcanic Zone in New Zealand (which occupies a similar area to the Madagi-Natron Basin) contributes an order of magnitude less, ~0.1 Mt C/a (96). Propagating fluxes along the entire EAR lead Lee et al. (93) to propose a total rift degassing of 14 to 26 Mt C/a, roughly equivalent to degassing from the entire mid-ocean ridge system, although this range was re-evaluated by Hunt et al. (94) to a more conservative 1 to 8 Mt C/a. Werner et al. (95) suggested fluxes from other rifts could potentially add another 8 to 11 Mt C/a, for a total global flux of 13 to 16 Mt C/a [using the mean EAR value in (94)], although cautioned that many more measurements are needed. Last, Brune et al. (97) compiled a database of rift lengths for the Mesozoic-Cenozoic and estimated global fluxes of CO₂ from rifts at present-day. They proposed two estimates of total contribution to C from present-day rifts, a mean estimate of 54 Mt of C/a, and a conservative estimate of 13.5 Mt of C/a. The primary difference of the two measurements reflects uncertainty around the large fluxes coming from parts of the EAR relative to other rifts. As the conservative estimate of Brune et al. (97) falls reasonably well in the center of the tentatively proposed global range of Werner et al. (95), we use their range as a conservative uncertainty estimate of rift degassing.

To scale this present-day value back through time, we use a rift-length database (98), as organized and processed in Brune et al. (97) (450 to 0 Ma) and extended by Merdith et al. (99) (450 to 1000 Ma). However, one issue with the transition between the two databases is that the Neoproterozoic and Early Paleozoic database predominantly identifies rifts that form during continental breakup and/or terrestrial detachment and migration. Comparatively, the database in (98) contains not only these rifts but also rifts formed in convergent settings, such as back-arc basins. This is likely a function of preservation bias, and to counteract this issue, we calculated the fraction of rifts formed in convergent settings (98), which for the period 450 to 0 Ma is 0.44 [i.e., on average, 44% of total rift length in the database (98) formed in convergent settings and is unaccounted for in the estimate in (99)]. Assuming that this proportion is constant back to 600 Ma, we applied a rift length multiplier

$$length_{adjusted} = length \times (1 - 0.56)^{-1}$$

to the rift lengths in (99).

Last, the rift database (98) only provides start times of each rift in the form of a geologic time period or epoch, and so an implicit assumption of our quantification is that the rifts were only active for the duration bounded by their identified geologic timescale age. This creates an issue for very young rifts, initiated since the Neogene, as their duration age (i.e., time we consider them “active” and degassing) becomes very short relative to older rifts that might have been attributed a start age of a longer epoch. Our estimate of rift lengths drops quite abruptly from the Pliocene to present-day (~25,000 to ~8000 km) because rifts given ages of “Pliocene” or “Pleistocene,” for example, do not contribute to the present-day estimate. Therefore, when we smooth our degassing curve, we do not fix the present-day estimate, and the result is that the smoothed curve increases at present-day and reduces the influence of rifts back through time.

## Peri-continental arc derived carbon

Much work was been done on disentangling and isolating different sources to solid-Earth degassing from continental arcs (4, 100, 101). Typically, the breakdown of continental arc carbon degassing is done on the basis of temperature, separating deeply sourced magmatic carbon (>450°C) from more shallow hydrothermally derived carbon (102–104). Here, we separate carbon sourced from continental arcs into two different categories: magmatic carbon sourced deeply from either the (ii) subducting sediments or (iii) continental crust and (iv) diffuse degassing along continental arcs (such as through soils and springs).

## Magmatic degassing

Within subduction zones, slab-derived carbon is liberated due to the warming and devolatilization of the slab and/or sediments being subducted and liberation of carbon stored in lower crustal reservoirs (105). Carbon pathways in subduction zones are poorly constrained, and most of the estimates of bulk carbon fluxes from these sources are based on mass-balance estimates between carbon entering a subduction zone and what is degassed (4, 100). In addition to improvements in observation of volcanic degassing over the past 20 years [e.g. (104)], advancements have also involved identifying isotopic signatures of different sources (106, 107). Analyzing data from a global compilation of fluxes, Aiuppa et al. (106) isolated high temperature
volcanic degassing. They argue that based on the S/C and Ba/La ratios, and LILE enrichment patterns of volcanic arcs, that most of the carbon emitted from these volcanoes come from the subducting slab (68 to 81%), particularly, from the sediments of the subducting slab. The remaining 19 to 31% of magmatic global volcanic arc emissions are thought to be sourced from assimilation of carbon reservoirs in the lower crustal column [e.g., arc magma-limestone interaction; (105)] based on their heavy enrichment in LILEs.

Volk (108) first speculated that an increase in carbon degassing through the Cretaceous present-day must have been driven by an increase in deep-ocean carbon burial [see also 109]. This hypothesis has been part of all biogeochemical models since the original GEOCARB and COPSE models (21, 29), typically referred to as the “B-factor” (for carbon burial), and scale down degassing before 140 Ma to 75% of present-day. Müller et al. (6) provided thermodynamic modeling of carbon pathways to estimate the volume of carbon sourced from subducting slabs over the Mesozoic and Cenozoic and argued that most of slab-derived carbon comes from subducting sediments (total influx into overriding plate at present-day of 31 to 70 Mt C/a). This study, which focused on the oceanic conveyer belt component of the deep carbon cycle, did not consider overriding crustal sources of carbon, and the atmospheric outflux (i.e., a mass balance approach between inputs into the overriding plate and outputs measured along volcanic arcs) was calibrated by Bekaert et al. (110), who estimated 14 to 24 Mt C/a based on carbon/sulfur (C/S) ratios from (103). Müller et al. (6) ascribed the entire reservoir from subducted sediments to this amount being outgassed; however, here, as we also consider assimilated CO₂ from overriding carbon storage, we partition this total amount into two sources, 68 to 81% from the sediment estimation of Müller et al. (6) and the remainder from the overriding plate. In effect, we take the exact outgassing flux of Muller et al. from 250 to 0 Ma and scale it down to 68 to 81% [i.e., (ii) 10 to 19 Mt C/a at present-day]. This represents the component of carbon sourced from pelagic sediments and replaces the B-factor in our model. Second, using the present-day estimate of 14 to 24 Mt C/a, we take the remainder 19 to 32% as the present-day flux of carbonate assimilation in the overriding plate [(iii) 2 to 7 Mt C/a]. This later range is scaled back in time through our pericontinental arc flux curve compiled above.

## Diffuse or passive degassing along continental arcs

The contribution of diffuse degassing along continental arcs—such as CO₂ released through soils, spring lakes, and vents on the flanks of volcanoes—is thought to be a large, but uncertain, contribution to tectonic degassing. We use the present-day estimates of a recent compilation (104), which suggest 14 ± 1.5 Mt C/a is released. This estimate only pertains to diffuse, low-temperature degassing on continents and does not include contributions from other submarine sources [e.g., (111)]. We use our pericontinental arc flux to scale these values back in time (fig. S5).

## Mid-ocean ridges

Mid-ocean ridge degassing is driven by the release of CO₂, and other volatiles, from the mantle beneath a spreading ridge. In these environments, the amount of CO₂ released during the formation of newly formed oceanic crust is typically calculated by analyzing the ratio of CO₂ to another element, chiefly ³He or Ba (112). Most estimates using these methods place the present-day flux at mid-ocean ridges between 10 and 100 Mt C year⁻¹ [e.g., (113)], corresponding

to a mantle concentration of around 100 ppm. Keller et al. (114) used a geodynamic model to calculate the spreading rate-dependent degassing of mid-ocean ridges, given a mantle concentration of C. This was used to estimate (6) carbon degassing at mid-ocean ridges over the past 250 Ma, with a resultant present-day flux of 13 ± 4 Mt C year⁻¹, using a mantle carbon concentration of 100 ± 20 ppm [e.g., (115)]. In our analysis, we use the present-day derived estimate of Müller et al. (6) and compile a composite tectonic-forcing curve back to 550 Ma to scale the relative MOR degassing back in time. This is a simplification of the method of Müller et al. (6) who calculated the flux dynamically based on known spreading rates at known mid-ocean ridges (114). However, the bulk of our analysis covers the time interval pre-Pangea breakup, thus making it impossible to use the same method unilaterally with the same degree of confidence. Instead, we scaled the present-day estimate of Müller et al. (6) using our tectonic forcing curve.

## Degassing of the Siberian Traps and other (s)LIPs

To investigate how the rapid release of large amounts of lithospheric carbon via large igneous provinces (LIPs) may affect our model, we performed a model run including estimates of LIP degassing (a table this data is included in the supplementary model files). Our model run is a nonexhaustive compilation of Phanerozoic LIPs (35) but does include key events including the Siberian Traps, Central Atlantic Magmatic Province, Ferrar-Karoo, and Deccan Traps. We also include the Whitsunday Volcanic Province, a ca. ~15-Ma long-lived sLIP that formed in the Cretaceous on the eastern margin of Australia (116). The results of this model run are shown in fig. S14. We found that LIPs do not perturb the model climate on >10 Ma timeframes— which are the focus of this study. This is because LIPs are either shorter lived than 10 Ma, or if they do have longer activities (e.g., Whitsunday Province), then this leads to a lower rate of CO₂ injection. Some LIPs do perturb the model climate on <10 Ma timescales, and the Siberian Traps causes a multidegree warming over several million years in our model. This is because carbon release rate from this LIP is thought to be up to an order of magnitude higher than other LIPs, up to 1 × 10¹³ mol C/a (38) over ~2 Ma. For this reason, we include that the Siberian Traps in our default model runs as an endmember (maximal) example. Our analysis does not preclude LIPs acting as important climate drivers on shorter timescales (e.g., through pulsed degassing), but these events are not the focus of this work.

The implementation of these LIPs into pySCION is done in a similar manner to previous work (117), where we include LIP degassing as a forcing that is described by a Gaussian function, with the midpoint defined halfway through the LIPs duration and the width set to half the duration time. The calculated pCO₂ contribution from this Gaussian is then added to the degassing flux. We assume this carbon has the same isotopic composition (−5‰) and model the carbon mass balance as such.

## Isolating mechanisms

To test the influence and strength of individual mechanisms, we isolate their contributions to be consistent with what is observed at present-day.

## Sutures, arcs, and degassing

For degassing, we assign a value of ‘1’ for each time step, thus ensuring that degassing remains constant through time at the present-day
value. For sutures and arcs, we simply set their enhancement factors to 1, so that they are considered to have the average weatherability of continental crust (equivalent to the default SCION model).

## Paleogeography

Isolating the effects of paleogeography are more difficult and to properly test this would require alternative paleogeographic (and their dependent paleoclimate) models. As a simple measure to test the effect of paleogeography on hydrology and weathering, we alter the values of the slope and runoff to be equal to the mean value of the globe at each time step and CO₂ step. This acts to reduce any change in gross silicate weathering but instead spreads it across the entire land mask at any time step (rather than being concentrated in certain areas). The homogenization occurs before the calculation of any contribution of arcs, sutures, or relict arcs to the silicate weathering flux.

## Supplementary Materials

### The PDF file includes:

Figs. S1 to S15

Table S1

Legends for data S1 and S2

### Other Supplementary Material for this manuscript includes the following:

Data S1 and S2

## REFERENCES AND NOTES

1. R. A. Berner, *The Phanerozoic Carbon Cycle: CO₂ and O₂* (Oxford Univ. Press, 2004).

2. D. L. Royer, R. A. Berner, I. P. Montañez, N. J. Tabor, D. J. Beerling, CO₂ as a primary driver of Phanerozoic climate. *GSA Today* **14**, 4–10 (2004).

3. B. J. W. Mills, A. J. Krause, C. R. Scotese, D. J. Hill, G. A. Shields, T. M. Lenton, Modelling the long-term carbon cycle, atmospheric CO₂, and Earth surface temperature from late Neoproterozoic to present day. *Gondwana Res.* **67**, 172–186 (2019).

4. P. B. Kelemen, C. E. Manning, Reevaluating carbon fluxes in subduction zones, what goes down, mostly comes up. *Proc. Natl. Acad. Sci. U.S.A.* **112**, E3997–E4006 (2015).

5. N. R. McKenzie, B. K. Horton, S. E. Loomis, D. F. Stockli, N. J. Planavsky, C.-T. A. Lee, Continental arc volcanism as the principal driver of icehouse-greenhouse variability. *Science* **352**, 444–447 (2016).

6. R. D. Müller, B. Mather, A. Dutkiewicz, T. Keller, A. S. Merdith, C. M. Gonzalez, W. Gorczyk, S. Zahirovic, Evolution of Earth's tectonic carbon conveyor belt. *Nature* **605**, 629–639 (2002).

7. C.-T. A. Lee, S. Thurner, S. Paterson, W. Cao, The rise and fall of continental arcs: Interplays between magmatism, uplift, weathering, and climate. *Earth Planet. Sci. Lett.* **425**, 105–119 (2015).

8. T. M. Gernon, T. K. Hincks, A. S. Merdith, E. J. Rohling, M. R. Palmer, G. L. Foster, C. P. Bataille, R. D. Müller, Global chemical weathering dominated by continental arcs since the mid-Palaeozoic. *Nat. Geosci.* **14**, 690–696 (2021).

9. F. A. Macdonald, N. L. Swanson-Hysell, Y. Park, L. Lisiecki, O. Jagoutz, Arc-continent collisions in the tropics set Earth's climate state. *Science* **364**, 181–184 (2019).

10. J. C. G. Walker, P. B. Hays, J. F. Kasting, A negative feedback mechanism for the long-term stabilization of Earth's surface temperature. *J. Geophys. Res.* **86**, 9776–9782 (1981).

11. R. G. Hilton, A. J. West, Mountains, erosion and the carbon cycle. *Nat. Rev. Earth Environ.* **1**, 284–299 (2020).

12. D. J. Burdige, Burial of terrestrial organic matter in marine sediments: A re-assessment. *Global Biogeochem. Cycles* **19**, 1–7 (2005).

13. C. France-Lanord, L. A. Derry, Organic carbon burial forcing of the carbon cycle from Himalayan erosion. *Nature* **390**, 65–67 (1997).

14. R. A. Berner, The rise of plants and their effect on weathering and atmospheric CO₂. *Science* **276**, 544–546 (1997).

15. T. J. Algeo, R. A. Berner, J. B. Maynard, S. E. Scheckler, Others, late devonian oceanic anoxic events and biotic crises: “rooted” in the evolution of vascular land plants. *GSA Today* **5**, 45–66 (1995).

16. M. Paffre, Y. Godderis, A. Pohl, Y. Donnadieu, S. Carretier, G. Le Hir, The complex response of continental silicate rock weathering to the colonization of the continents by vascular plants in the Devonian. *Am. J. Sci.* **322**, 461–492 (2022).

17. C. M. Marcilly, T. H. Torsvik, M. Domeier, D. L. Royer, New paleogeographic and degassing parameters for long-term carbon cycle models. *Gondwana Res.* **97**, 176–203 (2021).

18. Y. Goddéris, Y. Donnadieu, G. Le Hir, V. Lefebvre, E. Nardin, The role of palaeogeography in the Phanerozoic history of atmospheric CO₂ and climate. *Earth-Sci. Rev.* **128**, 122–138 (2014).

19. Y. Park, P. Maffre, Y. Goddéris, F. A. Macdonald, E. S. C. Anttila, N. L. Swanson-Hysell, Emergence of the Southeast Asian islands as a driver for Neogene cooling. *Proc. Natl. Acad. Sci. U.S.A.* **117**, 25319–25326 (2020).

20. V. Lefebvre, Y. Donnadieu, Y. Goddéris, F. Fluteau, L. Hubert-Théou, Was the Antarctic glaciation delayed by a high degassing rate during the early Cenozoic? *Earth Planet. Sci. Lett.* **371–372**, 203–211 (2013).

21. N. M. Bergman, T. M. Lenton, A. J. Watson, COPSE: A new model of biogeochemical cycling over Phanerozoic time. *Am. J. Sci.* **304**, 397–437 (2004).

22. R. A. Berner, GEOCARBSULF: A combined model for Phanerozoic atmospheric O₂ and CO₂. *Geochem. Cosmochim. Acta* **70**, 5653–5664 (2006).

23. B. J. W. Mills, Y. Donnadieu, Y. Goddéris, Spatial continuous integration of Phanerozoic global biogeochemistry and climate. *Gondwana Res.* **100**, 73–86 (2021).

24. Y. Donnadieu, Y. Goddéris, R. Pierrehumbert, G. Dromart, F. Fluteau, R. Jacob, A GEOCLIM simulation of climatic and biogeochemical consequences of Pangea breakup. *Geochem. Geophys. Geosyst.* **7**, 10.1029/2006GC001278 (2006).

25. Y. Goddéris, Y. Donnadieu, S. Carretier, M. Aretz, G. Dera, M. Macouin, V. Regard, Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering. *Nat. Geosci.* **10**, 382–386 (2017).

26. M. E. Raymo, W. F. Ruddiman, Tectonic forcing of late Cenozoic climate. *Nature* **359**, 117–122 (1992).

27. T. M. Lenton, S. J. Daines, B. J. W. Mills, COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time. *Earth-Sci. Rev.* **178**, 1–28 (2018).

28. S. Finnegan, K. Bergmann, J. M. Eiler, D. S. Jones, D. A. Fike, I. Eisenman, N. C. Hughes, A. K. Tripati, W. W. Fischer, The magnitude and duration of late ordovician–early silurian glaciation. *Science* **331**, 903–906 (2011).

29. R. A. Berner, A model for atmospheric CO₂ sub 2 over phanerozoic time. *Am. J. Sci.* **291**, 339–376 (1991).

30. A. S. Merdith, S. E. Williams, A. S. Collins, M. G. Tetley, J. A. Mulder, M. L. Blades, A. Young, S. E. Armistead, J. Cannon, S. Zahirovic, R. D. Müller, Extending full-plate tectonic models into deep time: Linking the Neoproterozoic and the Phanerozoic. *Earth-Sci. Rev.* **214**, 103477 (2021).

31. M. T. Jones, D. A. Jerram, H. H. Svensen, C. Grove, The effects of large igneous provinces on the global carbon and sulphur cycles. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **441**, 4–21 (2016).

32. I. M. Fendley, J. Frieling, T. A. Mather, M. Ruhl, S. P. Hesselbo, H. C. Jenkyns, Early Jurassic large igneous province carbon emissions constrained by sedimentary mercury. *Nat. Geosci.* **17**, 241–248 (2024).

33. S. V. Sobolev, A. V. Sobolev, D. V. Kuzmin, N. A. Krivolutskaya, A. G. Petrunin, N. T. Arndt, V. A. Radko, Y. R. Vasiliev, Linking mantle plumes, large igneous provinces and environmental catastrophes. *Nature* **477**, 312–316 (2011).

34. P. B. Wignall, Large igneous provinces and mass extinctions. *Earth Sci. Rev.* **53**, 1–33 (2001).

35. B. J. W. Mills, S. J. Daines, T. M. Lenton, Changing tectonic controls on the long-term carbon cycle from Mesozoic to present. *Geochem. Geophys. Geosyst.* **15**, 4866–4884 (2014).

36. H. Svensen, S. Planke, A. G. Polozov, N. Schmidbauer, F. Corfu, Y. Y. Podladchikov, B. Jamtveit, Siberian gas venting and the end-Permian environmental crisis. *Earth Planet. Sci. Lett.* **277**, 490–500 (2009).

37. C. Ganino, N. T. Arndt, Climate changes caused by degassing of sediments during the emplacement of large igneous provinces. *Geology* **37**, 323–326 (2009).

38. J. Dal Corso, B. J. W. Mills, D. Chu, R. J. Newton, T. A. Mather, W. Shu, Y. Wu, J. Tong, P. B. Wignall, Permo–Triassic boundary carbon and mercury cycling linked to terrestrial ecosystem collapse. *Nat. Commun.* **11**, 2962 (2020).

39. D. L. Royer, Y. Donnadieu, J. Park, J. Kowalczyk, Y. Goddéris, Error analysis of CO₂ and O₂ estimates from the long-term geochemical model GEOCARBSULF. *Am. J. Sci.* **314**, 1259–1283 (2014).

40. P. Maffre, J.-B. Ladant, J.-S. Moquet, S. Carretier, D. Labat, Y. Goddéris, Mountain ranges, climate and weathering. Do orogens strengthen or weaken the silicate weathering carbon sink? *Earth Planet. Sci. Lett.* **493**, 174–185 (2018).

41. J. Gaillardet, B. Dupré, P. Louvat, C. J. Allègre, Global silicate weathering and CO₂ consumption rates deduced from the chemistry of large rivers. *Chem. Geol.* **159**, 3–30 (1999).

42. A. J. Boucot, C. Xu, C. R. Scotese, R. J. Morley, Phanerozoic paleoclimate: An atlas of lithologic indicators of climate (SEPM Society for Sedimentary Geology, 2013).

43. M. Rogov, V. Ershova, C. Gaina, O. Vereshchagin, K. Vasileva, K. Mikhailova, A. Krylov, Glendonites throughout the Phanerozoic. *Earth Sci. Rev.* **241**, 104430 (2023).

44. T. H. Torsvik, R. Van der Voo, U. Preeden, C. Mac Niocail, B. Steinberger, P. V. Doubrovine, D. J. J. van Hinsbergen, M. Domeier, C. Gaina, E. Tohver, J. G. Meert, P. J. A. McCausland, L. R. M. Cocks, Phanerozoic polar wander, palaeogeography and dynamics. *Earth Sci. Rev.* **114**, 325–368 (2012).
45. C. R. Scotese, H. Song, B. J. W. Mills, D. G. van der Meer, Phanerozoic paleotemperatures: The earth's changing climate during the last 540 million years. *Earth Sci. Rev.* **215**, 103503 (2021).

46. D. G. van der Meer, C. R. Scotese, B. J. W. Mills, A. Sluijs, A.-P. van den Berg van Saparoea, R. M. B. van de Weg, Long-term Phanerozoic global mean sea level: Insights from strontium isotope variations and estimates of continental glaciation. *Gondwana Res.* **111**, 103–121 (2022).

47. E. J. Judd, J. E. Tierney, D. J. Lunt, I. P. Montañez, B. T. Huber, S. L. Wing, P. J. Valdes, A 485-million-year history of Earth's surface temperature. *Science* **385**, eadk3705 (2024).

48. T.sson, S. Rauzi, Oxygen isotope ensemble reveals Earth's seawater, temperature, and carbon cycle history. *Science* **383**, 666–670 (2024).

49. J. Veizer, D. Ala, K. Azmy, P. Bruckschen, D. Buhl, F. Bruhn, G. A. F. Carden, A. Diener, S. Ebneth, Y. Godderis, T. Jasper, C. Korte, F. Pawellek, O. G. Podlaha, H. Strauss, $^{87}$Sr/$^{86}$Sr, δ¹³C and δ¹⁸O evolution of Phanerozoic seawater. *Chem. Geol.* **161**, 59–88 (1999).

50. J. Veizer, A. Prokoph, Temperatures and oxygen isotopic composition of Phanerozoic oceans. *Earth Sci. Rev.* **146**, 92–104 (2015).

51. E. L. Grossman, M. M. Joachimski, Ocean temperatures through the Phanerozoic reassessed. *Sci. Rep.* **12**, 8938 (2022).

52. B. B. Barney, E. L. Grossman, Supplemental Material: Reassessment of ocean paleotemperatures during the Late Ordovician (The Geological Society of America, 2022); https://doi.org/10.1130/GEOL.S.18822950.v1.

53. M. J. Cramwinckel, M. Huber, I. J. Kocken, C. Agnini, P. K. Bijl, S. M. Bohaty, J. Frielig, A. Goldner, F. J. Hilgen, E. L. Kip, F. Peterse, R. van der Ploeg, U. Röhl, S. Schouten, A. Sluijs, Synchronous tropical and polar temperature evolution in the Eocene. *Nature* **559**, 382–386 (2018).

54. G. A. Henkes, B. H. Passey, E. L. Grossman, B. J. Shenton, T. E. Yancey, A. Pérez-Huerta, Temperature evolution and the oxygen isotope composition of Phanerozoic oceans from carbonate clumped isotope thermometry. *Earth Planet. Sci. Lett.* **490**, 40–50 (2018).

55. D. Evans, N. Sagoo, W. Renema, L. J. Cotton, W. Müller, J. A. Todd, P. K. Saraswati, P. Stassen, M. Ziegler, P. N. Pearson, P. J. Valdes, H. P. Affek, Eocene greenhouse climate revealed by coupled clumped isotope-Mg/Ca thermometry. *Proc. Natl. Acad. Sci. U.S.A.* **115**, 1174–1179 (2018).

56. B. J. W. Mills, Hot and cold Earth through time. *Science* **385**, 1276–1278 (2024).

57. G. L. Foster, D. L. Royer, D. J. Lunt, Future climate forcing potentially without precedent in the last 420 million years. *Nat. Commun.* **8**, 14845 (2017).

58. C. R. Witkowski, J. W. H. Weijers, B. Blais, S. Schouten, J. S. Sinninghe Damsté, Molecular fossils from phytoplankton reveal secular Pco2 trend over the Phanerozoic. *Sci. Adv.* **4**, eaat4556 (2018).

59. Cenozoic CO2 Proxy Integration Project (CenCO2PIP) Consortium, B. Hönisch, D. L. Royer, D. O. Breecker, P. J. Polissar, G. J. Bowen, M. J. Henehan, Y. Cui, M. Steinhorsdottir, J. C. McElwain, M. J. Kohn, A. Pearson, S. R. Phelps, K. T. Uno, A. Ridgwell, E. Anagnostou, J. Austermann, M. P. S. Badger, R. S. Barclay, P. K. Bijl, T. B. Chalk, C. R. Scotese, E. de la Vega, R. M. DeConto, K. A. Dyez, V. Ferrini, P. J. Franks, C. F. Giulivi, M. Gutjahr, D. T. Harper, L. L. Haynes, M. Huber, K. E. Snell, B. A. Keisling, W. Konrad, T. K. Lowenstein, A. Malinverno, M. Guillermic, L. M. Mejia, J. N. Milligan, J. J. Morton, L. Nordt, R. Whiteford, A. Roth-Nebelsick, J. K. C. Rugenstein, M. F. Schaller, N. D. Sheldon, S. Sosdian, E. B. Wilkes, C. R. Witkowski, Y. G. Zhang, L. Anderson, D. J. Beerling, C. Bolton, T. E. Cerling, J. M. Cotton, J. Da, D. D. Ekart, G. L. Foster, D. R. Greenwood, E. G. Hyland, E. A. Jagniecki, J. P. Jasper, J. B. Kowalczyk, L. Kunzschner, C. E. Lawrence, C. H. Lear, M. A. Martínez-Botí, D. P. Maxbauer, P. Montagna, B. D. A. Naafs, J. W. B. Rae, M. Raitzsch, G. J. Retallack, S. J. Ring, O. Seki, J. Sepúlveda, A. Sinha, T. F. Tesfamichael, A. Tripati, J. van der Burgh, J. Yu, J. C. Zachos, L. Zhang, Toward a Cenozoic history of atmospheric CO₂. *Science* **382**, eadi5177 (2023).

60. D. O. Breecker, Z. D. Sharp, L. D. McFadden, Atmospheric CO₂ concentrations during ancient greenhouse climates were similar to those predicted for A.D. 2100. *Proc. Natl. Acad. Sci. U.S.A.* **107**, 576–580 (2010).

61. J. W. B. Rae, Y. G. Zhang, X. Liu, G. L. Foster, H. M. Stoll, R. D. M. Whiteford, Atmospheric CO₂ over the past 66 million years from marine archives. *Annu. Rev. Earth Planet. Sci.* **49**, 609–641 (2021).

62. K. Caldeira, J. F. Kasting, Susceptibility of the early Earth to irreversible glaciation caused by carbon dioxide clouds. *Nature* **359**, 226–228 (1992).

63. M. V. Caputo, R. O. B. dos Santos, Stratigraphy and ages of four Early Silurian through Late Devonian, Early and Middle Mississippian glaciation events in the Parnaíba Basin and adjacent areas, NE Brazil. *Earth-Sci. Rev.* **207**, 103002 (2020).

64. P. E. Isaacson, E. Díaz-Martínez, G. W. Grader, J. Kalvoda, O. Babek, F. X. Devuyst, Late Devonian–Earliest Mississippian glaciation in Gondwanal洲 and its biogeographic consequences. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **268**, 126–142 (2008).

65. M. Streel, M. V. Caputo, J. H. G. Melo, M. Perez-Leyton, What do latest Famennian and Mississippian miospheres from South American diamicites tell us? *Palaeobiodiversity and Palaeoenvironments* **93**, 299–316 (2013).

66. A. J. Boucot, X. Chen, C. R. Scotese, 'Chapter 8: Permian' in *Phanerozoic Paleoclimate: An Atlas of Lithologic Indicators of Climate* (SEPM Society for Sedimentary Geology, 2013).

67. Y. Goddéris, Y. Donnadieu, A sink-or a source-driven carbon cycle at the geological timescale? Relative importance of palaeogeography versus solid Earth degassing rate in the Phanerozoic climatic evolution. *Geol. Mag.* **156**, 355–365 (2019).

68. A. Young, N. Flament, S. E. Williams, A. S. Merdith, X. Cao, R. D. Müller, Long-term Phanerozoic sea level change from solid Earth processes. *Earth Planet. Sci. Lett.* **584**, 117451 (2022).

69. C. M. Marcilly, T. H. Torsvik, C. P. Conrad, Global Phanerozoic sea levels from paleogeographic flooding maps. *Gondwana Res.* **110**, 128–142 (2022).

70. C. R. Scotese, An atlas of phanerozoic paleogeographic maps: The sea comes in and the seas go out. *Annu. Rev. Earth Planet. Sci.* **49**, 679–728 (2021).

71. W. Cao, S. Zahirovic, N. Flament, S. Williams, J. Golonka, R. D. Müller, Improving global paleogeography since the late Paleozoic using paleobiology. *Biogeosciences* **14**, 5425–5439 (2017).

72. T. H. Torsvik, L. R. M. Cocks, *Earth History and Palaeogeography* (Cambridge Univ. Press, 2016); www.cambridge.org/core/books/earth-history-and-palaeogeography/076D4FCF CEE4CA6EF8EF3B449594DD50.

73. D. Zheng, A. Merdith, Y. Goddéris, Y. Donnadieu, K. Gurung, B. J. Mills, Using deep learning to integrate paleoclimate and global biogeochemistry over Phanerozoic time. *Geosci. Model Dev. Discuss.* **2024**, 1–20 (2024).

74. M. Domeier, E. Font, N. Youbi, J. Davies, S. Nemkin, R. Van der Voo, M. Perrot, M. Benabbou, M. A. Boumehdi, T. H. Torsvik, On the Early Permian shape of Pangea from paleomagnetism at its core. *Gondwana Res.* **90**, 171–198 (2021).

75. S. Wu, Z. Liu, J. Du, Y. Liu, Change of global ocean temperature and decadal variability under 1.5 °C warming in FOAM. *J. Mar. Sci. Eng.* **10**, 1231 (2022).

76. J. R. Zondervan, R. G. Hilton, M. Dellinger, F. J. Clubb, T. Roylands, M. Ogrić, Rock organic carbon oxidation CO₂ release offsets silicate weathering sink. *Nature* **623**, 329–333 (2023).

77. Z. Li, Y. G. Zhang, M. Torres, B. J. W. Mills, Neogene burial of organic carbon in the global ocean. *Nature* **613**, 90–95 (2023).

78. K. Gurung, K. J. Field, S. A. Batterman, S. W. Poulton, B. J. W. Mills, Geographic range of plants drives long-term climate change. *Nat. Commun.* **15**, 1805 (2024).

79. M. Domeier, A plate tectonic scenario for the lapetus and Rheic oceans. *Gondwana Res.* **36**, 275–295 (2016).

80. M. Domeier, V. Magni, M. W. Hounslow, T. H. Torsvik, Episodic zircon age spectra mimic fluctuations in subduction. *Sci. Rep.* **8**, 17471 (2018).

81. A. Young, N. Flament, K. Maloney, S. Williams, K. Matthews, S. Zahirovic, R. D. Müller, Global kinematics of tectonic plates and subduction zones since the late Paleozoic Era. *Geosci. Front.* **10**, 989–1013 (2019).

82. S. Zahirovic, A. Eleish, S. Doss, J. Pall, J. Cannon, M. Pistone, M. G. Tetley, A. Young, P. Fox, Subduction and carbonate platform interactions. *Geosci. Data J.* **9**, 371–383 (2022).

83. M. Domeier, Early Paleozoic tectonics of Asia: Towards a full-plate model. *Geosci. Front.* **9**, 789–862 (2018).

84. J. Hartmann, M. Moosdorf, The new global lithological map database GLIM: A representation of rock properties at the Earth surface. *Geochem. Geophys. Geosyst.* **13**, Q12004 (2012).

85. P. Mann, A. Taira, Global tectonic significance of the Solomon Islands and Ontong Java Plateau convergent zone. *Tectonophysics* **389**, 137–190 (2004).

86. K. Porkoláb, T. Duretz, P. Yamato, A. Auzemery, E. Willingshofer, Extrusion of subducted crust explains the emplacement of far-travelled ophiolites. *Nat. Commun.* **12**, 1499 (2021).

87. J. de Bremond d'Ars, C. Jaupart, R. S. J. Sparks, Distribution of volcanoes in active margins. *J. Geophys. Res.* **100**, 20421–20432 (1995).

88. J. A. West, Thickness of the chemical weathering zone and implications for erosional and climatic drivers of weathering and for carbon-cycle feedbacks. *Geology* **40**, 811–814 (2012).

89. R. A. Berner, GEOCARB II: A revised model of atmospheric CO₂ over Phanerozoic time. *Am. J. Sci.* **294**, 56–91 (1994).

90. S. Gaffin, Ridge volume dependence on seafloor generation rate and inversion using long term sealevel change. *Am. J. Sci.* **287**, 596–611 (1987).

91. K. J. Matthews, K. T. Maloney, S. Zahirovic, S. E. Williams, M. Seton, R. D. Müller, Global plate boundary evolution and kinematics since the late Paleozoic. *Glob. Planet. Change* **146**, 226–250 (2016).

92. M. W. Hounslow, M. Domeier, A. J. Biggin, Subduction flux modulates the geomagnetic polarity reversal rate. *Tectonophysics* **742–743**, 34–49 (2018).

93. H. Lee, J. D. Muirhead, T. P. Fischer, C. J. Ebinger, S. A. Kattenhorn, Z. D. Sharp, G. Kianji, Massive and prolonged deep carbon emissions associated with continental rifting. *Nat. Geosyst.* **9**, 145–149 (2016).

94. J. A. Hunt, A. Zafu, T. A. Mather, D. M. Pyle, P. H. Barry, Spatially variable CO₂ degassing in the main Ethiopian rift: Implications for magma storage, volatile transport, and rift-related emissions. *Geochem. Geophys. Geosyst.* **18**, 3714–3737 (2017).

95. C. Werner, T. P. Fischer, A. Aiuppa, M. Edmonds, C. Cardellini, S. Carn, G. Chiodini, E. Cottrell, M. Burton, H. Shinohara, P. Allard, "Carbon dioxide emissions from subaerial
volcanic regions" in *Deep Carbon: Past to Present* (Cambridge Univ. Press, 2019), pp. 188-236.

96. T. M. Seward, D. M. Kerrick, Hydrothermal CO₂ emission from the Taupo Volcanic Zone, New Zealand. *Earth Planet. Sci. Lett.* **139**, 105–113 (1996).

97. S. Brune, S. E. Williams, R. D. Müller, Potential links between continental rifting, CO₂ degassing and climate change through time. *Nat. Geosci.* **10**, 941–946 (2017).

98. A. M. C. Şengör, B. A. Natal'in, Rifts of the world, in *Mantle plumes: Their identification*, R. E. Ernst, K. L. Buchan, Eds. (The Geological Society of America, 2001).

99. A. S. Merdith, S. E. Williams, S. Brune, A. S. Collins, R. D. Müller, Rift and plate boundary evolution across two supercontinent cycles. *Glob. Planet. Change* **173**, 1–14 (2019).

100. R. Dasgupta, M. M. Hirschmann, The deep carbon cycle and melting in Earth's interior. *Earth Planet. Sci. Lett.* **298**, 1–13 (2010).

101. P. J. Gorman, D. M. Kerrick, J. A. D. Connolly, Modeling open system metamorphic decarbonation of subducting slabs. *Geochem. Geophys. Geosyst.* **7**, Q04007 (2006).

102. A. Aiuppa, T. P. Fischer, T. Plank, P. Bani, CO₂ flux emissions from the Earth's most actively degassing volcanoes, 2005–2015. *Sci. Rep.* **9**, 5442 (2019).

103. T. P. Fischer, Fluxes of volatiles (H₂O, CO₂, N₂, Cl, F) from arc volcanoes. *Geochem. J.* **42**, 21–38 (2008).

104. T. P. Fischer, S. Arellano, S. Carn, A. Aiuppa, B. Galle, P. Allard, T. Lopez, H. Shinohara, P. Kelly, C. Werner, C. Cardellini, G. Chiodini, The emissions of CO₂ and other volatiles from the world's subaerial volcanoes. *Sci. Rep.* **9**, 18716 (2019).

105. E. Mason, M. Edmonds, A. V. Turchyn, Remobilization of crustal carbon may dominate volcanic arc emissions. *Science* **357**, 290–294 (2017).

106. A. Aiuppa, T. P. Fischer, T. Plank, P. Robidoux, R. Di Napoli, Along-arc, inter-arc and arc-to-arc variations in volcanic gas CO₂/ST ratios reveal dual source of carbon in arc volcanism. *Earth-Sci. Rev.* **168**, 24–47 (2017).

107. T. Kagoshima, Y. Sano, N. Takahata, T. Maruoka, T. P. Fischer, K. Hattori, Sulphur geodynamic cycle. *Sci. Rep.* **5**, 8330 (2015).

108. T. Volk, Sensitivity of climate and atmospheric CO₂ to deep-ocean and shallow-ocean carbonate burial. *Nature* **337**, 637–640 (1989).

109. K. Caldeira, Enhanced Cenozoic chemical weathering and the subduction of pelagic carbonate. *Nature* **357**, 578–581 (1992).

110. D. V. Bekaert, S. J. Turner, M. W. Broadley, J. D. Barnes, S. A. Halldórsson, J. Labidi, J. Wade, K. J. Walowski, P. H. Barry, Subduction-driven volatile recycling: A global mass balance. *Annu. Rev. Earth Planet. Sci.* **49**, 37–70 (2021).

111. K. A. Campbell, J. D. Farmer, D. Des Marais, Ancient hydrocarbon seeps from the Mesozoic convergent margin of California: Carbonate geochemistry, fluids and palaeoenvironments. *Geofluids* **2**, 63–94 (2002).

112. B. Marty, I. N. Tolstikhin, CO₂ fluxes from mid-ocean ridges, arcs and plumes. *Chem. Geol.* **145**, 233–248 (1998).

113. M. M. Hirschmann, Comparative deep Earth volatile cycles: The case for C recycling from exosphere/mantle fractionation of major (H₂O, C, N) volatiles and from H₂O/Ce, CO₂/Ba, and CO₂/Nb exosphere ratios. *Earth Planet. Sci. Lett.* **502**, 262–273 (2018).

114. T. Keller, R. F. Katz, M. M. Hirschmann, Volatiles beneath mid-ocean ridges: Deep melting, channelised transport, focusing, and metasomatism. *Earth Planet. Sci. Lett.* **464**, 55–68 (2017).

115. J. M. Tucker, S. Mukhopadhyay, H. M. Gonnermann, Reconstructing mantle carbon and noble gas contents from degassed mid-ocean ridge basalts. *Earth Planet. Sci. Lett.* **496**, 108–119 (2018).

116. S. E. Bryan, A. Ewart, C. J. Stephens, J. Parianos, P. J. Downes, The Whitsunday Volcanic Province, Central Queensland, Australia: Lithological and stratigraphic investigations of a silicic-dominated large igneous province. *J. Volcanol. Geotherm. Res.* **99**, 55–78 (2000).

117. J. Longman, B. J. W. Mills, Y. Donnadieu, Y. Goddéris, Assessing volcanic controls on Miocene climate change. *Geophys. Res. Lett.* **49**, e2021GL096519 (2022).

118. F. Crameri, G. E. Shephard, P. J. Heron, The misuse of colour in science communication. *Nat. Commun.* **11**, 5444 (2020).

119. Y. Goddéris, Y. Donnadieu, V. Lefebvre, G. Le Hir, E. Nardin, Tectonic control of continental weathering, atmospheric CO₂, and climate over Phanerozoic times. *C. R. Geosci.* **344**, 652–662 (2012).

120. R. C. Blakey, "Gondwana paleogeography from assembly to breakup—A 500 m.y. odyssey" in *Special Paper 441: Resolving the Late Paleozoic Ice Age in Time and Space* (The Geological Society of America, 2008), vol. 441, pp. 1–28.

**Acknowledgments:** We thank A. Farnsworth, C. Goldblatt, and the reviewers for constructive criticism that improved many aspects of this study. **Funding:** This work was supported by the EU MSCA-IF 893615 (NEOEARTH) (to A.S.M. and B.J.W.M.), ARC DECRA DE230101642 (to A.S.M.), WoodNext Foundation (to T.G.), and NERC project NE/X011208/1 (RIFT-CC) (to A.S.M. and B.J.W.M.). **Author contributions:** Writing—original draft: A.S.M., B.J.W.M., Y.G, and P.M. Investigation: T.M.G. and A.S.M. Writing—review and editing: Y.D., R.D.M., T.M.G., A.S.M., J.L., B.J.W.M., and Y.G. Methodology: Y.D., T.M.G., A.S.M., J.L., B.J.W.M., and Y.G. Resources: Y.D., A.S.M., and B.J.W.M. Funding acquisition: T.M.G., A.S.M., and B.J.W.M. Data curation: A.S.M. Validation: R.D.M., T.M.G., A.S.M., B.J.W.M., Y.G., and P.M. Supervision: B.J.W.M. Formal analysis: A.S.M. and J.L. Software: A.S.M., J.L., and B.J.W.M. Project administration: A.S.M. and B.J.W.M. Visualization: A.S.M. and B.J.W.M. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** The data used and generated in our study are available at https://zenodo.org/record/7940113 or from the original referenced publications discussed in text. Figures 1 to 4 were produced using perceptually uniform color maps (118). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

Submitted 15 November 2023

Resubmitted 9 July 2024

Accepted 14 January 2025

Published 14 February 2025

10.1126/sciadv.adm9798