Understanding the early Paleozoic carbon cycle balance and climate change from modelling

The Ordovician global cooling trend observed by several temperature proxies, which coincides with one of the most signiﬁcant evolutionary diversiﬁcations on Earth, is yet to be fully understood. This study presents new simulations of p CO 2 and surface temperatures using a spatially resolved climate-carbon cycle Earth system model fed with reﬁned continental reconstructions and new estimates of solid Earth degassing. First, in order to quantify the respective roles of paleogeography and degassing, the impact of continental drift alone is investigated. This is done by calculating its imprint on continental weathering rates throughout the Ordovician under a constant p CO 2 , and by testing different topography inputs. Secondly, the sensitivity of the Earth’s climate to paleogeography and carbon degassing changes is investigated by coupling the climate and carbon cycles in the GEOCLIM model. Based on our experiments we show that, although Early Ordovician high temperatures can be replicated within error margins, our new constraints cannot explain the intense cooling over the Mid to Late Ordovician, even if a progressive enhancement in Earth surface weatherability is taken into account. By using GEOCLIM in an inverse modelling approach, we calculate that the theoretical degassing necessary to reach proxy-derived temperatures for the Early Ordovician is three-to-ﬁve times higher compared to modern values. Further, in order to simulate the following Ordovician cooling trend, the solid Earth degassing must be reduced to modern-day values in only 30 Myrs. We conclude that, if accepting the veracity of the high Early Ordovician temperatures, alternative sources but also sinks of carbon must be considered to explain the climatic shift over the period.


Introduction
The Ordovician long-term global cooling trend is a key feature of the early Paleozoic climate (Trotter et al., 2008), and is also observed in recent studies (e.g.Edwards et al., 2021;Goldberg et al., 2021;Song et al., 2019).Global cooling initiated in the late Cambrian, and proceeded rapidly during the Early Ordovician, before tropical sea-surface temperatures (SSTs) reached * Corresponding author.
Recent studies suggest an intricate coupling between climate and marine biological evolution throughout the Ordovician.The Great Ordovician Biodiversification Event (GOBE) marked the radiation of the marine animal phyla that had emerged during the Cambrian Explosion (Cocks and Torsvik, 2020).Biodiversification may have been driven by cooling to more habitable temperatures alone (Trotter et al., 2008;Rasmussen et al., 2016), or in conjunction with the ensuing ocean oxygenation (Edwards et al., 2017), both being known to act synergistically on marine ectotherm habitability (Deutsch et al., 2015;Stockey et al., 2021).Ordovician cooling was not only associated with a significant increase in the number of different complex animal forms during the Mid Ordovician (∼470-455 Ma; GOBE), but the later Hirnantian cooling (with SSTs averaging 22 ± 6 • C) coincided with the only known mass extinction to have occurred during icehouse conditions (Torsvik et al., 2021).
The cause(s) of the long-term Ordovician cooling trend remain(s) poorly understood due to the sparsity of data and the large spread of temperature estimates (Fig. 1a).First, temporal trends from oxygen isotope data (δ 18 O) may be distorted by diagenesis, although meticulous screening for pristine samples and trace element analyses and the use of material supposedly resistant to diagenesis (such as biogenic apatite) permit gaining confidence in reconstructed temperatures (Hearing et al., 2018).Second, the heterogeneous spatial sampling of proxy data can create an underestimation of climatic perturbations, which has been shown for the Ordovician (Jones and Eichenseer, 2021).Absolute temperature values are even more challenging to constrain compared to temporal trends, due to large error margins associated with functions used to convert isotopic values to temperature (Lécuyer et al., 2013), as well as uncertain constraints on early Paleozoic seawater isotopic compositions.Some studies suggest that the seawater isotopic composition has remained unchanged (e.g.Henkes et al., 2018) while others follow the secular trend idea of the Veizer hypothesis (Veizer and Prokoph, 2015) of a gradual enrichment of the 18 O isotope (e.g., Galili et al., 2019;Hodel et al., 2018;Defliese, 2020).
Several hypotheses have been proposed to explain the Ordovician long-term cooling, either by increasing the weatherability of the continents or by decreasing the amplitude of the solid Earth degassing.On the one hand, Earth's surface weatherability increases when continental drift brings continents into the warm and humid equatorial area (Nardin et al., 2011, Fig. 1b) in response to early land-plant colonization (Lenton et al., 2012;Porada et al., 2016), or, for example, through increased weathering of mafic silicate rock linked to physical erosion through changes in topography (Kump et al., 1999;Young et al., 2009;Swanson-Hysell and Macdonald, 2017;Macdonald et al., 2019).On the other hand, volcanic outgassing rates may have decreased through the early Paleozoic (McKenzie et al., 2016).The search for a unique causative mechanism for the long-term Ordovician cooling trend is elusive and these various factors may all have played a role.Unfortunately, these different hypotheses have been proposed in isolation using different methods and tools, and the importance and contribution of each of these factors remain to be quantified in an integrated numerical framework.
Here, we use the climate-carbon model GEOCLIM (Donnadieu et al., 2006-updated) to simulate the long-term Ordovician climatic trend in response to changes in paleogeography and associated silicate weathering activity as well as changes in volcanic degassing rates.
The model is fed by climatic fields that have been simulated using the atmospheric general circulation model LMDz6 (see Supple-  2021) and represents the average of the least diagenetically altered samples.The purple shaded area represents the spread of data of the mid-to-least diagenetically altered data presented in Fig. 4 of Goldberg et al. (2021).Finally, the green dots are the temperatures reconstructed by Trotter et al. (2008) with their error band in shaded green.(b) Model simulations of atmospheric CO 2 levels.The green line is simulated using the forward model COPSE (Mills et al., 2019) whereas the stippled blue curve is obtained using GEOCARBSULF with updated paleogeography and degassing parameters (Marcilly et al., 2021).The dotted black curve is from Nardin et al. (2011) using GEOCLIM with the GCM FOAM.The numbers refer to the theoretical glacial inception levels with 3 and 6 • C climate sensitivity assuming 500 ppm inception level for modern times and a fainter sun in deep time (Marcilly et al., 2021).(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)mentary material) based on improved continental reconstructions (Marcilly et al., 2021).The GEOCLIM model accounts for the competing impacts of continental weathering and volcanic outgassing on the atmospheric partial pressure of CO 2 (pCO 2 ) and climate.This allows us to quantify the magnitude and temporal trend in global cooling, resulting from these competing mechanisms, and evaluate their relative contributions.In order to quantify the respective roles of paleogeography and degassing, we first explore the impact of continental drift alone, by calculating its imprint on continental weathering rates throughout the Ordovician under a constant pCO 2 .Then, the sensitivity of the Earth's climate to paleogeography and carbon degassing variations is investigated by coupling the climate and carbon cycles in the GEOCLIM model.This study highlights the current limits of our understanding of factors driving the early Paleozoic climate and presents conceivable directions for future studies.

