# Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering

Yves Goddéris, Yannick Donnadieu, Sébastien Carretier, Markus Aretz, Guillaume Dera, Mélina Macouin, Vincent Regard

► To cite this version:

Yves Goddéris, Yannick Donnadieu, Sébastien Carretier, Markus Aretz, Guillaume Dera, et al.. Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering. Nature Geoscience, 2017, 10 (5), pp.382 - 386. 10.1038/NGEO2931. hal-01584138

HAL Id: hal-01584138

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

Submitted on 4 Oct 2025

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 Attribution 4.0 International License
# Onset and ending of the late Palaeozoic ice age triggered by tectonically paced rock weathering

Yves Goddéris<fnref n="1" />*†, Yannick Donnadieu²,³*†, Sébastien Carretier<fnref n="1" />, Markus Aretz<fnref n="1" />, Guillaume Dera<fnref n="1" />, Mélina Macouin<fnref n="1" /> and Vincent Regard<fnref n="1" />

The onset of the late Palaeozoic ice age about 340 million years ago has been attributed to a decrease in atmospheric CO₂ concentrations associated with expansion of land plants, as plants both enhance silicate rock weathering—which consumes CO₂—and increase the storage of organic carbon on land. However, plant expansion and carbon uptake substantially predate glaciation. Here we use climate and carbon cycle simulations to investigate the potential effects of the uplift of the equatorial Hercynian mountains and the assembly of Pangaea on the late Palaeozoic carbon cycle. In our simulations, mountain uplift during the Late Carboniferous caused an increase in physical weathering that removed the thick soil cover that had inhibited silicate weathering. The resulting increase in chemical weathering was sufficient to cause atmospheric CO₂ concentrations to fall below the levels required to initiate glaciation. During the Permian, the lowering of the mountains led to a re-establishment of thick soils, whilst the assembly of Pangaea promoted arid conditions in continental interiors that were unfavourable for silicate weathering. These changes allowed CO₂ concentrations to rise to levels sufficient to terminate the glacial event. Based on our simulations, we suggest that tectonically influenced carbon cycle changes during the late Palaeozoic were sufficient to initiate and terminate the late Palaeozoic ice age.

The late Palaeozoic ice age (LPIA) started around 340–330 million years ago (Ma)<fnref n="1" />. It is one of the most prominent glacial events in Earth's history, now seen as a dynamic succession of ice advances and retreats over Gondwana<fnref n="1" />⁻³. Continental ice sheets are hypothesized to have persisted until early Permian times. This late Palaeozoic cooling has been attributed to a CO₂ decline driven by the colonization of the continents by vascular plants. Land plants are thought to have increased continental rock weatherability (defined as the product of the factors that affect rock chemical weathering other than climate change⁴), which, in turn, should have promoted CO₂ removal⁵. In numerical models of the global carbon cycle, the effect of land plants on weathering is simulated by introducing a factor which increases by a factor of six the global silicate rock weatherability between 380 and 350 Ma (ref. 5). In addition, the production of lignin in the postulated absence of decomposing fungi may have lead to an efficient storage of organic carbon, witnessed by an increase in the ¹³C content of marine carbonates. This latter effect has been accounted for in global carbon models through the inversion of the δ¹³C seawater composition⁵,⁶. Consequently, previous CO₂ models for the late Palaeozoic are primarily dominated by biotic forcing.

Although commonly accepted, this biotic hypothesis for the LPIA must be re-evaluated. First, the colonization of the continents by vascular plants was achieved in the Late Devonian, several tens of million years before the glacial event⁷. Second, the hypothesis of the absence of lignin decomposers has been invalidated recently: lignin-degrading fungi may have been present before the Carboniferous⁸. Third, an increase in terrestrial organic carbon burial centred on the LPIA is observed⁹, but it is poorly correlated with the ¹³C/¹²C signal, which displays uniform high values from 360 to 260 Ma (Fig. 1), and with the timing of the colonization by plants. The biotic hypothesis has left aside the role of the tectonic context (continental

drift and mountain uplift) on the late Palaeozoic climatic evolution. This period is characterized by an important geodynamic event that reshaped the Earth's surface in the equatorial area: the uplift of the Hercynian range.

Physical erosion is one of the factors controlling the consumption of atmospheric CO₂ through silicate weathering. In the present-day large flat tropical areas, low physical erosion induces the development of thick saproliths (defined as the surficial ‘soil’ layer resulting from the progressive chemical transformation of the bedrock, and extending from the surface to the impervious bedrock). Those thick saproliths maintain the system in a supply-limited weathering regime characterized by low CO₂ consumption¹⁰. The uplift of the Hercynian range may have had a huge impact on the CO₂ sink by increasing physical erosion and removing the thick saprolith through the development of steep slopes. Various mathematical formulations have been proposed to describe the tight link between erosion and chemical weathering at various spatial scales: from local vertical profiles¹¹,¹² to the regional scale¹³. However, although this suite of processes has been intensively discussed in the Cenozoic context, no long-term coupled geochemical-climate model accounting for these formalisms has been applied to the global scale.

In this paper, we quantify the contribution of the large Hercynian orogeny straddling the equator to the long-term CO₂ evolution from 350 Ma to 272 Ma. Five time slices spanning the period have been investigated. For each time slice, three sets of simulations have been first performed, one neglecting the role of saprolith thickness and the two others that account for it using two different formulations.

## Model description and calibration

We first use the spatially resolved GEOCLIM model¹⁴,¹⁵. In this case, the specific flux of CO₂ consumption by silicate weathering

<footnote n="1">¹Géosciences Environnement Toulouse, CNRS—Université Paul Sabatier - IRD, 31400 Toulouse, France. ²Laboratoire des Sciences du Climat et de l'Environnement, LSCE-IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France. ³Aix-Marseille Univ, CNRS, IRD, CEREGE, 13545 Aix-en-Provence, France. †These authors contributed equally to this work. *e-mail: yves.godderis@get.omp.eu; donnadieu@cerege.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1** | Time evolution of the isotopic composition of seawater, atmospheric CO₂ levels and the organic carbon burial over the late Palaeozoic. **a**, Seawater ⁸⁷Sr/⁸⁶Sr. The grey shaded area stands for the ±1σ envelope around the 10 Myr moving average of the data³³. Dots and diamonds stand for model simulations: GEOCLIM (no shielding effect on weathering); GEOCLIM-REG accounting for the shielding effect; GEOCLIM-REG (plateau) assuming that the Hercynian orogeny is a high-altitude plateau; GEOCLIM-WEST using the West model to calculate weathering; GEOCLIM-SENS accounting for the shielding effect and a doubled sensitivity of climate to CO₂. All simulation results coincide at 350 Ma. **b**, Calculated time evolution of atmospheric CO₂ level (dots and diamonds). Green horizontal lines represent the CO₂ threshold below which a continental-scale glaciation is initiated²⁵. The blue vertical bars represent the major pulses of the late Palaeozoic ice age¹,⁹. The grey shaded area is the ±1σ envelope around the moving average of the proxy data (10-Myr time window)³⁴. Crosses stand for most recent CO₂ estimates³. All simulation results coincide at 350 Ma. **c**, The burial flux of organic carbon-rich sediments (compilation for North America⁸). **d**, Seawater δ¹³C measured on sedimentary carbonates (crosses)³³. The solid red line is a moving average of the data (window width of 10 Myr).</figcaption></figure>

(mol m⁻² yr⁻²) for each continental grid cell *i* is calculated according to the following common equation¹⁶:

$$W_i^{\text{CO}_2} = k_G^{\text{CO}_2} R_i \exp \left[ \frac{-E_a}{\rho} \left( \frac{1}{T_i} - \frac{1}{288.15} \right) \right] \quad (1)$$

where $R_i$ (cm yr⁻¹) and $T_i$ (K) are, respectively, the mean annual runoff and the mean annual surface temperature for the grid element *i*. $E_a$ (J mol⁻¹) is an apparent activation energy and $\rho$ (J mol⁻¹ k⁻¹) the ideal gas constant. $k_G^{\text{CO}_2}$ is a calibration constant. For each time slice, we use the climate model FOAM¹⁷ to calculate $R_i$ and $T_i$ (see Methods and Supplementary Information). Equation (1) predicts maximum weathering in the equatorial warm and humid region.

In our second case, the weathering module of the GEOCLIM numerical model is improved by explicitly introducing the effect of the saprolith thickness on the efficiency of weathering (GEOCLIM-REG). Indeed, despite warm temperatures characterizing low latitudes, continental weathering is not optimal there, because of the development of thick weathering profiles, limiting the downward circulation of water and shifting the system from a

weathering-limited regime (where the kinetics of the dissolution reactions matter) to a supply-limited regime (where the accessibility of fresh minerals to the flowing waters becomes the limiting factor)¹⁸,¹⁹. Very thin soils as found in rapidly eroding areas are also not optimal because the short residence time of the water limits the weathering processes. To account for these processes, we introduce a shielding factor in equation (1) (ref. 13):

$$W_i^{\text{CO}_2} = k_{Gr}^{\text{CO}_2} R_i \exp \left[ \frac{-E_a}{\rho} \left( \frac{1}{T_i} - \frac{1}{288.15} \right) \right] \times \left[ \exp \left( \frac{-H_i}{d_1} \right) - k_1 \exp \left( \frac{-H_i}{d_2} \right) \right] \quad (2)$$

where $H_i$ is the thickness of the weathering profile from the surface to the bedrock (m) for the grid element *i*. $k_1$, $d_1$ and $d_2$ are calibration constants (see Methods and Supplementary Information). The shielding factor introduces a dependence of weathering²⁰ on saprolith thickness (Supplementary Fig. 1). On each grid cell, saprolith is produced at depth by the chemical attack of the bedrock (equation (2)) and removed by physical erosion at the surface (which depends on the water discharge and local slope, see Methods and Supplementary Fig. 2). In tropical flat environments, this shielding factor decreases silicate weathering by a factor of 10 to 1,000 (Supplementary Fig. 3). In a latter case called GEOCLIM-WEST, another mathematical formulation²¹ allowing for the simulation of the shielding effect is used to test the model dependency of our results.

Prior to any deep time simulation each mathematical formulation of weathering fluxes is tested with modern boundary conditions. The GEOCLIM-REG and GEOCLIM-WEST models are able to reproduce the CO₂ consumption by silicate weathering for a variety of continental watersheds. Conversely, the GEOCLIM model largely overestimates weathering fluxes in tropical areas, demonstrating the importance of the shielding factor (Supplementary Figs 4 and 5 and Methods).

## Atmospheric CO₂ in the late Palaeozoic

Using GEOCLIM and versions including the shielding effect, GEOCLIM-REG and GEOCLIM-WEST, atmospheric CO₂ levels are simulated (Fig. 1b) for five time slices spanning the late Palaeozoic²² (see Methods). For each palaeogeography (Supplementary Fig. 6), the calculation of the atmospheric CO₂ level arises from the equality between the sink—that is, the spatial distribution of continental silicate weathering—and the source—that is, the Earth degassing flux (see Methods). In both sets of simulations, we fix the degassing rate such that the calculated CO₂ at 350 Ma (prior to the glaciation) falls in the middle of the CO₂ proxy envelope²³, which is slightly above four times the pre-industrial CO₂ level. Sensitivity tests to the starting CO₂ value for 350 Ma have been run and show that results presented below remain valid within CO₂ proxy uncertainties (see Supplementary Information).