Response of weathering to the continental configuration
To explore the interplays between Ordovician paleogeography and chemical weathering, we ran the LMDz6 climate model (see Supplementary Information) with late Cambrian to early Silurian paleogeographic maps ranging from 490 Ma to 430 Ma at 10 Myr intervals (7 time slices from Marcilly et al., 2021, see Supplementary Fig. S1).Several models for the tectonic and paleogeographic evolution exist for the Ordovician (e.g.Scotese, 2016;Blakey, 2009;Swanson-Hysell and Macdonald, 2017), but we here chose to work with the global model developed by Torsvik andcollaborators (2012, 2016), which forms the basis for the exposed land maps produced in Marcilly et al. (2021; see details in supplementary section G).
In order to isolate the paleogeographical imprint, the pCO 2 in our simulations is fixed at 12 PAL (3360 ppmv) (i.e., Pre-industrial Atmospheric Level; 1 PAL = 280 ppmv) to avoid strongly nonlinear climate responses associated with sea-ice feedbacks encountered at lower pCO 2 levels (Pohl et al., 2014(Pohl et al., , 2016)).We vary the model solar luminosity with geological age (from ∼1313 W m −2 at 490 Ma to ∼1320 W m −2 at 430 Ma) according to the model of stellar evolution (Gough, 1981).
At the global scale, the global mean surface temperature (GMST) remains relatively constant at ∼22.5 • C (Fig. 2a) in our simulations.However, the continental climate (i.e., temperature, runoff, and discharge) is more sensitive to paleogeographical changes (Fig. 2a, b; Fig. 3).The most important observations include significant continental warming between 490 Ma and 480 Ma (Fig. 2a), as well as local maxima of global dryness at 480 Ma and 450 Ma (Fig. 2b, Fig. S2 in SI).From 490 to 480 Ma, the global mean temperature over the continents increases by 4 • C (from 14.5 • C to 18.6 • C), followed by fluctuations between 17 • C and 19.3 • C during the rest of the investigated time period (Fig. 2a).
The warmest episodes occur at 480 Ma, 460 Ma, and 430 Ma when the Gondwana continent extends less over the South Pole region.
In the tropical band (i.e., from 21 • S to 21 • N), the annual mean temperature increases by 2 • C (from 30.1 • C to 32.1 • C) from 480 to 430 Ma (Fig. 2a), mainly due to the solar constant increase.Despite the large fluctuations in emerged land surface, reaching a peak of ∼36 × 10 6 km 2 in the tropics by 450 Ma (40% of the total land exposed, Table S1 in SI), the paleogeography has apparently little influence on the modelled terrestrial temperature in the tropics (Fig. 2a).
In the subtropics, arid zones are larger at 480 Ma and 450 Ma compared to other periods (Fig. 3, Supplementary Fig. S2), which is when the surface of the Gondwana continent located below the descending branch of the Hadley circulation is at its maximum.The increase in aridity (i.e., areas where mean precipitation minus evaporation is inferior to 0) is also favoured by the drift of smaller blocks (Southern part of South-China and Baltica) entering into the subtropics between 460 and 440 Ma (SI Figs.S1,  S2).
Precipitation and runoff reach maximum values over scattered islands in the equatorial belt and equatorial Gondwana (Fig. 3).A second belt with significant runoff is shown by all model runs at mid to high latitudes (i.e., between 40 • S and 80 • S), reaching maximum values over the western coast of Gondwana and the western edge of mountain ranges (Fig. 3; see also Pohl et al., 2017).The slight decoupling between the mean continental runoff and globally cumulated continental discharge at 450-440 Ma is due to an increase in the total area of emerged (exposed) lands between 490  Overall, these results demonstrate a long-term increase in continental weathering in response to continental rearrangement during the Ordovician rather than a sudden rise at the end of the Ordovician (Fig. 2c & 4).This progressive enhancement of Earth's surface weatherability (see also Marcilly et al., 2021) should lead to a pCO 2 reduction, a point explored in the following section.(Herrmann et al., 2004;Lowry et al., 2014;Pohl et al., 2016).The numbers refer to the theoretical glacial inception levels with 2, 3 and 6 • C climate sensitivity assuming 500 ppm inception level from modern times and fainter sun activity in deeper time (Marcilly et al., 2021).b) Global mean surface temperature (GMST) of GEOCLIM, run with constant degassing (black) or variable degassing from SBF and Z (pink and turquoise).The temperature evolution at fixed CO 2 (12 PAL or 3360 ppmv) is added (orange) for comparison.Climatic trends shown with background colours in panel b are extracted from proxy compilations of different studies (Trotter et al., 2008;Song et al., 2019;Scotese et al., 2021;Goldberg et al., 2021) and simplified as follows: cooling until 480 Ma followed by rapid cooling until 450-440 Ma followed by a plateau.

Climate-carbon response to continental drift
We explore the climatic evolution over the entire Ordovician assuming a constant degassing rate of the Earth.We couple the LMDz6 climate model with a long-term carbon cycle model to compute the steady-state climate, based on the climate-carbon negative feedback (Walker et al., 1981).The resulting climatecarbon model, GEOCLIM (Donnadieu et al., 2006-updated), is described in the Supplementary Information and available online.The calculated CO 2 consumption by silicate rock weathering is calibrated on the present-day state: this means that under present-day conditions, the model predicts the silicate weathering flux of the main world rivers (Park et al., 2020).The corresponding global silicate weathering (3.3 Tmol of CO 2 consumed each year) is close to the estimate by Gaillardet et al. (1999) (2.5 Tmol/yr).This implies that the solid Earth degassing must be close to the CO 2 consumption by silicate weathering.To reach a steady state at present day, we assume that the magmatic degassing of the present day Earth is equal to 3.3 Tmol/yr of CO 2 .This value is close to the CO 2 degassing estimates by Gerlach (2011) (3.4 to 5.9 Tmol CO 2 /yr).Therefore, to highlight the effect of paleogeography alone, the solid Earth degassing rate for all experiments was set to 7.1 Tmol C/yr and kept constant.This value is equivalent to about 2 times that of the present-day, which corresponds to the average degassing level over the studied period from proxies for seafloor spreading (the proxies and corresponding values are presented in section 3.1).As shown in Fig. 5a (black line), the pCO 2 drops from 30 to 18 PAL (i.e., 8400-5040 ppm) between 490 and 430 Ma (Table S3 in SI), in response to the increasing surface weatherability (Fig. 2c).
The substantial drop in pCO 2 observed between 480 Ma and 470 Ma induces an important global cooling event, with GMST decreasing by 2.5 • C (Fig. 5b; black line).Between 470 Ma and 430 Ma, although pCO 2 slightly decreases (Fig. 5a), pCO 2 -driven cooling is counteracted by the paleogeographical evolution accompanied by an increase in solar luminosity which, everything else unchanged, impose a warming from 470 Ma to 430 Ma (Fig. 5b, orange curve).As a result, the modelled global temperature slightly increases from 470 Ma to 430 Ma (Fig. 5b; black curve).Therefore, in our model, the paleogeographical forcing alone cannot explain the long-term Ordovician cooling trend reported based on geological data (Goldberg et al., 2021;Song et al., 2019;Trotter et al., 2008).