Neglecting the shielding effect, GEOCLIM simulates a significant increase in atmospheric CO₂ levels from 1,000 to 3,000 ppm between 350 Ma and 272 Ma. This is the consequence of the development of large arid areas resulting from the assembly of the supercontinent Pangaea²⁴. The formation of the supercontinent modifies the source of moisture advection from the ocean to the continents. Rainout of air masses from coastal to inland regions results in very dry continental interiors. Less water flowing over the continents implies a decrease of the carbon sink acting through the silicate weathering. This aridity is also due to latitudinal drift: at 350 Ma, there were only 18 × 10⁶ km² of emerged continents within arid subtropical divergence zones (16% of the global continental area, compared to 30% today). The continental surface in divergence areas then rose to reach 31% of the continental surface at 290 Ma.

Atmospheric CO₂ simulated by GEOCLIM-REG (and-WEST) differs from GEOCLIM not only in terms of absolute values but also
in terms of temporal trend (Fig. 1b). It now shows a large decrease, crossing the atmospheric CO₂ threshold defined by climate ice-sheet modelling²⁵, between 350 and 323 Ma, whatever the shielding model used. Atmospheric CO₂ remains low across the time window corresponding to most of the late Palaeozoic ice age. Then, CO₂ concentrations rise rapidly at 290 Ma (Sakmarian), and converge to values similar to those simulated with GEOCLIM. To illustrate how physical erosion affects our results, we focus on simulations performed at an atmospheric CO₂ concentration of 1,120 ppm using results from GEOCLIM-REG (same behaviour occurs with GEOCLIM-WEST, see Supplementary Fig. 7). The contribution of the Hercynian orogen to the global CO₂ consumption by silicate weathering fluctuates through time (Fig. 2). It starts from 0% at 350 Ma and culminates to 39% at 308 Ma. It rapidly decreases in the Permian to end at about 19% at 272 Ma. Conversely, the contribution of continental areas located outside the Hercynian orogenic domain continuously decreases. This is due to the increase in aridity during the Pangaea assembly, which explains the trend of increasing CO₂ found in GEOCLIM. Figure 2b displays the contribution of the silicate weathering in the Hercynian domain using simulated steady-state CO₂ levels (GEOCLIM-REG simulation, black diamonds in Fig. 1b). In this case, the global silicate weathering is constant through time, balancing the solid Earth degassing flux, but the Hercynian contribution fluctuates from one time slice to the other.

The calculated values for atmospheric CO₂ are in agreement with the reconstructions based on proxy data from 350 to 308 Ma. In the Early Permian, modelled CO₂ levels reach 6,000 to 13,000 ppmv, above the 1σ-envelope of the proxy data (Fig. 1a), which culminates at about 2,500 ppmv. However, proxy data saturate at high CO₂, and their reliability decreases with CO₂ level (see Supplementary Information). Furthermore, model value depends on the sensitivity of the climate model to CO₂. The climate model used in this study displays a low sensitivity of 2 °C global warming for a CO₂ doubling. We performed an additional simulation (GEOCLIM-SENS) in which we doubled the climate sensitivity to CO₂ (see Supplementary Information). This simulation still predicts calculated CO₂ below the glacial threshold at 323 and 308 Ma, whereas Permian calculated values are now ranging from 2,600 to 4,000 ppmv, closer to proxy estimates. Beyond the discussion about absolute values, the model output mimics the overall trend defined by the proxy data: CO₂ decrease after 350 Ma, low CO₂ value during the LPIA, and CO₂ increase in the Permian.

The Tournaisian (350 Ma) weathering pattern differs widely from one set of simulations to the other (Fig. 3). In the GEOCLIM simulation, weathering fluxes are maximal at equatorial latitudes. In the GEOCLIM-REG simulation, the development of thick saprolith profiles, especially in the flat area along the equator, limits the contribution of this zone to the global weathering flux (Fig. 3). At 308 Ma, the uplift of the Hercynian mountain range located along the equator results in a very high runoff, in excess of 4 m yr⁻¹ at several places. Together with the increase in mean slope, these climatic conditions favour the onset of physical erosion and induce a 20% increase in the total CO₂ consumption flux (Fig. 2) despite the global trend towards an arid world seen before. Our simulation shows that restricted continental areas such as the Hercynian range may become the major climatic driver if optimal conditions occur there (tectonically active and high runoff). Once the range collapses, runoff and erosion along the equator decrease (272 Ma), thicker saproliths start to accumulate, and the weathering in this area can no longer compensate for the global decrease in runoff.

## Coeval evolution of the seawater strontium isotopic ratio

Accounting explicitly for physical erosion, the GEOCLIM-REG simulation predicts that the atmospheric CO₂ level falls below the

<figure><img src="image_3.png" /><figcaption>**Figure 2** | Calculated contribution of the weathering of the Hercynian mountains to the global weathering flux. **a**, Atmospheric CO₂ fixed at four times the pre-industrial level. **b**, Atmospheric CO₂ using simulated steady-state values (GEOCLIM-REG, black diamonds of Fig. 1b). Green bars stand for the Hercynian contribution and blue bars for the whole continental area excluding the Hercynian orogeny. All fluxes are normalized to the global weathering flux at 350 Ma (in %).</figcaption></figure>

glacial threshold just after 350 Ma until about 300 Ma. This is in good agreement with the timing of major phases of the LPIA¹. Accounting for the shielding effect in weathering equations also has implications concerning our understanding of the Sr isotopic cycle (see Supplementary Information for a full description). ⁸⁷Sr/⁸⁶Sr ratios measured in marine carbonates have often been interpreted in terms of fluctuation in the relative importance of the radiogenic continental Sr versus unradiogenic mantle Sr input. Yet, because the budget is calculated with a background CO₂ at steady state assuming a constant degassing rate, the continental versus oceanic strontium fluxes ratio should be constant through time²⁶. Another way to explain changes in seawater ⁸⁷Sr/⁸⁶Sr is to change the isotopic composition of the river discharge into the oceans²⁷. Here, the development of steep slopes in the Hercynian range promotes the weathering of crustal rocks (high ⁸⁷Sr/⁸⁶Sr ratio), while the coeval drop in temperature due to the decrease in CO₂ inhibits the weathering of young volcanic rocks²⁸ (low ⁸⁷Sr/⁸⁶Sr ratio). Both processes result in a temporary rise in the riverine seawater Sr isotopic ratio synchronous with the apex of the uplift phase, in agreement with data (Fig. 1a).