Climate-carbon response to enhanced topography
In order to test the effect(s) of changing topography, we ran three sets of simulations with varying topographic inputs.The first simulation (i.e., run #1, "revised topo" in Fig. 6) aimed to test the effect of higher altitudes, and so the altitudes were doubled compared to our previous simulations.Secondly, the influence of more substantial island arcs and obducted arcs were tested (Macdonald et al., 2019).To do so, new maps with increased relief areas and same doubled altitudes as in run #1 were created (i.e., run #2, "increased topo" in Fig. 6; see SI Fig. S7 for maps).The third set of simulations (i.e., run #3, "increased topo + slope" in Fig. 6) are similar to the second, but with increased slopes to accentuate the effect of increased topography on physical erosion, and hence on chemical weathering (see SI Fig. S8).All tests were first performed with constant degassing in order to isolate the topography effect (Fig. 6a).
The higher altitude simulations (run #1) led to similar temperatures compared to the previous simulations with low topography, even slightly higher between 460-440 Ma, and no intensification of the cooling can be observed (Fig. 6a).As the influence of topography is increased by taking into account larger accretion zones and island arcs (run #2), the temperature at 440 Ma is slightly lowered (Fig. 6a), but almost negligible as showing only a −0.18 • C difference.Even with steepened slopes (run #3), the difference only reached −0.51 • C (Fig. 6a).Contrarily to our first simulations with constant degassing and lower topography (Section 2.2), which suggest a warming from 450 Ma, all the simulations with enhanced topography suggest cooling from 460 Ma and lasting until 440 Ma.
These results show that changing topography does not have a large influence on the simulated Ordovician climatic trend.Taking into account larger topography and steeper slope does make the Earth slightly cooler (by ∼1 • C compared to the revised but not increased topography), but quasi uniformly on the whole period, since there is no significant temporal variation of global topography (growth nor decay) over the studied period.Therefore, the climatic trend is virtually unchanged, except for a slight prolonged cooling at 450-430 Ma in the increased topography and slope simulation.However, it is important to note that only modest changes in topography have been tested, as this parameter is poorly constrained through time.

Testing a variable degassing throughout the Ordovician
The Ordovician long-term cooling trend could potentially be attributed to fluctuations in the major long-term carbon sources.However, as direct measurements for solid Earth degassing in the Ordovician do not exist, we must rely on proxies.The long-term solid Earth degassing mainly occurs at spreading ridges, volcanic arcs and continental rifts (Kelemen and Manning, 2015).Reconstructions of degassing fluctuations through time have traditionally been based on the past activity of seafloor spreading, and such estimates are commonly used in long-term carbon cycle models (Berner, 2006;Goddéris and Donnadieu, 2019;Mills et al., 2019;Pohl et al., 2020).
Here we use two proxies for seafloor production rate, which in turn act as a proxy for solid Earth degassing (Fig. 7).The first is based on estimates of subduction fluxes from the full-plate model of Domeier and Torsvik (2014) extended back to the early Paleozoic (Marcilly et al., 2022) (see SI for more information).Note that the early Paleozoic estimates have large uncertainties since all the oceanic lithosphere prior to the Jurassic has been subducted and absolute longitudes in the Paleozoic are strongly debated.We here estimate uncertainties associated with the subduction flux based on previous versions of our subduction flux estimates (Fig. 7).The second proxy is the age distribution of zircon presented in Marcilly et al. (2021).Considering the increasing uncertainty of subduction flux estimates with deeper time, the age-frequency relationship of detrital zircons, mimicking variations in the subduction flux, both regionally and globally has been suggested as a more objective proxy for seafloor production (Domeier et al., 2018;McKenzie et al., 2016;Marcilly et al., 2021).Uncertainties for this proxy are based on its difference with subduction flux estimates in more recent time (Fig. 7).This zircon-based proxy suggests solid Earth degassing values in the same order of magnitude compared to other proxies commonly used in long-term carbon cycle modelling (e.g., Berner, 2006;Mills et al., 2019;Goddéris and Donnadieu, 2019;Merdith et al., 2021; Fig. 7), although it suggests lower estimates towards the end-Ordovician and a slight decrease between 480 Ma and 430 Ma.In contrast, the subduction-flux-based proxy we use here suggests higher values, especially between 490 and 470 Ma.
By introducing an up-to-date variable degassing reconstruction from subduction fluxes in the GEOCLIM simulations, the GMST temporal trend is amplified (Fig. 5b; purple curve) compared to simulations performed with a solid degassing held constant (Fig. 5b; black curve).While the Early Ordovician becomes warmer compared to the constant degassing scenario (T max of 28.5 • C vs. 26 • C at 480 Ma; Fig. 5b), the model now predicts a cooler Late Ordovician climate (22 • C vs. 24.5 • C at 440 Ma; Fig. 5b).However, even by combining paleogeographical evolution and solid Earth degassing variations, our model fails to reproduce the full extent of the observed tropical SST drop derived from proxy data (Song et al., 2019;Goldberg et al., 2021;Fig. 8).The lower degassing scenario derived from the zircon-based proxy produces simulations generally more in line with proxy temperatures (Fig. 8c).However, the high temperatures of the Early Ordovician are well out of reach, with a maximum SST modelled at 32.5 • C at 480 Ma against ∼40 • C from proxies (Fig. 8c).In summary, the simulation based on the subduction flux degassing matches more closely the proxy data between 490-480 Ma while the zircon proxy degassing curve is a better match from 470 Ma and onwards (Fig. 8b, & 8c), even if the total cooling intensity is only ∼2.5 • C (Fig. 8c).
Increasing the topography with varying degassing does not lead to a more pronounced global cooling (Fig. 6b,c), and contrarily to simulations with constant degassing, this trend is not changed nor intensified between 460 and 440 Ma.Here as well, the change in values at 440 Ma is negligible with a maximum change of only −0.42 • C with subduction flux degassing (Fig. 6b) and −0.53 • C with zircon degassing (Fig. 6c).
Overall, our simulations match the Early Ordovician cooling trend, but cannot explain why the cooling continues after the Middle Ordovician (∼460 Ma) and lasted for more than 40 Myrs (Fig. 8).It should be noted that, although various models and temperature reconstructions point towards a Mid-Ordovician cooling (Trotter et al., 2008;Rasmussen et al., 2016) there are no known Mid-Ordovician glacial deposits that can confirm this (Cocks and Torsvik, 2020).Furthermore, while the proxies suggest a total cooling of more than 10 • C (Fig. 8), our results show that the tectonic settings (from subduction fluxes) can drive a cooling close to 4 • C.This is in agreement with regional δ 18 O data from Baltica, suggesting a drop of 4-5 • C during the Middle Ordovician (Rasmussen et al., 2016).

Linking the δ 18 O record, paleogeography, and degassing
Regarding the conundrum of the Mid-Late Ordovician cooling, we complete our study by permitting the model parameters to reconcile modelled and proxy-derived SSTs.For that purpose, we adopt an inverse modelling approach by prescribing a tropical SST from the proxy reconstructions as a target in our model and calculating the degassing rate required to reproduce this Ordovician SST trend.We consider two alternative SST targets: after Goldberg et al. (2021) or after Song et al. (2019), both here averaged over 10 Myrs to fit our reconstruction intervals (see Fig. 1a & 8).The latter estimates have been checked to reflect only tropical temperatures (Marcilly et al., 2021).Further, Late Cambrian to early Silurian SSTs are almost exclusively derived from (apatite) conodonts that once lived from 12 • N to 23 • S (mean paleolatitude = 13± 5 • N/S).These SSTs should therefore reflect a reasonably well constrained cooling trend over the Ordovician with respect to the underestimation from biased geographical sampling described by Jones and Eichenseer (2021).
The two sets of temperature estimates differ mostly during the Late Cambrian-Early Ordovician, but following this time period, their trend and absolute values are similar and within error margins of each other (Figs.1a, 8).While the estimates by Song et al. (2019) are based on phosphatic fossils (conodonts), the ones from Goldberg et al. (2021) are from bulk rock isotopes.The scarcity of estimates in both studies may represent a limitation mostly for 490 Ma where only a few points are analyzed.In particular, the spread of the estimates is the largest for the Ordovician between 490-480 Ma in Goldberg et al. (2021), and therefore might represent the least constrained interval of their study as all data for this age is extracted from the same location (i.e., Newfoundland at the equator at the time).Similarly, the data reported in Song et al. (2019) become depleted of data coverage for 490 Ma after filtering and therefore reflect only an extrapolation of the continuity between the late Cambrian and 480 Ma.Therefore, we cannot rule out a local anomaly such as a temperature low as reconstructed in Goldberg et al. (2021).
The degassing rate scenarios obtained using the inverse approach differ significantly compared to those previously used in this study (Fig. 9a).For example, the inverse approach suggests that 10.7 to 11.8 Tmol C/yr would be required in order to match the observed mean tropical SST at 480 Ma (i.e., 38-41 • C; Song et al., 2019;Goldberg et al., 2021), which is significantly higher compared to the solid Earth degassing rate set in our previous experiments (averaging to 7.1 Tmol C/yr; see section 2.2).On the other hand, to simulate temperatures cooler than 30 • C at 450 Ma, the solid Earth degassing must be reduced to modern-day values (Fig. 9a).

Discussion
The existence of extremely high temperatures in the Late Cambrian-Early Ordovician remains uncertain as interpretations of observed δ 18 O data are debated.Two models for the early Paleozoic temperature exist; the "hot" and "cool" models (Scotese et al., 2021).Studies supporting the "cool" model suggest that the longterm trend in seawater δ 18 O is due to an enrichment of the heavier oxygen isotope ( 18 O) through geological time and does therefore not reflect change of the water temperature.This model, often referred as the "Veizer hypothesis" (Veizer and Prokoph, 2015), argues that early Paleozoic seawater was enriched in the light oxygen isotope ( 16 O) because of this evolution and that the high temperatures of the early Paleozoic are not accurate (Fig. 10).The main obstacle to this hypothesis is the lack of a primary mechanism that could trigger such a change, but also that other records that do not show the same feature.Those records support the "hot" model, which builds on the fact that the isotopic composition of seawater has remained constant so that the lighter δ 18 O values accurately reflect warmer paleo-temperatures during the early Paleozoic (e.g.Henkes et al., 2018;Bergmann et al., 2018) (Fig. 10a,b).Moreover, the peak in SST during this interval followed by a sharp cooling is well in agreement with the Cambrian radiation and mostly the GOBE demonstrating a strong (negative) correlation between biodiversity expansion and temperature (Cocks and Torsvik, 2020).If the "hot" model is correct, the state of the carbon mass balance during the Early Ordovician must have been largely different from the modern one.Indeed, the SST reconstructed from sedimentary δ 18 O records considering the "hot" hypothesis are critically high for the Early Ordovician and difficult to reproduce with a 3D-climate model.Tropical temperatures above 40 • C require very high atmospheric CO 2 .Those levels will necessarily depend on the sensitivity of the climate model to atmospheric CO 2 (SI Fig. S5).But regardless of what this sensitivity was in the deep past, which has been estimated to be ∼3 • C per CO 2 doubling during ice-free periods (greenhouse climates) and ∼6 • C per CO 2 doubling in the presence of land ice (Royer, 2016), or fluctuating from 2 to 5 • C (IPCC working group 1: Solomon et al., 2007), simulating temperatures above 40 • C in the tropical ocean will necessarily require very high levels of CO 2 .To generate a very warm tropical ocean, CO 2 levels as high as 300 times the present-day partial pressure are required in our models.By comparison, such levels are close to the post-snowball glaciation levels (Le Hir et al., 2008).We consider such CO 2 levels to be unrealistically high for the early Ordovician, however, they cannot be ruled out, because no reliable CO 2 proxies are extended further in the past than the Late Ordovician.
By targeting the proxy-based temperatures in our model, we have shown that the high CO 2 values prevailing in the early Ordovician would require a release flux of carbon into the oceanatmosphere system of about 3 to 5 times that of the present-day volcanic outgassing (Fig. 9a).A factor of 3 to 5 decreasing to only 1 over 30 Myrs (Fig. 9a) would require drastic, and perhaps unrealistic, changes in the geodynamic settings since the largest reduction in degassing reconstructed from proxies is by a factor of 2.2 (observed over the past 120 Myrs; Marcilly et al., 2021).
If we consider the proxy record as being biased (such as in the Veizer hypothesis; Fig. 10), then the calculated temperatures become significantly lower for the Ordovician.If our models with both degassing scenarios fit this "cool" model until 470 Ma, the simulated temperatures plot well above the data of Veizer following this time period (Fig. 10b).In order to simulate such temperatures, an even sharper decrease in degassing needs to be invoked, which seems unrealistic considering the tectonic activity of the time period and that the zircon degassing proxy already suggests the lowest estimates among different reconstruction methods.
Reconstructing actual CO 2 outgassing levels of the early Paleozoic is challenging.Considering the tectonic settings of the early Ordovician (490-480 Ma), with the initial rifting of the Rheic Ocean between Gondwana and Avalonia, a high rate of seafloor spreading is plausible.A change in volcanic outgassing has previously been argued to represent the principal driver for Ordovician climate change (McKenzie et al., 2016).However, no full-plate models, nor other methods to assess seafloor spreading (e.g., age-zircon distribution, sea level inversion, subduction zone lengths, continental rift lengths, or combinations of these) have ever suggested such high rates during the Phanerozoic or even the Neoproterozoic (Goddéris and Donnadieu, 2019;Marcilly et al., 2021;Mills et al., 2019).
Adding fluxes from continental rifts (Brune et al., 2017) and explicitly considering the intersection of continental magmatic arcs with carbonate-filled basins (e.g., Lee and Lackey, 2015) could certainly increase solid Earth degassing during the early Ordovician, but an additional release of carbon required to sustain high temperatures may come from other sources than solid Earth degassing.Imbalances in the organic carbon cycle may theoretically release CO 2 into the ocean and atmosphere if the oxidation of reduced sediments exposed on the continents exceeds the burial of organic carbon.The overall increase in the carbonate δ 13 C from the base of the Ordovician to the Upper Ordovician sections suggests a progressive predominance of burial versus oxidation (Veizer et al., 1999;Edwards and Saltzman, 2016;Hu et al., 2021) which is compatible with a decreasing carbon release from the Early to the Late Ordovician.However, the amplitude of the decrease needed to explain the Ordovician cooling is extremely large.If we assume that the "missing" Early Ordovician carbon source originates from the organic carbon sub-cycle instead of the outgassing flux, then the early Paleozoic must have been characterized by a very low organic carbon burial rate, or by significant oxidation of reduced sediments on the continents.Simple mass balance calculations show that if this were the case, the δ 13 C of seawater would be as low as −20 PDB in the Early Ordovician, far below the observed value of about −2 PDB.Finally, another process able to release carbon into the ocean-atmosphere system is the oxidative weathering of sedimentary pyrite exposed in the uplifted regions.This process generates sulfuric acid, which dissolves carbonate rocks, releasing carbon in the rivers and ultimately in the ocean (Torres et al., 2014).However, it has been shown that this carbon source is only transitory and cannot drive the long-term evolution of the carbon cycle (Maffre et al., 2020).Similar to the high temperatures of the Early Ordovician, our simulations do not replicate the cooling trend reconstructed based on proxy data (Figs. 1 & 8b,c).However, it should be noted that only the mean tropical SST has been considered, while proxy data have been collected at various geographical locations.Accounting for the spread of simulated SSTs at various longitudes within the tropical band reduces the model-data discrepancy (Fig. 8b,c).
However, the corresponding ±1 σ SST envelope, more representative of common SST conditions found in the simulated climatic states, remains far from capturing the proxy-derived temperatures (Fig. 8b,c).Moreover, when invoking changing degassing, the correlation between our models and the latest temperature reconstructions (i.e., Song et al., 2019;Scotese et al., 2021& Goldberg et al., 2021) is reasonable with an average Pearson correlation coefficient of 0.85 for both scenarios (Table 1).
Our simulations reveal that the paleogeographical evolution on its own cannot explain observed Ordovician cooling, and if no change in solid Earth degassing is accounted for, the continen-tal drift and masses distribution imposes a warming from 470 Ma onward (Fig. 5b).This is highlighted by the very low Pearson coefficient correlation between the latest temperature reconstructions and our model (0.2; Table 1).This observation differs slightly from the findings of Nardin et al. (2011), even if, globally, both studies show a similar flat trend in temperature (i.e., maximum variation of ∼3 • C).However, these authors used the GEOCLIM model with a fixed degassing, as well as a simpler representation of continental weathering.Their climatic fields were derived from a previous generation of the model, based on a lower spatial resolution (FOAM) of different continental reconstructions (Nardin et al., 2011).Further, the temporal resolution of this study was also lower compared to our study, with one GEOCLIM simulation every 20 Myrs (vs. 10 Myrs).Discrepancies between Nardin et al. (2011) and the results presented here reflect therefore differences in the continental reconstructions, the climate model, and the representation of continental weathering (Park et al., 2020).The major difference between these two GEOCLIM studies is that the model Constraining pCO 2 thresholds for the glacial onset during the Ordovician has been the topic of a range of studies, going back to Crowley and Baum (1991).Some of these studies coupled climate models to land-ice models, which effectively simulate ice-sheet nucleation (Herrmann et al., 2004;Lowry et al., 2014;Pohl et al., 2016).Due to differences in the models and boundary conditions, values reported by these studies are widespread (1120 ppm to 3360 ppm; Herrmann et al., 2004;Lowry et al., 2014;Pohl et al., 2016).Pohl et al. (2016) used a similar climate modelling setup as our study (LMDz6 forced with FOAM-derived SSTs, see Supplementary Information), so these results might be the most comparable.The simulations presented here yield CO 2 levels below the 3360 ppm threshold as suggested by Pohl et al. (2016) by 440 Ma, at the end of a 40-Myr cooling (Fig. 5a).Therefore, we demonstrate that CO 2 levels for the Ordovician are compatible with the onset of a glaciation such as the short-lived Hirnantian glaciation (445-444 Ma), even though the intense global cooling as shown by proxies cannot be replicated.For comparison, other carbon cycle models, such as GEOCARBSULFvolc, estimate CO 2 levels almost similar to the simulations presented here when applying the same degassing parametrization (variable age-zircon distribution degassing; Marcilly et al., 2021).This study estimated theoretical glacial inception levels at 1035 ppm and 2142 ppm, based on climate sensitivities of 6 • C and 3 • C. Assuming an even lower climatic sensitivity (2 • C) yields a threshold as high as 4434 ppm (Fig. 5a).
Model-data disagreement in our simulations may come, at least in part, from the limitations of our models.Regarding the climatic component, results are necessarily model-dependent.More particularly, the double intertropical convergence zone, which represents a very general bias of ocean-atmosphere general circulation models, and is amplified in our simulations by the continental configuration featuring vast ocean basins at low latitudes, reduces precipitation rates over tropical emerged lands (Figs. 2, 3).This bias alters the climate-weathering relationship.In addition, due to the large number of simulations required to force the carbon cycle component, we did not run the ocean-atmosphere version of our CMIP6-class model.Instead, the atmospheric component (LMDz6) was forced with SSTs derived from an ocean-atmosphere model of a previous generation run under identical boundary conditions (FOAM see Supplementary Information).Such setup, used in previous work (Ladant et al., 2014;Licht et al., 2014;Pohl et al., 2016), may also lead to bias in our results, although this impact remains difficult to quantify.Importantly, we emphasize that uncertainty in the model climatic sensitivity would only alter the simulated pCO 2 and not affect simulated temperatures, which constitute our focus.Another source of uncertainty could come from the computation of weathering rates.The recent developments embedded in the version of GEOCLIM used here (Park et al., 2020) largely rely on continental topography, which is used to derive local slope and physical erosion rates (see Supplementary Information).Therefore, uncertainties in topography on land, inherent in every deep-time study due to the lack of robust and widespread proxies for paleoelevation, constitute a limitation in our work.
The inability of our climate-carbon cycle model to simulate the long-term Ordovician cooling trend may also reflect missing mechanisms.In particular, including the impact of land plant colonization on weathering could help reconcile models and proxies (Porada et al., 2016;Lenton et al., 2012).However, the impact of non-vascular land plants on weathering is debated (Edwards et al., 2015;Porada et al., 2016).In addition, recent studies pinpoint the impact of the colonization scenario on weathering rates (Maffre et al., 2022), and this colonization remains poorly documented in the Ordovician.Consequently, the impact of non-vascular plants on Ordovician climate evolution remains elusive.
Another possible mechanism to explain the intensification of climate change during the Ordovician could reside in the weathering of obducted volcanic arcs as recently demonstrated by Conwell et al. (2022).Indeed, this study argues, using 87 Sr/ 86 Sr seawater records, that a shift in the weathering source lithology is synchronous with the Ordovician cooling, and this can be linked to enhanced weathering of uplifted mafic oceanic crust associated with the Taconic orogeny.This hypothesis has also been suggested by several authors (i.e., Swanson-Hysell and Macdonald, 2017); Macdonald et al., 2019).Macdonald et al. (2019) placed Laurentia about 10 • more northward between 480 and 444 Ma (Fig. 11a) compared with our reconstructions (see SI Fig S6), motivated by exhuming mafic and ultramafic lithologies during arc-continent collisions within the warm and wet tropical latitudes (±10 • ).Swanson-Hysell and Macdonald (2017) argued that a Taconic-related increase in silicate weatherability contributed to long-term Ordovician cooling, but this was not quantified in their study.One way of estimating changes in silicate weatherability is to examine the amount of exposed land where chemical weathering is assumed to be the most intense.Exposed land within ± 10 • was steadily increasing during the Ordovician (SI Table S1), peaking at around 450 Ma with ∼22 × 10 6 km 2 or ∼24% of total land exposed between 10 • S and 10 • N. In the Late Ordovician, we estimate about 3000 km globally of ophiolite-bearing sutures in our paleogeographic model or ∼5500 km in the Swanson-Hysell and Macdonald (2017) model.These estimates are based on Macdonald et al. ( 2019) but here revised since no Ordovician ophiolites are known in Greenland and all Ordovician ophiolites in Baltica (Fig. 11a) were obducted in the Silurian (e.g., Andersen and Andresen, 1994;Jakob et al., 2022).Fig. 1 in Macdonald et al. (2019) at 445 Ma, therefore, requires a major revision in suture lengths.In order to compare suture lengths with exposed land we (generously) assume an obduction width of 100 km, and as a simple theoretical exercise we multiply the resulting area by a factor of 10 to account for potentially higher weatherability of mafic/ultramafic rocks and elevated topography (Dessert et al., 2003;West et al., 2005).Note, however, that hydrothermal basaltic alteration and ultramafic serpentinization would have lowered their weathering efficiency compared to direct weathering of juvenile volcanics.This exercise shows that the potentially added weathering effect of obducted rocks in the tropics (compared to the amount of exposed land, mostly 'granitic') amounts to 27% at most, and as low as 3% considering the Scotese (2016) reconstructions (Fig. 11).
Adding the potential effect of obduction and increased weathering using the GEOCARBSULF model (e.g., increasing the f Aw / f A parameter, i.e. "fraction of land area undergoing chemical weath- ering") results in a theoretical CO 2 reduction of 150 ppm in our preferred paleogeographic model (Fig. 11b).This exercise however, assumes that the climate gradients during the Ordovician were similar to modern times, but the occurrences of evaporites in Laurentia at low latitudes suggest a dry climate, rather than warm and wet tropics (e.g., Boucot et al., 2013;Torsvik & Cocks 2016;Torsvik et al., 2021).
When integrating island arcs and increased topography linked to obducted arcs as suggested by Macdonald et al. (2019) (Fig. 6, SI Fig. S6 and SI Fig. S7), our simulations did not yield any large effect linked to topography.Moreover, considering the trend in cooling enhancement seen in the simulations with a shift from low to moderate topography, only a considerable change towards very high relief can trigger large fluctuations.

Conclusion
By using a spatially resolved climate-carbon cycle Earth system model fed with updated continental reconstructions and new estimates of volcanic degassing, we show that the Early Ordovician high temperatures can be modelled within error margins of the proxy data.However, our new constraints cannot explain the intense cooling over the Mid to Late Ordovician even if revealing a progressive enhancement in Earth's surface weatherability.The calculated, theoretical degassing scenario necessary to reproduce the proxy-derived Ordovician long-term cooling trend features degassing rates reaching 3 to 5 times that of the current value, which is not supported by geodynamical models nor geological proxies.Therefore, we interpret the inability of our up-to-date climatecarbon cycle model to reproduce Ordovician cooling as the existence of a missing parameter or a faulty influence of a present one.Accepting the veracity of the high Early Ordovician temperatures, alternative sources but also sinks of carbon must be considered.The contribution of the organic carbon sub-cycle and oxidative weathering of pyrite is currently not supported by proxy data.The impact of the colonization of land surfaces by non-vascular plants may help bring models and data into a closer agreement, but the overall impact of vegetation on climate remains difficult to quantify.Our simulations with increased topography did not intensify the Ordovician cooling, but note that the model used here holds large limitations considering the integration of relief and slopes.A more detailed analysis on this subject is necessary to draw any firm conclusions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.(a) Compilation of temperature reconstructions covering the Late Cambrian to early Silurian.These datasets are from Song et al. (2019) in red, filtered to reflect only fossils found within ±30 • latitude with error margins shaded.The stippled purple curve is from Goldberg et al. (2021) and represents the average of the least diagenetically altered samples.The purple shaded area represents the spread of data of the mid-to-least diagenetically altered data presented in Fig. 4 of Goldberg et al. (2021).Finally, the green dots are the temperatures reconstructed by Trotter et al.(2008)  with their error band in shaded green.(b) Model simulations of atmospheric CO 2 levels.The green line is simulated using the forward model COPSE(Mills et al., 2019) whereas the stippled blue curve is obtained using GEOCARBSULF with updated paleogeography and degassing parameters(Marcilly et al., 2021).The dotted black curve is fromNardin et al. (2011) using GEOCLIM with the GCM FOAM.The numbers refer to the theoretical glacial inception levels with 3 and 6 • C climate sensitivity assuming 500 ppm inception level for modern times and a fainter sun in deep time(Marcilly et al., 2021).(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. (a) Modelled Global Mean Surface Temperatures (GMST), continental mean, tropical mean and tropical (21 • N-21 • S) sea-surface temperature (SST) though time.(b) Mean continental runoff (solid line, left axis), which corresponds to precipitation minus evaporation and globally cumulated continental discharge (dashed line, right axis) corresponding to the area-integrated runoff through time.(c) Total silicate weathering fluxes predicted by GEOCLIM (relative to pre-industrial control weathering flux).Climatic variables from LMDz simulations conducted at 12 PAL pCO 2 (3360 ppmv) were used to force the weathering model.

Fig. 4 .
Fig. 4.Weathering patterns for each time slice with a constant pCO 2 of 12 PAL (3360 ppmv).Simulated weathering patterns are only driven by continental temperature and runoff, both depending on tectonic settings.B, Baltica, L, Laurentia, S, Siberia, La, Laurussia.Same projection and graticules as in Fig.3.

Fig. 5
Fig.5.a) Equilibrium pCO 2 with constant (black), and variable degassing from subduction fluxes (SBF; pink) and age-zircon distribution (Z; turquoise).The purple arrow and band display the glaciation threshold range estimated from climate models coupled with land-ice models(Herrmann et al., 2004;Lowry et al., 2014;Pohl et al., 2016).The numbers refer to the theoretical glacial inception levels with 2, 3 and 6 • C climate sensitivity assuming 500 ppm inception level from modern times and fainter sun activity in deeper time(Marcilly et al., 2021).b) Global mean surface temperature (GMST) of GEOCLIM, run with constant degassing (black) or variable degassing from SBF and Z (pink and turquoise).The temperature evolution at fixed CO 2 (12 PAL or 3360 ppmv) is added (orange) for comparison.Climatic trends shown with background colours in panel b are extracted from proxy compilations of different studies(Trotter et al., 2008;Song et al., 2019;Scotese et al., 2021;Goldberg et al., 2021) and simplified as follows: cooling until 480 Ma followed by rapid cooling until 450-440 Ma followed by a plateau.