Among potential uncertainties in our results, the timing of the orogeny and the shape of the Hercynian range are the most important. Regarding the timing, we have defined a scenario relying on several published studies (see Supplementary Information). Regarding the shape of the range, our first case uses a Himalayan type of range, with high altitudes and steep slope (around 20%). We also tested the case where the central part of the Hercynian orogeny is a high-altitude plateau, with reduced slopes (around 10%, see Supplementary Information, GEOCLIM-REG plateau simulation). The overall evolution of the weathering fluxes is not affected: the CO₂ consumption by silicate weathering rises from 350 Ma to 308 Ma, and then declines. But the amplitude of the changes in atmospheric
<figure><img src="image_4.png" /><figcaption>b GEOCLIM_REG</figcaption></figure>

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

**Figure 3** | Palaeogeographic pattern of the CO₂ consumption by silicate weathering. Log₁₀ of the CO₂ consumption (Gmol yr⁻¹) for three of the five simulated time slices (atmospheric CO₂ fixed at four times the pre-industrial level). **a**, The three globes display the results of the GEOCLIM simulation (no saprolith shielding effect). Maximum weathering occurs in warm and humid areas. **b**, When the shielding effect is accounted for (GEOCLIM-REG), weathering is inhibited by the development of thick saprolith along the equator, except when steep slopes and high runoff maintain a high erosive flux (308 Ma). The grey continental surfaces represent arid zones with no runoff and no weathering.

CO₂ is slightly reduced in the plateau configuration (Fig. 1b), as the highest elevations reached during the maximal phase of the orogeny are no longer accompanied by a change in the mean slope which saturates at 10%.

**An integrated scenario for the late Palaeozoic glaciation**

The classical hypothesis linking the LPIA to the colonization of the continents by land plants must be revisited. Modern environments characterized by warm and humid climatic conditions favouring a dense vegetation cover display high weathering rates and associated CO₂ consumption only if physical erosion periodically removes the thick soil cover²⁹,³⁰. Moreover, the transport of the terrestrial organic carbon by rivers and its subsequent burial in sediments are both positively correlated to physical erosion rates of continents³¹,³², more physical erosion favouring the transfer of organic carbon to the ocean and its burial. This suggests that the modification of erosion/weathering systems by the Hercynian range uplift may have had a major control on the pumping of CO₂ through organic carbon burial (Fig. 1c). The tectonic scenario is not competing with the ‘biotic’ hypothesis: our results support a holistic picture of the carbon cycle where physical erosion promotes silicate weathering and organic carbon burial.

In summary, our model demonstrates that both horizontal (continental drift) and vertical (mountain uplift) tectonics must be coupled for a model to accurately simulate the Earth's climatic and geochemical history. If the rise of the Hercynian equatorial mountain range may plunge the Earth into a glacial stage, the assembly of a supercontinent may end the glaciation. Going one step further, our study shows that expansive tropical lowlands, and not only high altitude mountains, are important drivers of the carbon

cycle and have to be accounted for in any attempts aiming to reconstruct long-term climatic trend.

**Methods**

Methods, including statements of data availability and any associated accession codes and references, are available in the online version of this paper.

Received 30 November 2016; accepted 10 March 2017; published online 10 April 2017

**References**

1. Montanez, I. P. & Poulsen, C. J. The Late Paleozoic ice age: an evolving paradigm. *Annu. Rev. Earth Planet. Sci.* **41**, 629–656 (2013).

2. Isbell, J. L. et al. Glacial paradoxes during the late Paleozoic ice age: evaluating the equilibrium line altitude as a control on glaciation. *Gondwana Res.* **22**, 1–19 (2012).

3. Montanez, I. P. et al. Climate, $p_{CO_2}$, and terrestrial carbon cycle linkages during late Paleozoic glacial–interglacial cycles. *Nat. Geosci.* **9**, 824–828 (2016).

4. Kump, L. R. & Arthur, M. A. in *Tectonic Uplift and Climate Change* (ed. Ruddiman, W. F.) 399–426 (Springer, 1997).

5. Berner, R. A. *The Phanerozoic Carbon Cycle* (Oxford Univ. Press, 2004).

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

7. Davies, N. S. & Gibling, M. The sedimentary record of Carboniferous river: continuing influence of land plant evolution on alluvial process and Palaeozoic ecosystems. *Earth Sci. Rev.* **120**, 40–79 (2013).

8. Nelsen, M. P., Dimichele, W. A., Peters, S. E. & Boyce, C. K. Delayed fungal evolution did not cause the Paleozoic peak in coal production. *Proc. Natl Acad. Sci. USA* **9**, 2442–2447 (2016).

9. Montanez, I. P. A late Paleozoic climate window of opportunity. *Proc. Natl Acad. Sci. USA* **113**, 2334–2336 (2016).
10. Stallard, R. F. Relating chemical and physical erosion. *Rev. Mineral.* **31**, 543–564 (1995).

11. Gabbet, E. J. & Mudd, S. M. A theoretical model coupling chemical weathering rates with denudation rates. *Geology* **37**, 151–154 (2009).

12. Hilley, G. E., Chamberlain, C. P., Moon, S., Porder, S. & Willet, S. D. Competition between erosion and reaction kinetics in controlling silicate-weathering rates. *Earth Planet. Sci. Lett.* **293**, 191–199 (2010).

13. Carretier, S., Godderis, Y., Delannoy, T. & Rouby, D. Mean bedrock-to-saprolite conversion and erosion rates during mountain growth and decline. *Geomorphology* **209**, 39–52 (2014).

14. Donnadieu, Y. et al. A GEOCLIM simulation of climatic and biogeochemical consequences of Pangea breakup. *Geochem. Geophys. Geosyst.* **7**, Q11019 (2006).

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

16. Oliva, P., Viers, J. & Dupré, B. Chemical weathering in granitic crystalline environments. *Chem. Geol.* **202**, 225–256 (2003).

17. Donnadieu, Y., Puceat, E., Moiroud, M., Guillocheau, F. & Deconinck, J. F. A better-ventilated ocean triggered by Late Cretaceous changes in continental configuration. *Nat. Commun.* **7**, 10316 (2016).

18. Godderis, Y., Donnadieu, Y., Tombozafy, M. & Dessert, C. Shield effect on continental weathering: implication for climatic evolution of the Earth at the geological timescale. *Geoderma* **145**, 439–448 (2008).

19. Stallard, R. F. & Edmond, J. M. Geochemistry of the Amazon: 3. Weathering chemistry and limits to dissolved inputs. *J. Geophys. Res.* **92**, 8293–8302 (1987).

20. Strudley, M., Murray, A. & Haff, P. Emergence of pediments, tors, and piedmont junction from a bedrock weathering-regolith thickness feedback. *Geology* **34**, 805–808 (2006).

21. West, J. A. 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).

22. Golonka, J. in *Phanerozoic Reef Patterns* Vol. 72 (eds Kiessling, W., Flügel, E. & Golonka, J.) 21–75 (SEPM Special Publications, 2002).

23. Royer, D. L. in *The Atmosphere History* Vol. 6 (ed. Farquhar, J.) 251–267 (Elsevier, 2014).

24. Tabor, N. J. & Poulsen, C. J. Paleoclimate across the Late Pennsylvanian–Early Permian tropical paleolatitudes: a review of climate indicators, their distribution, and relation to paleophysiographic climate factors. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **268**, 293–310 (2008).

25. Lowry, D. P., Poulsen, C. J., Horton, D. E., Torsvik, T. & Pollard, D. Thresholds for Paleozoic ice sheet initiation. *Geology* **42**, 627–630 (2014).

26. Goddéris, Y. & François, L. M. The Cenozoic evolution of the strontium and carbon cycles: relative importance of continental erosion and mantle exchanges. *Chem. Geol.* **126**, 169–190 (1995).

27. Galy, A., France-Lanord, C. & Derry, L. A. The strontium isotopic budget of Himalayan rivers in Nepal and Bangladesh. *Geochem. Cosmochim. Acta* **63**, 1905–1925 (2002).

28. Li, G. & Elderfield, H. Evolution of carbon cycle over the past 100 million years. *Geochem. Cosmochim. Acta* **103**, 11–25 (2013).

29. Millot, R., Gaillardet, J., Dupré, B. & Allègre, C. J. The global control of silicate weathering rates and the coupling with physical erosion: new insights from rivers of the Canadian Shield. *Earth Planet. Sci. Lett.* **196**, 83–98 (2002).

30. West, A., Galy, A. & Bickle, M. Tectonic and climatic controls on silicate weathering. *Earth Planet. Sci. Lett.* **235**, 211–228 (2005).

31. Galy, V. et al. Efficient organic carbon burial in the Bengal fan sustained by the Himalayan erosional system. *Nature* **450**, 407–410 (2007).

32. Galy, V., Peucker-Ehrenbrink, B. & Eglinton, T. Global carbon export from the terrestrial biosphere controlled by erosion. *Nature* **521**, 204–207 (2015).

33. Veizer, J. et al. $^{87}$Sr/$^{86}$Sr, $\delta^{13}$C and $\delta^{18}$O evolution of Phanerozoic seawater. *Chem. Geol.* **161**, 59–88 (1999).

34. Royer, D. in *Treatise on Geochemistry* Vol. 6 (eds Holland, H. D. & Turekian, K. K.) 251–267 (Elsevier, 2014).

## Acknowledgements

CNRS INSU/SYSTER program and the ANR projects TERRES and Anox-Sea have provided funding for this work. We thank the CEA/CCRT for providing access to the HPC resources of TGCC under the allocation 2015-012212 made by GENCI.

## Author contributions

Y.G. and Y.D. conceived the study and wrote the paper with contributions from all co-authors. The physical erosion model was developed by Y.G., S.C. and Y.D.; Y.G. and Y.D. performed the simulations. G.D., M.A., M.M. and V.R. reconstructed the palaeogeographic settings of the late Palaeozoic (five maps).

## Additional information

Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Correspondence and requests for materials should be addressed to Y.G. or Y.D.

## Competing financial interests

The authors declare no competing financial interests.
## Methods

**Modelling principle.** The residence time of carbon in the exospheric system (ocean+atmosphere) is about 200 kyr (ref. 35). This short residence time allows one to describe the geologic carbon cycle evolution as a succession of quasi steady states. Given that simulated time slices are separated by 20 to 30 million years, we write the following equation³⁶:

$$F_{\text{vol}} = W + \Delta\text{org} \quad (3)$$

where $F_{\text{vol}}$ is the total CO₂ outgassing of the solid Earth (mol yr⁻¹), W the total consumption of CO₂ by continental rock weathering (mol yr⁻¹), and Δorg is the net consumption of CO₂ by the geologic organic carbon cycle (organic carbon burial minus oxidation of sedimentary organic carbon, mol yr⁻¹). Equation (3) is the consequence of the short residence time of carbon compared to the geologic timescale, and of the existence of a negative feedback between climate and silicate weathering³⁷⁻³⁹. Although previous arguments for the existence of this feedback were mostly theoretical, a recent study by Zeebe and Caldeira⁴⁰ provides the first observational evidence for the operation and efficacy of this feedback over the past 610,000 years. Δorg is a source of oxygen for the ocean-atmosphere system when positive, and a net sink when negative. The residence time of oxygen in the exospheric system is about 3 million years. Working with a time step at least six times longer, we set Δorg to zero, assuming that the oxygen level is at steady state for each time slice. As we aim to reconstruct the long-term steady state of the exospheric system, the atmospheric CO₂ level is fully determined by the balance between outgassing and silicate weathering. The O₂ level adapts so that Δorg is equal to 0. Equation (3) becomes:

$$F_{\text{vol}} = W \quad (4)$$

where W is a complex function of climate (temperature and runoff), and hence of CO₂. We use a three-dimensional climate model to simulate the continental temperature and runoff patterns for each continental configuration, as a function of the atmospheric CO₂ level. These climatic patterns are used to build the mathematical function W(CO₂). Equation (4) can thus be solved for the multi-million years CO₂ through a simple numerical bisection method.

**Climate model and experiments.** Climatic maps are first generated using FOAM¹⁷,⁴¹. This model combines a low spectral resolution R15 (48 × 40 grid) atmosphere model counting 18 altimetric levels with a highly efficient medium-resolution (128 × 128 grid) ocean module composed of 24 bathymetric levels. FOAM successfully simulates many aspects of the present-day climate and compares well with other contemporary medium-resolution climate models. It was previously used to investigate numerous past climate changes, ranging from Neoproterozoic glaciations¹² to the onset of the ACC at the Eocene–Oligocene boundary⁴³.

The five sets of simulations are based on existing palaeogeographic reconstructions²² for the Tournaisian (350 Ma), the Mississippian–Pennsylvanian boundary (323 Ma), the Moscovian (308 Ma), the end of the Sakmarian (290 Ma), and the Cisuralian–Guadalupian boundary (272 Ma) (Supplementary Fig. 6 and Supplementary Information). For the sake of simplicity, simulations are referred to in the main text by age. Shared boundary conditions for all continental configurations are that surface types are set to average model surface characteristics by latitudinal bands (that is, tropical forest for the equatorial band, arid area for the tropical bands, deciduous forest for the mid-latitude and so on) and that the Earth's orbit around the Sun is circular (eccentricity = 0) and the Earth's obliquity is 23.5° (this setting leads to an equal annual insolation for both hemispheres). Solar luminosity is assumed to evolve through time according to the stellar evolution models from 3.5% reduction at 350 Ma to 2.7% at 272 Ma (ref. 44). For each continental configuration, climate simulations have been done using a large range of atmospheric CO₂ levels from 2 to 64 times the pre-industrial concentration (280 ppm).

The climatic experiments were integrated for 2,000 years without flux corrections or deep ocean acceleration. During the last 100 years of model integration, there is no apparent drift in the upper ocean (between the surface and 300 m depth) and <0.1 °C per century change in globally averaged ocean temperature. All model results have been averaged over the last 50 model years.

**CO₂ calculations.** Using climate simulations, the equation (4) can be rewritten:

$$F_{\text{vol}} = \sum_{i=1}^{\text{ngrid}} W_i^{\text{CO}_2} \quad (5)$$

where ngrid is the total number of continental grid elements for each time slice, and $W_i^{\text{CO}_2}$ is the consumption of CO₂ by silicate weathering (mol yr⁻¹) for the grid element i. Equation (5) can thus be solved for the multi-million year CO₂. Below, we show the three mathematical formulations used for the term $W_i^{\text{CO}_2}$. The Supplementary Table 1 gives numerical values for parameters defined below.

**GEOCLIM case.** The following equation is solved for each time slice:

$$F_{\text{vol}} = \sum_{i=1}^{\text{ngrid}} k_{\text{G}}^{\text{CO}_2} \text{area}_i R_i \exp \left[ \frac{-E_a}{\rho} \left( \frac{1}{T_i} - \frac{1}{288.15} \right) \right] + \sum_{i=1}^{\text{ngrid}} k_{\text{G,b}}^{\text{CO}_2} \text{area}_i R_i \exp \left[ \frac{-E_{ab}}{\rho} \left( \frac{1}{T_i} - \frac{1}{288.15} \right) \right] \quad (6)$$

where the first sum represents the consumption of CO₂ by silicate weathering excluding weathering of young volcanics, and the second sum stands for weathering of young volcanic rocks. $k_{\text{G}}^{\text{CO}_2}$ and $k_{\text{G,b}}^{\text{CO}_2}$ are constants (mol m⁻³) calibrated so that the total CO₂ consumption by Ca- and Mg-silicate weathering matches 2.5 × 10¹² mol yr⁻¹ (ref. 45), with 14% due to the weathering of young volcanic rocks⁴⁶ and area_i is the surface of the continental grid cell i (m²). $T_i$ (K) and $R_i$ (m yr⁻¹) are the temperature and runoff field from the climate model simulations. $F_{\text{vol}}$ is calculated using the climate simulated for 350 Ma at four times the pre-industrial CO₂ level, just above the glacial threshold²⁵, in the middle range of CO₂ proxies. Sensitivity experiments to this starting CO₂ condition can be found in the Supplementary Information. The value of $F_{\text{vol}}$ is then kept constant for the other time periods. Then, using equation (6), the model looks for the atmospheric CO₂ level allowing equality between $F_{\text{vol}}$ and the right term of the equation (6). This procedure requires the interpolation of the climatic fields produced by the climate model for each time slice and for several fixed CO₂ levels, so that $T_i$ and $R_i$ are continuous functions of CO₂.

**GEOCLIM-REG case.** Time evolution of the saprolith thickness $H_i$ can be written as follows:

$$\frac{dH_i}{dt} = W_i - Q_i \quad (7)$$

where $W_i$ is the net production of saprolith by weathering of the bedrock (m yr⁻¹) and $Q_i$ the removal of saprolith by physical erosion (m yr⁻¹) for the grid element i. $W_i$ is the difference between production of saprolith by the chemical attack of the bedrock and the loss of saprolith as dissolved materials. Both fluxes are assumed proportional since weathering inside the saprolith is pending at first order on the supply of materials (see the validation section below). Our simulations explore the averaged state of the continental surfaces over periods spanning at least 10 million years. Consequently, we set the derivative term in equation (7) to 0 and get:

$$W_i = Q_i \quad (8)$$

Then, we use the following relationship¹⁷ where physical erosion is set proportional to the square root of the water discharge of the grid element i and to the slope $S_i$ (%):

$$Q_i = k_a (R_i \text{area}_i)^{0.5} S_i \quad (9)$$

Using equation (2) for $W_i$, equation (8) is rewritten as follows:

$$k_d R_i t \exp \left[ \frac{-E_a}{\rho} \left( \frac{1}{T_i} - \frac{1}{288.15} \right) \right] \left[ \exp \left( \frac{-H_i}{d_1} \right) - k_1 \exp \left( \frac{-H_i}{d_2} \right) \right] = k_a (R_i \text{area}_i)^{0.5} S_i \quad (10)$$

$W_i$ (m yr⁻¹) is assumed to be dependent on climate ($R_i$ and $T_i$) and on the saprolith shielding factor defined as:

$$\text{shielding factor} = \exp \left( \frac{-H_i}{d_1} \right) - k_1 \exp \left( \frac{-H_i}{d_2} \right)$$

Parameters for this factor are defined in the Supplementary Information. This saprolith shielding factor first rises with saprolith thickness $H_i$ until it reaches a maximum value (Supplementary Fig. 1). Then, when $H_i$ (m) increases above the optimum thickness (that is, 20 cm), the saprolith shielding factor rapidly decreases¹³,⁴⁸. Called the 'humped' law, this factor simulates the decrease in saprolith production when saprolith thickness is either too thin (residence time of water in contact with the bedrock is too low) or too thick (limiting the percolation of water). Note that for a saprolith thickness equal to 0, corresponding to bedrock landscapes, the saprolith production does not fall to zero and the weatherability is equal to one third of its optimum value. First proposed by Gilbert⁴⁹ in 1877 for local soil production, the mathematical expression used here to simulate chemical weathering (including the temperature and runoff effect) is similar to the mathematical law used by others²⁰. Gabet and Mudd⁵¹ and Strudley et al.²⁰ proposed a similar theoretical description of saprolith production rate by chemical dissolution. Equation (10) is solved for the steady-state saprolith thickness $H_i$ for each time slice, for each CO₂ level and for each grid element. The calibration constant $k_a$ is fixed so that the global delivery of suspended solids to the ocean⁵² reaches 20 Gt yr⁻¹ for the present-day control simulation, assuming a mean density
of the materials of 2,800 kg m⁻³. At first order, the CO₂ consumption by silicate weathering is assumed to be proportional to the production of dissolved solids, and thus to Wᵢ, representing the net supply of materials available for weathering. Consequently, the CO₂ balance for each time slice can be written as:

$$F_{\text{vol}} = \sum_{i=1}^{\text{ngrid}} \left\{ k_{\text{Gr}}^{\text{CO}_2} \text{area}_i R_i \exp\left[\frac{-E_a}{\rho}\left(\frac{1}{T_i} - \frac{1}{288.15}\right)\right] \times \left[\exp\left(\frac{-H_i}{T_i R_i S_i}\right) - k_1 \exp\left(\frac{-H_i}{T_i R_i S_i}\right)\right] + W_i^{\text{volc}} \right\} \quad (11)$$

where Hᵢ is a function of Rᵢ, Tᵢ and Sᵢ. kᵢ CO₂ is the calibration constant (mol m⁻³). Wᵢ volc is the consumption of CO₂ by weathering of young volcanic rocks depending on Rᵢ and Tᵢ (see GEOCLIM case). We do not apply the 'saprolith correction' to young volcanic rocks. This is based on field observations showing that basalt weathering is well described by Arrhenius' Law, which, in turn, suggests a minor role for saprolith thickness⁴⁶. Since Sᵢ is prescribed (see Supplementary Information), and since the climatic patterns Tᵢ and Rᵢ depends on CO₂, the above equation can be solved for CO₂ for each simulated time slice (Supplementary Fig. 14). Fᵢ vol is still fixed so that CO₂ is just above the glacial threshold at 350 Ma. Then, Fᵢ vol is held constant at this value for the other simulated periods.

**GEOCLIM-WEST case.** West²¹ proposes a parametric relationship between chemical denudation and total (physical+chemical) denudation. Assuming that CO₂ consumption is proportional to chemical denudation (proportionality constant k Gr CO₂ ), the West law²¹ becomes:

$$k_{\text{Grw}}^{\text{CO}_2} W_i^{\text{CO}_2} = (k_{\text{Grw}}^{\text{CO}_2} W_i^{\text{CO}_2} + Q_i) \times \left[ 1 - \exp \left( -K \left( 1 - \exp(-k_w R_i) \right) \exp \left( \frac{E_a}{RT_0} - \frac{E_a}{RT_i} \right) \frac{\left( \frac{z}{W_i^{\text{CO}_2} + Q_i} \right)^{\sigma+1}}{\sigma+1} \right) \right] \quad (12)$$

As in the GEOCLIM-REG simulation, the physical erosion Qᵢ is calculated for each grid element, for each continental configuration and for each CO₂ level (equation (9), right endmember). Equation (12) can then be solved for W i CO₂ . The existence of the negative feedback between climate and silicate weathering implies a close match between global silicate weathering and outgassing:

$$F_{\text{vol}} = \sum_{i=1}^{\text{ngrid}} (W_i^{\text{CO}_2} \{T_i, R_i, S_i\} + W_i^{\text{volc}}) \quad (13)$$

Tᵢ and Rᵢ being a function of the atmospheric CO₂, CO₂ can be calculated for each time slice. Fᵢ vol is still fixed so that calculated CO₂ is just above the glacial threshold, and then held constant for the younger time slices.

Finally, given the spatial resolution of the model, and whatever the case, we assume that all solid or dissolved materials produced by erosion or weathering on one grid element are directly transferred to the ocean by rivers.