Fig. 6 .
Fig. 6.Global mean average surface temperature (GMST) simulated in GEOCLIM.Simulations with a) constant degassing, b) variable degassing from subduction fluxes, andc) variable degassing from zircon proxy.Each panel includes 4 simulations with varying topography.The solid line is the topography used in our primary set of maps, the dashed line with doubled relief altitude, the dotted line with doubled relief altitude and increased topography fromMacdonald et al. (2019) and finally the dash-dotted line is the same as the dotted line but with increased slopes.

Fig. 7 .
Fig. 7. Degassing relative to present-day flux derived from subduction fluxes ("subduction flux deg.") from a full-plate model (pink solid line) and age zircon distribution ("Zircon deg.") (turquoise dashed line).The proxies are plotted with their estimated error range.For comparison purposes, other estimations are plotted: Mills et al. (2019) based on hybrid subduction and rift lengths, Berner (2004) using sealevel inversion and calculated subduction flux from the Merdith et al. (2021) full plate model.

Fig. 8 .
Fig. 8. Equilibrium tropical sea-surface temperature (SST) simulated in GEOCLIM.Simulations with a) constant degassing, b) variable degassing from subduction fluxes, and c) variable degassing from zircon proxy.Shaded pink areas represent standard deviation.Purple and red solid lines represent mean tropical SST estimates from Goldberg et al. (2021) and Song et al. (2019), respectively.