**Model calibration and validation.** Using the CRU-TS 2.1 database for temperature, the UNH/GRDC database for runoff and the SRTM database (30" initial resolution⁵³) for calculating the slope at a spatial resolution of 0.5° by 0.5°, we calculate the erosion fluxes and the silicate weathering fluxes on the modern land-sea configuration. kᵢ (equation (10)) is first fixed so that the model calculates a total discharge of sediments into the oceans⁵² equal to 20 Gt yr⁻¹. Once kᵢ is fixed, model outputs are compared to available data from a selected set of watersheds chosen for their contrasting climatic and tectonic conditions: the Congo river for warm and humid environments, the Mississippi river for temperate climates and the Mackenzie river for cold environments. In addition, we test the model for two watersheds strongly impacted by the presence of an orogen to check whether or not the model is reliable for active tectonic settings (the Amazon and Ganges-Brahmaputra). Erosion fluxes are in good agreement with field measurements, except for the Congo River, where the model overestimates the flux by a factor of ten (Supplementary Fig. 4). The model is particularly efficient for large catchments impacted by orogens, such as the Amazon and Ganges/Brahmaputra, suggesting its reliability for active orogenic areas.

The k Gr CO₂ constant (equation (10)) is then calculated so that the total CO₂ consumption by Ca- and Mg-silicate weathering, integrated over the whole continental surface, equals 2.15 × 10¹² molC yr⁻¹ (ref. 45). Then we check again the reliability of the model at the regional scale. The GEOCLIM-REG model is

able to reproduce the relative difference from one watershed to the other (Supplementary Fig. 5). Compared to the GEOCLIM simulation, the introduction of a weathering shielding factor considerably improves the data-model comparison, by limiting the weathering efficiency at low latitude.

Regarding the GEOCLIM-WEST simulation and using the best-fit parameters found in ref. 21, only a weak shielding effect is predicted in the equatorial region. The resulting weathering fluxes are overestimated for the Congo and Amazon watersheds. Acknowledging the spreading of the fitted data (fig. 2A in ref. 21), we found that the watershed data can be accurately reproduced by changing the numerical value of three of the parameters (Table 1 in ref. 21). Eₐ has been set to 14.6 kJ mol⁻¹, z to 20 t m⁻² and (σ + 1) to 1.13. All those new values are included within the 66% confidence range defined in the original publication. We used in the paper these latter sets of parameters when simulating the atmospheric CO₂ level across the late Palaeozoic.

**Code availability.** The key part of the code written for this study is made available at the following address: https://geoclimmodel.wordpress.com/download-area. A runnable version is available upon request to the first author.

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

References

35. Francois, L. M. & Godderis, Y. Isotopic constraints on the Cenozoic evolution of the carbon cycle. *Chem. Geol.* **145**, 177–212 (1998).

36. Kump, L. R. & Arthur, M. A. Interpreting carbon-isotope excursions: carbonates and organic matter. *Chem. Geol.* **161**, 181–198 (1999).

37. Berner, R. A., Lasaga, A. C. & Garrels, R. M. The carbonate–silicate geochemical cycle: an effects on atmospheric carbon dioxide over the past 100 million years. *Am. J. Sci.* **284**, 641–683 (1983).

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

39. Berner, R. A. & Caldeira, K. The need for mass balance and feedback in the geochemical carbon cycle. *Geology* **25**, 955 (1997).

40. Zeebe, R. & Caldeira, K. Close mass balance of long-term carbon fluxes from ice-core CO₂ and ocean chemistry records. *Nat. Geosci.* **1**, 312–315 (2008).

41. Jacob, R. L. *Low Frequency Variability in a Simulated Atmosphere Ocean System* PhD thesis, Univ. Wisconsin (1997).

42. Pierrehumbert, R. T. High levels of atmospheric carbon dioxide necessary for the termination of global glaciation. *Nature* **429**, 646–649 (2004).

43. Lefebvre, V., Donnadieu, Y., Sepulchre, P., Swingedouw, D. & Zhang, Z. S. Deciphering the role of southern gateways and carbon dioxide on the onset of the Antarctic Circumpolar Current. *Paleoceanography* **27**, PA002345 (2012).

44. Gough, D. O. Solar interior structure and luminosity variations. *Sol. Phys.* **74**, 21–34 (1981).

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

46. Dessert, C., Dupre, B., Gaillardet, J., Francois, L. M. & Allegre, C. J. Basalt weathering laws and the impact of basalt weathering on the global carbon cycle. *Chem. Geol.* **202**, 257–273 (2003).

47. Whipple, K. X. & Tucker, G. E. Dynamics of the stream-power river incision model: implications for height limits of mountain ranges, landscapes response timescales, and research needs. *J. Geophys. Res.* **104**, 17661–17674 (1999).

48. Heimsath, A. M., Hancock, G. R. & Fink, D. The 'humped' soil production function: eroding Arnhem Land, Australia. *Earth Surf. Process. Landf.* **34**, 1674–1684 (2009).

49. Gilbert, G. K. *Report on the Geology of the Henry Mountains (Utah)* (Department of the Interior, 1877).

50. Norton, K. P., Molnar, P. & Schlünegger, F. The role of climate-driven chemical weathering on soil production. *Geomorphology* **204**, 510–517 (2014).

51. Gabet, E. J. & Mudd, S. M. Bedrock erosion by root fracture and tree throw: a coupled biogeomorphic model to explore the humped soil production function and the persistence of hillslope soils. *J. Geophys. Res.* **115**, F04005 (2010).

52. Syvitski, J. P. M., Peckham, S. D., Hilberman, R. & Mulder, T. Predicting the terrestrial flux of sediment to the global ocean: a planetary perspective. *Sedim. Geol.* **162**, 5–24 (2003).

53. Farr, T. G. et al. The shuttle radar topography mission. *Rev. Geophys.* **45**, RG2004 (2007).