Fig. 9
Fig. 9. a) CO 2 degassing rate required by GEOCLIM to match the observed tropical sea-surface temperature (SST) of Goldberg et al. (2021) (purple), Song et al. (2019) (red line) and for comparison degassing from subduction fluxes proxy (pink line) and zircon proxy (turquoise line).b) pCO 2 of the corresponding GEOCLIM runs as in a.

Fig. 10 .
Fig. 10.(a) Oxygen isotope data from Veizer et al. (1999) in black.Changes following the secular trend from Veizer and Prokoph (2015) are indicated in grey.The dataset of the "warm" models from Bergmann et al. (2018) and Henkes et al. (2018) are represented with their error margins in shaded corresponding colours.(b)Reconstructed paleotemperatures from phosphatic fossils from Song et al. (2019), Bulk-rock δ 18 O from Goldberg et al. (2021) and carbonate fossil database ofVeizer et al. (1999).The grey data points represent the calculated temperature of the Veizer dataset following changes due to secular trend as suggested by the Veizer hypothesis.

Fig. 11
Fig. 11.a) Late Ordovician reconstruction of the Iapetus bordering continents.We place the southern margin of Laurentia at ∼20 • S (Model-I, Torsvik et al., 2012) whilst Swanson-Hysell and Macdonald (2017; see their figure DR1 at 445 Ma) place Laurentia ∼10 • more northward (Model-II, see supplementary section G).Thick dark brown lines in Model-II reconstruction were named 'ophiolite-bearing sutures' in Macdonald et al. (2019; see their Fig. 1 at 445 Ma) but there are no known ophiolites in Greenland and ophiolites in western Baltica were all obducted in the Silurian.This has been accounted for in Model-I reconstruction where no ophiolites in the Iapetus realm were obducted within 10 • from the equator (yellow band).Globally, reconstruction in Swanson-Hysell and Macdonald (2017) and Macdonald et al. (2019) are identical to those of Torsvik & Cocks (2016), except for locating Laurentia 10 • more northward during the Ordovician (Supplementary Fig. G).This result in a much wider Iapetus Ocean at 445 Ma, the Caroline Terrane has yet not collided with Laurentia and the Model-II requires fast latitudinal plate velocities (18 cm yr −1 ) to close the Iapetus at around 430 Ma.(b) Exposed land globally within 10 • from the equator based on our Model I, Model-II (as Model-I except modified Laurentia) and Scotese (2016; Model III).All models are based on exposed land mapped out in 10 Myrs intervals in Marcilly et al. (2021) and peaks in exposed land in all models occur at 450 Ma.Swanson-Hysell and Macdonald (2017) argued for a peak in obducted arc lengths in the Late Ordovician but did not attempt to quantify the effect in terms of atmospheric CO 2 drawdown.As a simple back-of-the-envelope calculation, we estimate the effect of increased weathering by considering an obduction width of 100 km, which have been multiplied by 10 in order to account for a higher theoretical weathering efficiency (mafic versus felsic rocks) and elevated topography.Compared with the total exposed land (10 • S-10 • N) the added weathering potential can be (generously) estimated to 27% (Model-I), 15% (Model-II) and 3% (Model II).In terms of CO 2 drawdown that amounts to 275 ppm (Model-I), 150 ppm (Model-II) and 30 ppm (Model-III) in the GEOCARBsulf model.The red stippled line is based on the original estimates in Macdonald et al. (2019).