* Mes activités de recherche reposent sur l'utilisation des outils de la **modélisation numérique du climat mais aussi des cycles biogéochimiques (notamment celui du carbone) ainsi que de l'altération continentale.** L'intérêt de ces outils réside dans leur pouvoir de quantification des processus et rétroactions contrôlant le climat terrestre.

ResearcherID: A-5676-2011

## Thèses co-encadrées

* **Alexandre Pohl** (Soutenance effectuée le 16/11/2016) sujet : *Modélisation des changements climatiques et cryosphériques à l'Ordovicien*. Taux d'encadrement : 50%. Autre encadrant : Y. Donnadieu (DR CNRS).

* **Delphine Tardif** (thèse débutée le 01/10/2017) sujet : *Les moussons Asiatiques à la limité Eocène-Oligocène*. Taux d'encadrement : 33%. Autre encadrant : Y. Donnadieu (DR CNRS) et F. Fluteau (professeur IPGP).

* **Salomé Hennequin** (thèse débutée le 01/10/2017) sujet : *Etude des couplages entre l'évolution géologique des cycles biogéochimiques, du climat et de la biodiversité océanique*, Taux d'encadrement : 33%. Autre encadrant : G. Goddéris (DR CNRS) et F. Fluteau (professeur IPGP).

* **Mickael Mussard** (thèse abandonnée en Septembre 2014 suite à des soucis familiaux). sujet : *Quantification des perturbations environnementales liées aux traps du Deccan*. Début de thèse 01/10/2011 -Durée de la thèse 3ans. Taux d'encadrement : 50%. Autre encadrant : F. Fluteau (professeur IPGP).

* **Yoram Teitler** (Soutenance effectuée le 21/12/2012) sujet : *Climat et composition de l'atmosphère au Précambrien : des paléosols aux simulations numériques*. Durée de la thèse 4ans et 3mois. Taux d'encadrement : 33%. Autres encadrants : P. Philippot (Professeur -CNRS-IPGP), F. Fluteau (professeur IPGP).

## Stages Masters

**Clément Dufour** (M1) sujet : « **évolution de l'albédo planétaire** ». 2011 (co-encadré avec Stéphane Jacquemoud, Université Dénis Diderot - Paris 7), **Nathan Jeiger** (M1) sujet : « **les perturbations climatiques liées aux changements de circulations océaniques au Cénozoïque** » (stages co-dirigé avec Y. Donnadieu, CNRS, LSCE centre de Saclay), **Pauline Galmiche** (M1) sujet : « **Analyse de la stabilité de la calotte Ouest Antarctique, quel rôle pour la surrection andine et le passage de Panama** » (stages co-dirigé avec Y. Donnadieu et P. Sepulchre, CNRS, LSCE centre de Saclay), **Alexandre Pohl** (M2) sujet : « **Modélisation des changements climatiques et cryosphériques à l'Ordovicien.** 2nd semestre de l'année 2012-2013. Taux encadrement 33% (co-encadré par Yannick Donnadieu (CNRS), Christophe Dumas (CEA) du LSCE centre de Saclay)**

## Stages Licences

**Léna Gurriaran** (L2) sujet : « le climat du Crétacé, l'importance de l'altération des silicates de la croûte océanique » Taux encadrement 100% (stage non obligatoire)
# RESPONSABILITES

* **Responsabilités au sein de l'université Paris 7 Denis Diderot - IPGP)**

* Depuis 2016 : Directeur de la Licence Sciences de la Vie et de Terre (univ. Paris 7)

* 2012 – 2017 : Membre élu de la commission Recherche de l'Université Paris 7

* 2012 – 2017 : Membre élu de la commission des moyens de l'Université Paris 7

* 2012 – 2017 : Membre du conseil de l'UFR (université Paris 7-IPGP)

* 2014-2015 : Directeur adjoint de l'UFR (université Paris 7-IPGP)

* **Responsabilités nationales**

* depuis 2019 : évaluateur HCERES pour les formations

* 2014-2017 : jury au concours de l'Agrégation externe de SVT

* depuis 2014 : membre du du GDR Climat Environnement Vie (CEV) : les apports du pré-Quaternaire et co-responsable de la thématique - Volcanisme et Perturbations extrêmes du système Terre

* depuis 2012 : expert au CT1 (Comité thématique sur le climat) au comité GENCI (grand équipement national de calcul intensif)

* **Responsabilités Scientifiques**

* 2018-2019 Responsable du projet INSU (comité SYSTER) « Paléomousson : de l'enregistrement continental au forçage orbital en Asie à l'Eocène »

* 2017-2018 Participant au projet INSU (comité SYSTER) : « Paléomousson » piloté par F. Fluteau (IPGP)

* 2017-2018 Participant au projet INSU (comité INTERVIE) : « Couplages des cycles biogéochimiques, du climat et de la biodiversité à l'échelle des temps géologiques » piloté par Y. Goddéris (GET)

* 2016-2017 Participant au projet INSU (comité SYSTER) : « Conséquences du grand bombardement tardif sur la Terre archéenne » piloté par F. Fluteau (IPGP)

* 2015-2016 Responsable (P.I.) du projet INSU (comité SYSTER) : les glaciations du Néoprotérozoïque : un test pour comprendre la réponse du système Terre aux événements extrêmes à l'échelle de temps géologiques.

* 2015-2016 Participant au projet INSU (comité SYSTER) : « Evolution de la production carbonatée pélagique à la limite K/Pg » piloté par B. Sucheras (CEREGE)

* 2013-2015 Participant au Svalbard project (PI. I Faichild et D. Benn) « étude de la Terre boule de neige » en tant que responsable du pôle modélisation
* 2013-2014 - Responsable (P.I.) du projet INSU (comité SYSTER) : Modélisation des effets des grandes provinces basaltiques sur l'environnement et la productivité marine (2013-2014)

* 2009-2012 Participant au projet Emergeance « quantifier et comprendre les variations des températures océaniques passées – projet modèle – données », responsable de la partie modélisation climat et cycles biogéochimiques (P.I. Bonifaci Magali, arrêt du projet : fin 2012)

* 2009-2012 - Participation à l'ANR eLife en tant que responsable de la partie modélisation des climats anciens (P.I. de l'ANR eLife : Philippot Pascal, IPGP, arrêt de l'ANR : fin 2012)

# ENSEIGNEMENT

J'enseigne essentiellement dans 3 grands domaines : **le climat, la biologie, et la géologie de terrain.**

i) les filières de formations aux futurs enseignants (40 %) - responsable d'une UE

ii) stage de terrain – modélisation 3D des structures géologiques (40%)

iii) en licence de géologie (20 %) - co-responsable de 2 UEs.

## Présentation synthétique des enseignements

<table><thead><tr><th rowspan="2">UE</th><th colspan="7">Heures Equivalents TD (CM=1,5TD-TP)</th></tr><tr><th>Filière</th><th>Niveau</th><th>Effectif</th><th>Formation</th><th>CM</th><th>TD</th><th>TP-Terrain</th><th>Heures Totales</th></tr></thead><tbody><tr><td>Stage de terrain (Languedoc)</td><td>LVT</td><td>L2</td><td>40</td><td>initiale</td><td>0</td><td>0</td><td>30</td><td>30</td></tr><tr><td>Module de cartographie</td><td>LVT</td><td>L3</td><td>40<br>(2 groupes)</td><td>initiale</td><td>0</td><td>0</td><td>24</td><td>24</td></tr><tr><td>Préparation aux écrits</td><td>MEEF</td><td>M1</td><td>40</td><td>initiale</td><td>7</td><td>15</td><td>0</td><td>25,5</td></tr><tr><td>Préparation aux oraux</td><td>MEEF</td><td>M1</td><td>20<br>(2 groupes)</td><td>initiale</td><td>0</td><td>12</td><td>0</td><td>12</td></tr><tr><td>Exploitation géologique du sous-sol de Paris</td><td>MEEF</td><td>M1-M2</td><td>40<br>(2 groupes)</td><td>initiale</td><td>0</td><td>0</td><td>6</td><td>6</td></tr><tr><td>Stage de cartographie 3D (12 jours)</td><td>Masters<br>géophysique</td><td>M1-M2</td><td>16</td><td>Initiale-<br>continue</td><td>0</td><td>0</td><td>72</td><td>72</td></tr><tr><td>Biologie, Biochimie pour les géosciences</td><td>STEP</td><td>L1</td><td>72</td><td>Initiale</td><td>12</td><td>0</td><td>0</td><td>18</td></tr><tr><td>L'Homme et la Planète</td><td>STEP</td><td>L2</td><td>72</td><td>Initiale</td><td>10</td><td>0</td><td>0</td><td>15</td></tr></tbody></table>

<table><thead><tr><th>Année</th><th>Heures équivalents TD</th></tr></thead><tbody><tr><td>2015-2016</td><td>202</td></tr><tr><td>2016-2017</td><td>261,5</td></tr><tr><td>2017-2018</td><td>199,5</td></tr><tr><td>2018-2019</td><td>216,5</td></tr></tbody></table>

<footnote>* STEP Sciences de la Terre, de l'Environnement et des Planètes
Responsable de la Licence VT (décharge 32h)</footnote>
• Pédagogie

Ayant été moi-même enseignant dans le secondaire entre 2002-2003, je me suis impliqué, depuis
mon recrutement à l'université Paris VII, dans les filières de formation aux futurs enseignants de
SVT (ex prépa-CAPES, actuelles ESPE). J'ai eu l'opportunité d'être jury au concours de
l'agrégation externe de SVT.

De ces expériences, j'ai décidé de m'investir en Licence (~ 60 % de mon temps – hors stage de terrain). Avec la réforme LMD, la licence est devenue le niveau crucial d'acquisition, à une période où les étudiants perçoivent encore mal l'intérêt des principes étudiés - spécialement si ces derniers sont présentés frontalement. Pour cette raison j'ai choisi de développer des approches pédagogiques :

* En supprimant les présentations type conférence et en réduisant les supports visuels à un nombre très réduit (3-4/2h de cours). En effet les cours magistraux sous forme de présentation type conférence (30-50 slides/heure) donnent un nombre d'informations visuelles bien trop grand pour un esprit qui n'est pas encore formé à sélectionner l'information. L'étudiant éprouve des difficultés à hiérarchiser l'importance des notions évoquées.

* Pour remplacer ce mode de présentation, j'ai choisi de repartir de la craie et du tableau noir (ou d'une mise en situation) afin de faire émerger une grande question, à partir de laquelle les notions d'importance sont construite. Cette approche, classique du secondaire, n'est possible que si le nombre d'étudiants reste modéré et devient délicate à appliquer dans un amphithéâtre.

* Lorsque cela est opportun, j'ai aussi opté pour l'insertion de l'étudiant dans la société. Cette approche repose sur la discussion du bien fondé scientifique de questions abordés dans l'actualité. La finalité de cette approche est de montrer comment est construit la savoir scientifique afin de permettre à l'étudiant de pouvoir discuter d'un point de vue scientifique l'actualité (ex : le réchauffement climatique et les derniers communiqués de presse du GIEC) et donc de s'affranchir des fakenews. Tous les cours ne sont pas propices à cette approche, cette approche est tout particulièrement développée dans l'UE « L'Homme et la Planète » où il est questions des problèmes environnementaux.

CONFERENCES (dernières années)

2017. Séminaire-invité à l'université de Lausanne (Suisse) en 2017 : « l'effondrement des producteurs primaires à la limite Crétacé-Paléogène »

2019. Séminaire-invité iSTEP (UMPC Paris) en 2019 « L'Apport de la modélisation biogéochimique à la compréhension de la crise K-Pg »

2019. Séminaire au GDR Climat Environnement Vie (CEV) – Lyon : « L'Apport de la modélisation biogéochimique à la compréhension de la crise K-Pg »

2019. Séminaire d'équipe au LSCE – Gif Sur Yvette - « L'Apport de la modélisation biogéochimique à la compréhension de la crise K-Pg »

Communications

2005-2019: Nombreuses communications orales, posters en conférences internationales (Goldschmidt, EGU, GSA, AGU) et nationales (RST 2104, 2016)
Publications de rang A

1. 2007 Le Hir, G., Ramstein, G., Donnadieu, Y., Pierrehumbert, R. Investigating plausible mechanisms to trigger a deglaciation from a hard snowball Earth. *Comptes Rendus de l'Aca-démie des Sciences, snowball Earth special issue*, tome 339, n°3 et 4, p274-287,

2. 2008 Le Hir, G., Goddéris, Y., Ramstein, G., and Donnadieu, Y., Scenario for the evolution of atmospheric pCO2 during a snowball Earth. *Geology*, vol n36, n°1, p47-50,

3. 2008 Le Hir, G., Goddéris, Y., Donnadieu, Y., and Ramstein, G. A geochemical modeling study of the evolution of the chemical composition of seawater linked to a “snowball” gla-ciation. *Biogeosciences*, vol n°5, issue n°1, p253-267,

4. 2009 Le Hir, G., Goddéris, Y., Donnadieu, Y., Pierrehumbert, R, Halverson, G., Macouin M, Nedelec A., and Ramstein, G, The snowball Earth aftermath: exploring the limits of conti-nental weathering processes, *Earth and Planetary Science Letter*, vol 277, issue 3-4, p453-463,

5. 2010_Le Hir, G., Donnadieu, Y., Krinner, G., and Ramstein, G. Toward Toward the snowball earth deglaciation... *Climate Dynamics*, vol 35, p285-297,

6. 2011 Le Hir, G., Donnadieu, Y., Goddéris, Y., Meyer-Berthaud, B., Ramstein, G., Blakey, R., The climate change caused by the land plant invasion in the Devonian, Earth and Planetary Science Letter, vol 310, p203-212

7. 2011 Nardin, E., Goddéris, Y., Donnadieu, Y., Le Hir, G., Blakey, R., Pucéat, E., Aretz, M., Modeling the early Paleozoic long-termed climatic trend, Geological Society of America Bulletin. Vol. 123, p: 1181-1192, doi: 10.1130/B30364.1

8. 2012 Abbot, D., Voigt, A., Branson, M; Pierrehumbert R., Pollard, D., Le Hir, G, Koll, D. Clouds and Snowball Earth deglaciation, Geophysical Research Letters. Vol. 39. DOI: 10.1029/2012GL052861

9. 2012 Goddéris, Y., Donnadieu, Y. Lefebvre, V., Le Hir, G ; Nardin, E 2012, Tectonic control of continental weathering, atmospheric CO2, and climate over Phanerozoic times Comptes rendus Geosciences. Vol: 344 Issue: 11-12 Pages: 52-662 DOI: 10.1016/j.crte.2012.08.009

10. 2013 Abbot, D., Voigt, A, Li, D. Le Hir, G ; Pierrehumbert, RT ; Branson, M ; Pollard, D ; Koll, D. Robust elements of Snowball Earth atmospheric circulation and oases for life. Journal of Geophysical Research-Atmospheres. Vol: 118 Issue: 12 Pages: 6017-6027 DOI: 10.1002/jgrd.50540

11. 2013 Goddéris Y.; Donnadieu Y.; Le Hir G.; Lefebvre V.; Nardin E. The role of paleogeography in the Phanerozoic history of atmospheric CO2 and climate. Earth Science Reviews. 128,122-138, doi. 10.1016/j.earscirev.2013.11.004

12. 2014 *Mussard M, Le Hir G., Fluteau F., Lefebvre F., Goddéris Y. Modeling the carbon-sulfate interplays in climate changes related to the emplacement of Continental Flood
Basalts., Geological Society of America vol; 505, p339-352, 2014. doi: 10.1130/2014.2505(17).

13. 2014 *Teitler, Y; Le Hir, G; Fluteau, F; Philippot, P; Donnadieu, Y; Investigating the Paleoproterozoic glaciations with 3-D climate modeling, Earth and Planetary Science Letters, vol. 395, p:71-80, doi:10.1016/j.epsl.2014.03.044

14. 2014 Le Hir, G., *Teitler, Y., Fluteau, F., Donnadieu, Y., and Philippot, P., The faint young Sun problem revisited with a 3-D climate-carbon model-Part 1,Climate of the Past, vol. 10, p:697-713, doi:10.5194/cp-10-697-2014

15. 2014 *Pohl, A; Donnadieu, Y; Le Hir, G; Buoncristiani, J-F; Vennin, E,; Effect of the Ordovician paleogeography on the (in) stability of the climate, Climate of the Past, 10, 2053-2066, doi:10.5194/cp-10-2053-2014

16. 2014 Godderis, Y.; Donnadieu, Y.; Le Hir, G. Lefebvre V., Nardin, E,; The role of palaeogeography in the Phanerozoic history of atmospheric CO2 and climate, Earth-Science Reviews, 128, 122-138, DOI: 10.1016/j.earscirev.2013.11.004

17. 2015 *Teitler, Y., Philippot, P., Gérard, M., Le Hir G., Fluteau F., Ader M., Ubiquitous occurrence of basaltic-derived paleosols in the Late Archean Fortescue Group, Western Australia. PrecambrianResearch. doi:10.1016/j.precamres.2015.05.014

18. 2015 Benn D. Le Hir G., Bao H., Donnadieu Y, Dumas C., Fleming E.J., Hambrey M., McMillan E., Petronis M., Ramstein G., Stevenson C., Wynn P., Fairchild I. Orbitally forced ice sheet fluctuations during the Marinoan Snowball Earth glaciations. Nature Geoscience, vol.8, p704-707,, doi:10.1038/ngeo2502

19. 2015 Hoareau G., Bomou B., van Hinsbergen D., Carry N., Marquer D., Donnadieu Y., Le Hir G., Vrielynck G., Walter-Simonnet A.V. Did high Neo-Tethys subduction rates contribute to early Cenozoic warming? Climate of the Past, vol. 11, pp1751-1767, doi:10.5194/cp-11-1751

20. 2016 Pohl A., Donnadieu Y., Le Hir G.*, Ladant J.B., Dumas C., Alvarez-Solas J., Vandenbroucke T. Glacial onset predated Late Ordovician climate cooling. Paleoceanography, vol. 31, pp800-821. doi:10.1002/2016PA002928

21. 2016 Godderis, Y.; Le Hir, G.; Macouin, M.; Donnadieu Y., Hubert-Théou Lucie, Dera G., Aretz M., Fluteau F., Lie Z.X., Halverson G. Paleogeographic forcing of the strontium isotopic cycle in the Neoproterozoic Gondwana Research Volume: 42 Pages: 151-162 DOI: 10.1016/j.gr.2016.09.013

22. 2017 Charnay B., Le Hir G., Fluteau F., Forget F., Catling D.C., A warm or a cold early Earth? New insights from a 3-D climate-carbon model, Volume 474, Pages 97-109. doi.org/10.1016/j.epsl.2017.06.029

23. 2017 Van der Meer, D. G.; van Saparoea, A; van Hinsbergen, D. J. J.; van de Wega R.M.B, Godderis Y., LeHir G., Donnadieu Y. Reconstructing first-order changes in sea level during the Phanerozoic and Neoproterozoic using strontium isotopes. Gondwana Research Volume: 44 Pages: 22-34. DOI: 10.1016/j.gr.2016.11.002
24. 2017 Pohl, A.*, Donnadieu, Y., Le Hir, G. and Ferreira, D.. The climatic significance of Late Ordovician-early Silurian black shales. *Paleoceanography*, 32(4): 397-423. DOI10.1002/2016PA003064

25. 2017 Hoffman P., Abbot D., Ashkenazy Y., Benn D., Brocks J.J., Cohen P., Cox G., Creveling J., Donnadieu Y., Erwin D., Fairchild I., Ferreira D., Goodman J., Halverson G., Jansen M., Le Hir G., Love G., Macdonald F., Maloof A., Partin C., Ramstein G., Rose B., Rose C., Sadler P., Tziperman E., Voigt A., Warren S. Snowball Earth climate dynamics and Cryogenian geology-geobiology, Science Advances, vol 3(11): DOI10.1126/sciadv.1600983

26. 2018 Pohl, A.*, Harper D., Donnadieu, Y., Le Hir G., Nardin E., Servais T. Possible patterns of marine primary productivity during the Great Ordovician Biodiversification Event, Lethaia, DOI 10.1111/let.12247

27. 2019 Michaud-Dubuy A., Carazzo G., Tait, S., Le Hir G, Fluteau F., Kaminiski E. Impact of wind direction variability on hazard assessment in Martinique (Lesser Antilles): the example of the 14 ka cal BP Bellefontaine Plinian eruption of Mount Pelée volcano Journal of Volcanology and Geothermal Research vol 381, Pages 193-208, doi.org/10.1013/j.volgeores.2019.06.004

28. Le Hir G., Fluteau, F., Sucheras, B., Goddéris, Y., Amplifying Factors Leading to the Collapse of Primary Producers during the Chicxulub Umpact and Deccan Traps, accepted Geological Society of America

Autres Publications

Chapitres de livres

Goddéris, Y., Le Hir, G., Donnadieu, Y., 2011, Modelling the Snowball Earth, "Geological Society, London, Memoirs",36,1,151-161, Geological Society of London, editor Armaud E, Shields G., Halverson G.

Le Hir, G., Paleomagnetism and climate, 2011, CNFGG

Donnadieu, Y Goddéris, Y., Le Hir, G., 2012, Neoproterozoic atmospheres and glaciation; *Treatise of Geochemistry vol. 13*, Editor Donald E. Canfield, James Farquhar, James F. Kasting

Le Hir, G. et Fluteau, F., 2012, la Climatologie, Encyclopedia Universalis

Goddéris, Y., Ramstein G., Le Hir, G., 2013, Le climat du Précambrien, in *Paléoclimatologie, CNRS editions*

Godderis, Y; Donnadieu, Y; Ramstein, G., Pucéat, E., Le Hir G., 2013, Le climat de la Terre à l'échelle destemps géologiques in Le Climat à découvert, CNRS éditions

Godderis, Y; Donnadieu, Y; Ramstein, G., Pucéat, E., Le Hir G., 2013, Le cycle du carbone à l'échelle des temps, CNRS éditions
Le Hir, G., La vie sur la Terre glacée, *La Recherche*, 2005

Le Hir, G., Goddéris, Y., Planète Blanche, *Sciences et Vie, Hors-série*, 2006

Le Hir, G. Goddéris, Y., La Terre boule de neige, *Les Dossiers de la Recherche*, 2007

Le Hir G. « la Terre boule de neige » bulletin de l'APBG (Association des professeurs de Biologie et Géologie) 2015

P. SanJoffre et G. Le Hir « La Terre boule de Neige, «Cataclysme climatique et tremplin biologique » Pour la Science, 2018

Pohl, A, Donnadieu, Y, Le Hir G., la modélisation peut elle aider à comprendre le climat d'il y a 450 Ma ? La Météorologie, n° 105, 2019 (article invité suite à la remise du prix Prud'homme à A. Pohl pour sa thèse)
RECUEIL DES TRAVAUX

Une sélection de 10 papiers est présentée dans la suite de ce manuscrit et classée par ordre
chronologique de publication :

1. 2008 Le Hir, G., Goddéris, Y., Ramstein, G., and Donnadieu, Y., Scenario for the evolution of atmospheric pCO2 during a snowball Earth. *Geology*, vol n36, n°1, p47-50.

2. 2009 Le Hir, G., Goddéris, Y., Donnadieu, Y., Pierrehumbert, R, Halverson, G., Macouin M, Nedelec A., and Ramstein, G, The snowball Earth aftermath: exploring the limits of continental weathering processes, *Earth and Planetary Science Letter*, vol 277, issue 3-4, p453-463,

3. 2011 Le Hir, G., Donnadieu, Y., Goddéris, Y., Meyer-Berthaud, B., Ramstein, G., Blakey, R., The climate change caused by the land plant invasion in the Devonian, Earth and Planetary Science Letter, vol 310, p203-212

4 2014 Le Hir, G., *Teitler, Y., Fluteau, F., Donnadieu, Y., and Philippot, P., The faint young Sun problem revisited with a 3-D climate-carbon model-Part 1,Climate of the Past, vol 10, p:697-713, doi:10.5194/cp-10-697-2014

5. 2014 *Pohl, A; Donnadieu, Y; Le Hir, G; Buoncristiani, J-F; Vennin, E, Effect of the Ordovician paleogeography on the (in) stability of the climate, Climate of the Past, 10, 2053-2066, doi:10.5194/cp-10-2053-2014

6. 2014 Godderis, Y.; Donnadieu, Y.; Le Hir, G. Lefebvre V., Nardin, E, The role of palaeogeography in the Phanerozoic history of atmospheric CO2 and climate, Earth- Science Reviews, 128, 122-138, DOI: 10.1016/j.earscirev.2013.11.004

7. 2015 Benn D. Le Hir G., Bao H., Donnadieu Y, Dumas C., Fleming E.J., Hambrey M., McMillan E., Petronis M., Ramstein G., Stevenson C., Wynn P., Fairchild I. Orbitally forced ice sheet fluctuations during the Marinoan Snowball Earth glaciations. Nature Geoscience, vol.8, p704-707, doi:10.1038/ngeo2502

8. 2016 Pohl A*, Donnadieu Y., Le Hir G. Ladant J.B., Dumas C., Alvarez-Solas J., Vandenbroucke T. Glacial onset predated Late Ordovician climate cooling. Paleoceanography, vol. 31, pp800-821, doi:10.1002/2016PA002928

9. 2017 Charnay B., Le Hir G., Fluteau F., Forget F., Catling D.C., A warm or a cold early Earth? New insights from a 3-D climate-carbon model, Volume 474, Pages 97-109. doi.org/10.1016/j.epsl.2017.06.029

10 Le Hir G., Fluteau, F., Sucheras, B., Goddéris, Y., Amplifying Factors Leading to the Collapse of Primary Producers during the Chicxulub Impact and Deccan Traps, sous presse Geological Society of America
# Scenario for the evolution of atmospheric pCO₂ during a snowball Earth

Guillaume Le Hir
Gilles Ramstein
Yannick Donnadieu
UMR CEA/CNRS/UVSQ, 1572 Laboratoire des Sciences du Climat et de l'Environnement, CE Saclay,
Orme des Merisiers, Bât. 701, 91191 Gif sur Yvette Cedex, France

Yves Goddéris UMR CNRS, 5563 Laboratoire des Mécanismes et Transferts en Géologie, Observatoire Midi-Pyrénées,
14, avenue Edouard Belin, 31400 Toulouse, France

## ABSTRACT

The snowball Earth theory, initially proposed by J.L. Kirschvink to explain the Neoproterozoic glacial episodes, suggests that the Earth was globally ice covered at 720 Ma (Sturtian episode) and 640 Ma (Marinoan episode). The reduction of the water cycle and the growth of large ice sheets led to a collapse of CO₂ consumption through continental weathering and biological carbon pumping. As a consequence, atmospheric CO₂ built up linearly to levels allowing escape from a snowball Earth. In this contribution, we question this assumed linear accumulation of CO₂ into the atmosphere. Using a numerical model of the carbon-alkalinity cycles, we suggest that during global glaciations, even a limited area of open waters (10³ km²) allows an efficient atmospheric CO₂ diffusion into the ocean. This exchange implies that the CO₂ consumption through the low-temperature alteration of the oceanic crust persists throughout the glaciation. Furthermore, our model shows that rising CO₂ during the glaciation increases the efficiency of this sink through the seawater acidification. As a result, the atmospheric CO₂ evolution is asymptotic, limiting the growth rate of the atmospheric carbon reservoir. Even after the maximum estimated duration of the glaciation (30 m.y.), the atmospheric CO₂ is far from reaching the minimum deglaciation threshold (0.29 bar). Accounting for this previously neglected carbon sink, processes that decrease the CO₂ deglaciation threshold must be further explored.

**Keywords:** carbon cycle, seafloor weathering, snowball Earth, Neoproterozoic.

## INTRODUCTION

Although processes and conditions required to trigger a snowball Earth (Kirschvink, 1992) have been intensively explored (Donnadieu et al., 2002, 2004a, 2004b; Schrag et al., 2002; Warren et al., 2002), studies related to the deglaciation of a snowball Earth remain scarce (Caldeira and Kasting, 1992; Higgins and Schrag, 2003; Pierrehumbert, 2005). The amount of atmospheric CO₂ needed to trigger the deglaciation from a globally ice-covered Earth was first evaluated at 0.12 bar (Caldeira and Kasting, 1992) and then modified at 0.29 bar by Pierrehumbert (2004) to account for the decreased solar luminosity in Neoproterozoic time. Using a more complex and comprehensive model (an atmospheric general circulation model, i.e., AGCM), Pierrehumbert (2005) demonstrated that the threshold value of 0.29 bar calculated with an energy balance model (EBM) was underestimated, and suggested that this value is at the lower end of the CO₂ levels required to deglaciate the snowball Earth.

In this contribution we focus on the various processes that may affect the way CO₂ would accumulate in a snowball atmosphere. It has been assumed that the drastic reduction of the continental weathering and the primary productivity occurring in a snowball Earth would have left the CO₂ outgassed by volcanoes free

to accumulate in the atmosphere. Hoffman et al. (1998) calculated that atmospheric pCO₂ would reach 0.12 bar at the present outgassing rate after ~4 m.y. if no air-sea gas exchange took place, but that partial gas exchange through cracks in sea ice would increase this estimate. In contrast, Higgins and Schrag (2003) assumed that the ocean and atmosphere were in equilibrium with respect to CO₂ due to the air-sea gas exchange during the entire glacial event. However, in this later case an important additional factor able to modify the seawater chemistry, i.e., the low-temperature alteration of the oceanic crust, has to be considered if air-sea gas exchange took place. Because a large amount of the atmospheric CO₂ could be dissolved into the ocean, the ocean acidification promotes the removal of carbon through basaltic crust dissolution. This sink may have moved the simple linear atmospheric CO₂ accumulation toward an asymptotic evolution, raising the question of whether deglaciation melting can be reached.

## CARBON CYCLE AND SEAFLOOR WEATHERING DURING SNOWBALL EVENTS

To investigate and quantify the carbon cycle evolution in a snowball environment, we designed a simple numerical model of the

carbon-alkalinity cycles that includes the calculation of the seafloor weathering response to the acidification of the ocean (see Appendix). The model includes a box for the atmosphere, while the ocean is partitioned between a surface and a deep reservoir connected to a reservoir of water percolating through the oceanic crust, with a residence time of $2 \times 10^4$ yr (Elderfield and Schultz, 1996). The CO₂ exchange flux, F, at the ocean-atmosphere interface is calculated as a function of the gradient between the calculated atmospheric ($pCO₂^{atm}$) and dissolved oceanic $pCO₂$ ($pCO₂^{ocean}$):

$$F = \alpha(pCO_2^{atm} - pCO_2^{ocean}) \text{area}_{\text{open}}, \quad (1)$$

where α is the diffusion coefficient fixed at 16 mol/(yr m² PAL [Present Atmospheric Level]) (Etcheto et al., 1991), a value typical of high-latitude cold seawater, and areaopen is the area of open water, a value considered as a model parameter that can be tested. $pCO₂^{ocean}$ is estimated using Henry's law, fixing surface ocean temperature at 275 K during the snowball Earth. Carbonate speciation is calculated for each water reservoir at each time step. We assume an oceanic vertical mixing 100 times lower than the present-day thermohaline circulation (0.02 Sv) between the surface and deep ocean to account for the maximal expected overturning reduction due to the ice-covered ocean, because it seems physically improbable to completely shut down the general circulation (Zhang et al., 2001). Low-temperature seafloor alteration is calculated using laboratory kinetic laws (see Appendix). The dissolution of each mineral constitutive of the fresh basaltic crust (assumed to be a tholeiitic basalt) is modeled as a function of the temperature, T, of the water in contact with the basalt (fixed at 40 °C) (Brady and Gislason, 1997), and the pH of the percolating waters. The conclusions of this study partly depend on the model parameters, including the vertical mixing of the ocean, the present-day amplitude of the seafloor weathering, the water flux percolating the oceanic crust as well as its temperature, the seafloor spreading rate, the presence or absence of seafloor carbonates, and the area of open water. In the following we concentrate on the two last parameters, which are by far the most important. More details and
diagrams about sensitivity tests, regarding other important parameters such as the vertical mixing of the ocean and the efficiency of the seafloor weathering, can be found in GSA Data Repository Figures DR1 and DR2.<fnref n="1" />

## Efficient Exchange Between the Ocean and Atmosphere

We first tested the minimal surface of open waters allowing an efficient diffusion of the atmospheric CO₂ into the ocean. Patches of open water were proposed in the original snowball Earth theory as refugia for living organisms (Kirschvink, 1992). Their existence has been confirmed by the study of Neoproterozoic formations (Hoffman and Maloof, 1999). To test the efficiency of the CO₂ diffusion from the atmosphere into the ocean, we define the diffusion efficiency E:

$$E = 1 - \frac{pCO_2(\text{area}_{\text{open}})}{pCO_2(\text{hard snowball})}, \quad (2)$$

where $pCO_2(\text{area}_{\text{open}})$ is the CO₂ level reached after a simulation of 30 m.y., the maximum estimated duration of the snowball event (Hoffman et al., 1998), for a given free water area ($\text{area}_{\text{open}}$), and $pCO_2$ (hard snowball) is the same value assuming an ocean fully covered with ice. This E parameter is then normalized to its present-day value ($\text{area}_{\text{open}} = 363 \times 10^6 \text{ km}^2$). In our set of model simulations, wherein the ice-free surface is varied from the present-day ocean surface ($363 \times 10^6 \text{ km}^2$) to 0 km², we calculate that above 3000 km², atmospheric CO₂ has enough time at the geologic time scale (here several millions of years) to dissolve within the ocean, and the area of open water is not a limiting factor for diffusion (Fig. 1). Below 3000 km², atmospheric CO₂ mixing into the ocean is limited. Today the surface of emerged active hydrothermal areas largely exceeds the $10^3 \text{ km}^2$ scale, being equal to $1.5 \times 10^6 \text{ km}^2$ (Dessert et al., 2003). During global glaciations, emerged geothermal environments may be associated with small open water areas, maintained by heat production and lava flows. However, the existence of thick sea ice (Goodman, 2006) may reduce the contact between the patches of free waters and the whole ocean, limiting the CO₂ dissolution. For that reason, we cannot firmly state that the atmospheric CO₂ diffusion was fully efficient, but only 3000 km² of open waters in contact with the whole ocean are enough to achieve massive CO₂ dissolution.

<figure><img src="image_41.png" /><figcaption>Figure 1. Normalized CO₂ diffusion efficiency as function of the ice-free oceanic surface.</figcaption></figure>

## Implications

If CO₂ is allowed to dissolve within the ocean during the building up of extremely high atmospheric carbon levels, seawater pH should drastically decrease and promote both deep-sea carbonates and oceanic basaltic crust dissolution.

## Seafloor Weathering and Its Implication on the $pCO_2$

In the absence of deep seafloor sedimentary carbonates, the seawater acidification enhances the consumption of dissolved carbon through low-temperature alteration of the oceanic basaltic crust. Under the corrosive action of seawater percolating the crust at almost kilometric depth, basalts dissolve and release alkalinity, which is then removed through carbonate precipitation in veins from highly saturated pore waters (Alt and Teagle, 1999; Brady and Gislason, 1997; Spivack and Staudigel, 1994). The CO₂ consumption through seafloor basaltic weathering rises until the ocean is sufficiently acidified, so that consumption of carbon through seafloor weathering balances solid Earth outgassing.

The question is whether this theoretical steady-state CO₂ is higher than the required deglaciation threshold. In our scenario accounting for seafloor weathering, the atmospheric CO₂ build up does not rise linearly as in the classical snowball scenario (Fig. 2), but reaches an asymptotic value. Assuming a present-day total degassing rate of $7.15 \times 10^{12} \text{ mol/yr}$ (a value fixed to balance the present continental silicate and seafloor basalt CO₂ consumption rates) (Brady and Gislason, 1997; Gaillardet et al., 1999), we found that the atmospheric partial pressure does not exceed 0.17 bar within the estimated snowball duration of 12 m.y. (Bodiselitsch et al., 2005), because the seawater pH goes rapidly down to 6, increasing the seafloor weathering rate by a factor slightly above 4. Even after 30 m.y., the maximal duration for the snowball event (Hoffman et al., 1998), and with a volcanic outgassing rate enhanced by 20% (which is the maximum increase recorded over the past

<figure><img src="image_42.png" /><figcaption>Figure 2. CO₂ evolution as function of the snowball Earth duration; 30 m.y. is maximum estimate. One run was performed with 10 km² of ice-free surface, all others with ice-free surface of 3000 km². Snowball earth deglaciation threshold (blue line) is 0.29 bar, estimated through an energy balance model simulation for a Neoproterozoic insolation. With a modern outgassing rate, the difference between the run with or without the seafloor weathering (black thick line and black dashed line, respectively) shows the impact of the CO₂ consumption by the low-temperature weathering of the basaltic seafloor. Higher volcanic degassing rate (red line) ($1.2 \times$ present day; see text) would help at rising more rapidly the atmospheric CO₂ level but not sufficiently to overcome the deglaciation threshold. Because total degassing depends on the seafloor-spreading rate (Gaffin, 1987), the CO₂ consumption increases proportionally the surface of fresh seafloor basaltic rocks exposed to seawater. Below 10 km² of ice-free surface, the CO₂ diffusion becomes so weak that the atmosphere and ocean remain virtually isolated, the CO₂ accumulates in the atmosphere, and the seafloor weathering effect becomes negligible.</figcaption></figure>

160 m.y.; Cogne and Humler, 2004), the CO₂ threshold of 0.29 bar required to escape the snowball glaciation cannot be reached (Fig. 2).

## Buffering of Seawater pH through Carbonate Dissolution or Basaltic Weathering?

Seafloor sedimentary carbonate dissolution may have buffered seawater pH during the snowball glaciation, thus reducing the efficiency of seafloor weathering. However, maintain-

<footnote n="1">¹GSA Data Repository item 2008013, Figures DR1–DR3, supplementary information on the model and additional tests, is available online at www.geosociety.org/pubs/ft2008.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.</footnote>
ing carbonate saturation (Higgins and Schrag, 2003) under 0.12 bar of CO₂ would require a calcium concentration in seawater twice its present-day level, corresponding to the dissolution of 2.5 × 10¹⁹ mol of carbonates if preglacial calcium concentration was close to its present-day value. If carbonate were of shelfal origin (pelagic carbonate producers being virtually absent), and considering the maximum past carbonate shelf extent of 15 × 10⁶ km² (Walker et al., 2002), this would correspond to the dissolution of 66 m of carbonate sediments. This value would be even greater since the 0.12 bar threshold is largely underestimated (Pierrehumbert, 2004). Because sea level fell by more than 400 m during the glacial event (Hoffman et al., 2007), these shelf carbonates might not be available for dissolution.

However, carbonate slumps along continental slopes (Swart, 1992) or preglacial direct abiotic precipitation might have partly overlaid the deep seafloor with carbonates. We tested this hypothesis by running the model assuming that deep-sea carbonate dissolution occurs kinetically during the glacial event, forcing the seawater saturation index to 1. This corresponds to the maximum impact of carbonate dissolution on seawater chemistry, since dissolution stops when the saturation index reaches 1. We found that alkalinity rose to 18 equivalents/m³ when CO₂ reached 0.29 bar, as a result of Ca²⁺ release by the carbonate dissolution and Ca²⁺/Mg²⁺ exchange at mid oceanic ridges (Elderfield and Schultz, 1996). Seawater pH is calculated as 6.17, close to the 5.9 calculated in the absence of seafloor sedimentary carbonates. However, oceanic crust dissolution is sensitive to small pH changes. We calculate that, accounting for the maximum effect of carbonate buffering and because basaltic seafloor weathering is reduced, 0.29 bar of CO₂ can be reached in 19 m.y. This maximum effect corresponds to the dissolution of a 3-m-thick sedimentary carbonate layer spread over the global seafloor surface.

A possible anoxic state of the deep ocean may increase the pH of seawater during the snowball Earth. This would also result in decreased efficiency of the seafloor weathering. All along the glacial state, sulfate reduction may have produced alkalinity if the deep ocean was anoxic. The importance of this process depends upon the SO₄²⁻ level in seawater prior to the glacial event, and upon the availability of organic carbon in the water column and on the deep seafloor. A recent estimation shows that the SO₄²⁻ level might have been an order of magnitude below the present-day level (Hurtgen et al., 2002), thus largely below the dissolved inorganic carbon content of a snowball ocean, itself an order of magnitude above its present-day level in our calculations due to efficient atmospheric carbon dissolution.

For that reason, we suggest that the ocean pH was still largely controlled by carbonate speciation, even during a snowball event.

## CONCLUSION

If the snowball Earth occurred, it must have deglaciated. Our study shows that if the area of surface water in contact with the deep ocean is sufficient, the dissolution of the basaltic oceanic crust during the snowball glaciation may delay the atmospheric CO₂ rise required to reach the melting threshold, up to unacceptable time spans. The snowball theory thus requires either an extremely efficient dissolution of a large seafloor sedimentary carbonate reservoir whose existence remains hypothetical, or additional plausible processes able to decrease the CO₂ threshold through, for example, an albedo decrease, such as the deposition of dark fine volcanic ash on the ice and snow as suggested by Pierrehumbert (2004), or the existence of a thin equatorial ice (Goodman, 2006; Pollard and Kasting, 2005, 2006; Warren and Brandt, 2006).

## APPENDIX: GEOCHEMICAL MODELING PROCEDURE

To simulate the evolution of the carbon content before the global glaciation, we first run the model according to Proterozoic preglacial conditions until a steady state is reached (Fig. DR3; see footnote 1). We use a standard silicate weathering law, where Ca²⁺ and Mg²⁺ weathering rates are proportional to global runoff and depend on mean global air temperature (Oliva et al., 2004). Continental carbonate weathering is proportional to global runoff times the concentration of Ca²⁺ in equilibrium with calcite at the calculated mean air temperature and atmospheric CO₂. Both air temperature and runoff are assumed to be a function of atmospheric CO₂. Biological productivity in the surface ocean is made dependent on the input flux of phosphorus from the deep ocean (to simulate upwelling consequences) and from the weathering of continental surfaces. These dependencies are used for the calibration procedure, because continental silicate, carbonate, and phosphorus weathering are shut down during the snowball glaciation.

To apply snowball conditions we instantaneously shut down all continental weathering fluxes to reproduce the consequences of the rapid growth of continental ice sheets and the strong weakening of the hydrological cycle. Surface water temperature is forced to drop from 20 °C to 2 °C, and the vertical mixing between the surface and deep ocean decreases from 21 × 10⁶ m³/s to 0.2 × 10⁶ m³/s. The model then evolves freely.

The α value used in equation 1 depends on wind speed (Wanninkhof, 1992). The wind

action may be more efficient on one large surface than on several small. Since the total of ice-free surface is the sum of several small surfaces, K may be overestimated in our scenario. However in the extreme case wherein K is reduced by 10, the threshold surface required for an efficient exchange reaches 3 × 10⁴ instead of 3 × 10³ km², which remains very small. Note, however, that wind speed during the snowball event remains largely unknown, and the K we used is probably one of the best estimates that can be made.

## Seafloor Weathering Parameterization

The deep ocean reservoir is in contact with a water reservoir incorporated into the oceanic crust. The temperature of the deep ocean is held constant at 275 K, and based on a present-day estimation (Elderfield and Schultz, 1996), the exchange water flux between the reservoir of percolating waters and the deep ocean is 0.2 × 10⁶ m³/s, given that 735 × 10¹³ kg/yr of water are involved in the water circulation off axis at low temperature.

The reservoir of percolating waters is depleted in carbon through low-temperature dissolution of the basaltic crust, under the corrosive action of seawater. The release of alkalinity through this process is assumed to finally lead to the precipitation of carbonate minerals from highly saturated waters in veins, thus removing carbon from the ocean-atmosphere system at the million year time scale. The global sink of carbon is estimated as ~1.5–2.4 × 10¹² mol/yr (Alt and Teagle, 1999). Dissolution of seafloor basalts is calculated using laboratory kinetic laws to estimate the response of the H⁺, OH⁻, and H₂O promoted dissolution to change in seawater pH. The mineralogical composition of the crust is taken as the mineralogy of a tholeiitic basalt, containing 54% labradorite, 32% diopside, 9% basaltic glass, 4.3% apatite, and 1.3% forsterite. Kinetic constants are taken from the WITCH model database (Godderis et al., 2006; Drever, 1997). For each mineral, the overall dissolution rate R basalt (in moles of Ca²⁺ released per m² of mineral per year) is estimated as following:

$$R_{\text{basalt}} = \sum_{i} d_{i} \exp\left(\frac{-E_{i}}{R \cdot T_{p}}\right) C_{i}^{n_{i}}, \quad (\text{A1})$$

where the sum extends up to all the dissolving species (H⁺, OH⁻ and H₂O), while the index i stands for the dissolving species. Eᵢ and dᵢ are, respectively, the activation energy and dissolution constant for each mineral depending on the species i promoted dissolution. Tₚ is the temperature of the percolating waters into the oceanic crust at which dissolution occurs, here fixed at 313 K. The value Cᵢ stands for the concentration in the species i that promotes dissolution, while nᵢ is the order of the dissolution reaction. We assume that carbon consumption
through seafloor weathering (hereafter sfw) equals half of the alkalinity released through dissolution. The calculated overall CO₂ consumption flux F<sub>sfw</sub> is multiplied by a constant area, representing the total mineral surface area exposed to dissolution:

$$F_{\text{sfw}} = \text{SR} \cdot \text{area} \cdot R_{\text{basalt}}. \quad \text{(A2)}$$

This constant is in fact a calibration constant that was fixed so that the total CO₂ consumption rate through seafloor weathering equals 1.6 × 10¹² mol/yr under present-day conditions, within the lower range of recent estimations (Alt and Teagle, 1999; Spivack and Staudigel, 1994). SR is the seafloor spreading, normalized to its present-day value. It is kept constant at 1 in all simulations, except when assuming higher outgassing rate, itself corresponding to a higher seafloor-spreading rate. This overall kinetic formalism predicts a minimal basaltic dissolution rate at the present value of the ocean pH of 8.2.

## ACKNOWLEDGMENTS

This work was supported by the ECLIPSE program (Environnement et Climat du Passé: Histoire et Evolution) of the Institut National des Sciences de l'Univers (INSU) of the French Centre National de la Recherche Scientifique (CNRS). We thank G. Halverson and P. Regnier for comments and corrections of this manuscript, and P. Sepulchre for his help and encouraging comments. Constructive comments by S. Warren, P. Hoffman, and anonymous reviewers helped to improve the manuscript.

## REFERENCES CITED

Alt, J.C., and Teagle, D.A.H., 1999, The uptake of carbon during alteration of ocean crust: Geochimica et Cosmochimica Acta, v. 63, p. 1527-1535, doi: 10.1016/S0016-7037(99)00123-4.

Bodiselitsch, B., Koeberl, C., Master, S., and Reimold, W.U., 2005, Estimating duration and intensity of Neoproterozoic snowball glaciations from Ir anomalies: Science, v. 308, p. 239-242, doi: 10.1126/science.1104657.

Brady, P.V., and Gislason, S.R., 1997, Seafloor weathering controls on atmospheric CO₂ and global climate: Geochimica et Cosmochimica Acta, v. 61, p. 965-973, doi: 10.1016/S0016-7037(96)00385-7.

Caldeira, K., and Kasting, J.F., 1992, Susceptibility of the early Earth to irreversible glaciation caused by carbon dioxide clouds: Nature, v. 359, p. 226-228, doi: 10.1038/359226a0.

Cogne, J.P., and Humler, E., 2004, Temporal variation of oceanic spreading and crustal production rates during the last 180 My: Earth and Planetary Science Letters, v. 227, p. 427-439, doi: 10.1016/j.epsl.2004.09.002.

Dessert, C., Dupre, B., Gaillardet, J., Francois, L.M., and Allegre, C.J., 2003, Basalt weathering laws and the impact of basalt weathering on the global carbon cycle: Chemical Geology, v. 202, p. 257-273, doi: 10.1016/j.chemgeo.2002.10.001.

Donnadieu, Y., Ramstein, G., Fluteau, F., Besse, J., and Meert, J., 2002, Is high obliquity a plausible cause for Neoproterozoic glaciations?: Geophysical Research Letters, v. 29, 23, doi: 10.1029/2002GL015902.

Donnadieu, Y., Godderis, Y., Ramstein, G., Nedelec, A., and Meert, J., 2004a, A 'snowball Earth' climate triggered by continental break-up through changes in runoff: Nature, v. 428, p. 303-306, doi: 10.1038/nature02408.

Donnadieu, Y., Ramstein, G., Fluteau, F., Roche, D., and Ganopolski, A., 2004b, The impact of atmospheric and oceanic heat transports on the sea-ice-albedo instability during the Neoproterozoic: Climate Dynamics, v. 22, p. 293-306, doi: 10.1007/s00382-003-0378-5.

Drever, J.I., 1997, The geochemistry of natural waters (third edition): Englewood Cliffs, New Jersey, Prentice-Hall, Inc., 436 p.

Elderfield, H., and Schultz, A., 1996, Mid-ocean ridge hydrothermal fluxes and the chemical composition of the ocean: Annual Review of Earth and Planetary Sciences, v. 24, p. 191-224, doi: 10.1146/annurev.earth.24.1.191.

Etcheto, J., Boutin, J., and Merlivat, L., 1991, Seasonal-variation of the CO₂ exchange coefficient over the global ocean using satellite windspeed measurements: Tellus, ser. B, Chemical and Physical Meteorology, v. 43, p. 247-255, doi: 10.1034/j.1600-0889.1991.00017.x.

Gaffin, S., 1987, Ridge volume dependence on seafloor generation rate and inversion using long-term sealevel change: American Journal of Science, v. 287, p. 596-611.

Gaillardet, J., Dupre, B., Louvat, P., and Allegre, C.J., 1999, Global silicate weathering and CO₂ consumption rates deduced from the chemistry of large rivers: Chemical Geology, v. 159, p. 3-30, doi: 10.1016/S0009-2541(99)00031-5.

Godderis, Y., Francois, L.M., Probst, A., Schott, J., Moncoulon, D., Labat, D., and Viville, D., 2006, Modelling weathering processes at the catchment scale: The WITCH numerical model: Geochimica et Cosmochimica Acta, v. 70, p. 1128-1147, doi: 10.1016/j.gca.2005.11.018.

Goodman, J.C., 2006, Through thick and thin: Marine and meteoric ice in a "Snowball Earth" climate: Geophysical Research Letters, v. 33, p. L16701, doi: 10.1029/2006GL026840.

Higgins, J.A., and Schrag, D.P., 2003, Aftermath of a snowball Earth: Geochemistry, Geophysics, Geosystems, v. 4, doi: 10.1029/2002GC000403.

Hoffman, P.F., and Maloof, A.C., 1999, Glaciation: The snowball theory still holds water: Nature, v. 397, p. 384, doi: 10.1038/17006.

Hoffman, P.F., Kaufman, A.J., Halverson, G.P., and Schrag, D.P., 1998, A Neoproterozoic snowball earth: Science, v. 281, p. 1342-1346, doi: 10.1126/science.281.5381.1342.

Hoffman, P.F., Halverson, G.P., Domack, E.W., Husson, M.H., Higgins, J.A., and Schrag, D.P., 2007, Are basal Ediacaran (635 Ma) post-glacial "cap-dolostones" diachronous?: Earth and Planetary Science Letters, v. 258, p. 114-131.

Hurtgen, M.T., Arthur, M.A., Suits, N.S., and Kaufman, A.J., 2002, The sulfur isotopic composition of Neoproterozoic seawater sulfate: Implications for a snowball Earth?: Earth and Planetary Science Letters, v. 203, p. 413-429, doi: 10.1016/S0012-821X(02)00804-X.

Kirschvink, J.L., 1992, Late Proterozoic low-latitude glaciation: The snowball earth, in Schopf, J.W., and Klein, C., eds., The Proterozoic biosphere: Cambridge, Cambridge University Press, p. 51-52.

Oliva, P., Dupre, B., Martin, F., and Viers, J., 2004, The role of trace minerals in chemical weathering in a high-elevation granitic watershed (Estibere, France): Chemical and mineralogical evidence: Geochimica et Cosmochimica Acta, v. 68, p. 2223-2243, doi: 10.1016/j.gca.2003.10.043.

Pierrehumbert, R.T., 2004, High levels of atmospheric carbon dioxide necessary for the termination of global glaciation: Nature, v. 429, p. 646-649, doi: 10.1038/nature02640.

Pierrehumbert, R.T., 2005, Climate dynamics of a hard Snowball Earth: Journal of Geophysical Research, v. 110, D01111, doi: 10.1029/2004JD005162.

Pollard, D., and Kasting, J.F., 2005, Snowball Earth: A thin-ice solution with flowing sea glaciers: Journal of Geophysical Research-Oceans, v. 110, D07010, doi: 10.1029/2004JC002525.

Pollard, D., and Kasting, J.F., 2006, Snowball Earth: A thin-ice solution with flowing sea glaciers: Reply: Journal of Geophysical Research-Oceans, v. 111.

Schrag, D.P., Berner, R.A., Hoffman, P.F., and Halverson, G.P., 2002, On the initiation of a snowball Earth: Geochemistry, Geophysics, Geosystems, v. 3, doi: 10.1029/2001GC000219.

Spivack, A.J., and Staudigel, H., 1994, Low-temperature alteration of the upper oceanic crust and the alkalinity budget of seawater: Chemical Geology, v. 115, p. 239-247, doi: 10.1016/0009-2541(94)90189-9.

Swart, R., 1992, Facies analysis of Late Proterozoic carbonate turbidites in the Zerrissene Basin, Damara-Orogen, Namibia: Journal of African Earth Sciences, v. 14, p. 283-294, doi: 10.1016/0899-5362(92)90105-L.

Walker, L., Wilkinson, B.H., and Ivany, L.C., 2002, Continental drift and Phanerozoic carbonate accumulation in shallow-shelf and deep-marine settings: Journal of Geology, v. 110, p. 75-87, doi: 10.1086/324318.

Wanninkhof, R., 1992, Relationship between windspeed and gas-exchange over the ocean: Journal of Geophysical Research-Oceans, v. 97, p. 7373-7382.

Warren, S.G., and Brandt, R.E., 2006, Snowball earth: A thin-ice solution with flowing sea glaciers: Comment: Journal of Geophysical Research-Oceans, v. 111, C09016, doi:10.1029/2005JC003411.

Warren, S.G., Brandt, R.E., Grenfell, T.C., and McKay, C.P., 2002, Snowball Earth: Ice thickness on the tropical ocean: Journal of Geophysical Research-Oceans, v. 107, no. C10, 3167, doi: 10.1029/2001JC001123.

Zhang, R., Follows, M.J., Grotzinger, J.P., and Marshall, J., 2001, Could the Late Permian deep ocean have been anoxic?: Paleooceanography, v. 16, p. 317-329, doi: 10.1029/2000PA000522.

Manuscript received 29 May 2007

Revised manuscript received 30 August 2007

Manuscript accepted 31 August 2007

Printed in USA
# The snowball Earth aftermath: Exploring the limits of continental weathering processes

Guillaume Le Hir<sup>a</sup>, Yannick Donnadieu<sup>a</sup>, Yves Goddéris<sup>b,*</sup>, Raymond T. Pierrehumbert<sup>c</sup>, Galen P. Halverson<sup>d</sup>, Mélina Macouin<sup>b</sup>, Anne Nédélec<sup>b</sup>, Gilles Ramstein<sup>a</sup>

## ARTICLE INFO

Article history:

Received 22 May 2008

Received in revised form 10 October 2008

Accepted 10 November 2008

Available online xxxx

Editor: P. DeMenocal

Keywords:

snowball Earth

greenhouse

cap dolostones

weathering

modelling

## ABSTRACT

Carbonates capping Neoproterozoic glacial deposits contain peculiar sedimentological features and geochemical anomalies ascribed to extraordinary environmental conditions in the snowball Earth aftermath. It is commonly assumed that post-snowball climate dominated by CO₂ partial pressures several hundred times greater than modern levels, would be characterized by extreme temperatures, a vigorous hydrological cycle, and associated high continental weathering rates. However, the climate in the aftermath of a global glaciation has never been rigorously modelled. Here, we use a hierarchy of numerical models, from an atmospheric general circulation model to a mechanistic model describing continental weathering processes, to explore characteristics of the Earth system during the supergreenhouse climate following a snowball glaciation. These models suggest that the hydrological cycle intensifies only moderately in response to the elevated greenhouse. Indeed, constraints imposed by the surface energy budget sharply limit global mean evaporation once the temperature has warmed sufficiently that the evaporation approaches the total absorbed solar radiation. Even at 400 times the present day pressure of atmospheric CO₂, continental runoff is only 1.2 times the modern runoff. Under these conditions and accounting for the grinding of the continental surface by the ice sheet during the snowball event, the simulated maximum discharge of dissolved elements from continental weathering into the ocean is approximately 10 times greater than the modern flux. Consequently, it takes millions of years for the silicate weathering cycle to reduce post-snowball CO₂ levels to background Neoproterozoic levels. Regarding the origin of the cap dolostones, we show that continental weathering alone does not supply enough cations during the snowball melting phase to account for their observed volume.

© 2008 Elsevier B.V. All rights reserved.

## 1. Introduction

Geological evidence suggests the occurrence of two global glaciations during the Neoproterozoic. The older glaciation, the Sturtian event, is commonly inferred to have occurred at ca. 700–680 Ma, while the latter, end-Cryogenian glaciation ended at 635 Ma (Hoffman et al., 2004; Condon et al., 2005). The melting of these snowball Earth events requires very high atmospheric pCO₂, so that the greenhouse warming overcomes the freezing feedback induced by a high sea-ice albedo. The amount of CO₂ required to initiate melting of a snowball was first calculated to be 0.12 bar (400 PAL) using a simple energy balance model (Caldeira and Kasting, 1992) but more recent modelling suggests it was above 0.2 bar (660 PAL) (Pierrehumbert, 2004). Geochemical evidence seems to confirm exceptionally high pCO₂ following at least the 635 Ma snowball event (Kaseman

et al., 2005; Bao et al., 2008), but the behavior of this supergreenhouse climate under very high CO₂ levels is still poorly understood. Was it just a short pulse, lasting a few hundreds of kyrs, or was it longer-lived, with environmental consequences enduring for millions of years? According to the classic snowball Earth hypothesis the supergreenhouse drove a hyper-active hydrological cycle inducing intense continental weathering, such that CO₂ levels would have been drawn down from very high levels to preglacial levels in about 200 kyrs (Higgins and Schrag, 2003).

Up to now, this timespan subsequent to the snowball glacial event has never been investigated with a comprehensive climate model. In this contribution, four numerical models (Fig. 1) are used to investigate the evolution of the post glacial climate: (1) a 3D-climate model to characterize the super hothouse climate, (2) a continental weathering model to investigate the behaviour of the CO₂ weathering sink, (3) a simple carbon-alkalinity cycle model to simulate the evolution of the atmospheric CO₂ levels in the immediate aftermath of the snowball Earth, and (4) a simple eustatic model to investigate the

<footnote><sup>a</sup> LSCE, CNRS-CEA-UVSQ, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> LMTG, CNRS-Université de Toulouse, Observatoire Midi-Pyrénées, Toulouse, France</footnote>

<footnote><sup>c</sup> Dept of Geophysical Sciences, University of Chicago, Illinois, USA</footnote>

<footnote><sup>d</sup> School of Earth and Environmental Sciences, University of Adelaide, North Terrace, Adelaide, SA 5005, Australia</footnote>

<footnote>* Corresponding author.
E-mail address: godderis@lmtgobs-mip.fr (Y. Goddéris).</footnote>
<figure><img src="image_43.png" /><figcaption>**Fig. 1.** General design of the simulations. At each level, a question is added and addressed using a specific model. When going down through the diagram, any new specific question requires the coupling of a hierarchy of models (specified by the sign '+').</figcaption></figure>

consequences of post snowball continental weathering on cap dolostone accumulation. Our aim is to use numerical modelling to produce a first order reconstruction of the supergreenhouse climate and the timescale over which it is relaxed. Weathering processes are influenced by the air temperature controlling the dissolution rates of minerals, and by the continental runoff, which modulates dissolution through water-rock interactions and carries the weathering products to the oceans. We will first focus on climatic reconstructions of the supergreenhouse environment, after which the behaviour of the weathering sink under this climate will be discussed, and its consequences on the duration of the supergreenhouse will be investigated. Specifically, we will address the accumulation of cap dolostones which were deposited during the post glacial eustatic transgression and are a key element of the snowball hypothesis.

## 2. The climate of the snowball aftermath

### 2.1. Climate model description and simulation design

The climate of the post snowball is investigated with the general circulation climate model FOAM 1.5. The atmospheric component of FOAM is a parallelized version of NCAR's Community Climate Model 2 (CCM2) with the upgraded radiative and hydrologic physics incorporated in CCM3 v. 3.2 (Jacob, 1997). All simulations have been performed with an R15 spectral resolution (4.5°×7.5°) and 18 vertical levels. FOAM is used in mixed-layer mode, meaning the atmospheric model is linked to a 50-meter mixed-layer ocean with heat transport parameterized through diffusion. FOAM has been extensively applied to the study of ancient climate (Poulsen, 2003; Poulsen and Jacob, 2004; Donnadieu et al., 2006) and Pierrehumbert has demonstrated the validity of the FOAM radiative code in highly enriched CO₂ atmospheres such as the one hypothesized in the snowball Earth aftermath (Pierrehumbert, 2004). We have used a simplified paleogeographic model for the late Cryogenian (i.e. Marinoan) adapted from Macouin [pers.com] wherein the continental fragments of the former supercontinent Rodinia are clustered in low latitudes. An orography is assumed considering geological features (old orogenic belts and recent rifting) (Meert and Torsvik, 2003) [Macouin, pers. com.] (Fig. 2). We prescribe the presence of orogens and rifts (Meert and Torsvik, 2003) as elevated areas, between 500 and 1100 m (these relatively low altitudes are the consequence of the spectral smoothing implied by the model low surficial resolution), while the continental area outside uplifted surface has an elevation of 75 m. Solar luminosity is fixed at 94% of its present day value (1286 W/m²) according to stellar evolution models (Gough, 1981). Orbital parameters are set at their present day values. Regarding continental runoff, a present day simulation of FOAM linked to the mixed-layer ocean predict a global value of 27 cm/yr, a value very close to the observed 25 cm/yr (Berner and Berner, 1987) or 22 cm/yr (Gaillardet et al., 1999).

The melting phase of the snowball Earth has never been investigated. Indeed, modelling the climate under several hundred times the present day CO₂ pressure and with the presence of melting icecaps but also of melting sea-ice glaciers remains an unsolved challenge in terms of numerical tools and of computational times. The aim of this paper is not to model in detail the deglaciation: rather the

<figure><img src="image_44.png" /><figcaption>**Fig. 2.** End-Cryogenian (ca. 635 Ma) paleogeography used in the numerical modelling in this paper. Location of continental flood basalt provinces are shown in dark grey.</figcaption></figure>
focus here is on the supergreenhouse climate occurring just after the melting of the continental and sea ice. Considering the Neoproterozoic solar weakness and nonlinear relation between CO₂ concentration and the resultant greenhouse forcing, the following range of atmospheric CO₂ mixing ratios (400, 200, 100, 50, 25, 20, 15, and 10 PAL) has been used to model the influence of CO₂ concentrations on post-glacial temperature and runoff. The simulations reach equilibrium with the prescribed boundary conditions after 30 years. The analyses are done using the average of the next ten years. All the continents are assumed to be covered with a rocky regolith (visible albedo 0.24, infrared albedo 0.4, roughness 0.05 m, potential evaporation 0.1)

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

<figure><img src="image_46.png" /><figcaption>**Fig. 3.** a) Mean continental runoff and mean continental temperature calculated by the FOAM GCM under Neoproterozoic conditions and for various atmospheric CO₂ levels. The calculated present day values are also shown b) Energy fluxes calculated by the FOAM GCM under Neoproterozoic conditions and for various atmospheric CO₂ levels. *L* is the latent heat, *F*<sub>sens</sub> the sensible heat, IR the surface infra-red cooling and *S*<sub>abs</sub> the total incoming solar radiation adsorbed by the surface.</figcaption></figure>

<table><caption>**Table 1**
Water vapor mixing ratio at 5 different CO₂ levels for Neoproterozoic conditions, 1 PAL equals 280 ppm</caption><thead><tr><th>Atmospheric pCO<sub>2</sub> (PAL)</th><th>q<sub>s</sub></th><th>q</th><th>P (mm/mo)</th><th>Land ratio</th></tr></thead><tbody><tr><td>400</td><td>0.064</td><td>0.052</td><td>141</td><td>0.88</td></tr><tr><td>200</td><td>0.042</td><td>0.035</td><td>127</td><td>0.96</td></tr><tr><td>100</td><td>0.029</td><td>0.024</td><td>113</td><td>0.94</td></tr><tr><td>50</td><td>0.022</td><td>0.018</td><td>102</td><td>0.88</td></tr><tr><td>10</td><td>0.014</td><td>0.011</td><td>85</td><td>0.825</td></tr></tbody></table>

*q*<sub>s</sub> is the saturation water vapor mixing ratio at surface air temperature, in kg water per kg air. *q* is the actual water mixing ratio of the air, in kg water per kg air. *P* is the global mean precipitation in mm/month, and the “land ratio” is the ratio of the mean precipitation rate over land to the global mean precipitation rate.

## 2.2. Post snowball Earth climate: runoff and water cycle behaviour

Fig. 3a shows average terrestrial temperature and runoff evolution as a function of the atmospheric CO₂ mixing ratio. The mean continental temperature exhibits a roughly logarithmic response to CO₂ rise while runoff follows two pathways in our runs, a first one characterized by a logarithmic increase and a second one, above 100 PAL, characterized by a runoff saturating around 30–33 cm/yr which corresponds to a maximum increase of 22% when compared to the simulated runoff under present-day climate. This modest increase in runoff despite vastly higher pCO₂ sharply contrasts with the 7 fold increase in runoff assumed in the only other attempt to model silicate weathering rates under a post-snowball climate. In this previous study (Kump pers. com. in Higgins and Schrag, 2003), maximum runoff increase was estimated based on the inventory of the maximum amplitude of a variety of hydrological factors (Holland, 1978). As will be discussed further in Section 3, our result has important consequences for the rates of chemical weathering and geochemical evolution of the oceans following a snowball glaciation.

The weak increase in precipitation (and hence in runoff) with temperature in the post-Snowball world plays a key role in our results, and so merits some further discussion. In fact, a very similar phenomenon, although smaller, has attracted a lot of attention with regard to the effect of anthropogenic global warming on precipitation (notably Held and Soden, 2006). In particular, it has been found in virtually all climate simulations that increasing temperature by doubling or quadrupling CO₂ increases the precipitation much less than one would expect from the increase in saturation vapor pressure obtained from the Clausius-Clapeyron relation. Specifically, the Clausius-Clapeyron relation predicts that precipitation should increase by about 7% per degree of warming, whereas the actual figure for moderate warming is more like 2%. The presumed explanation for this disparity is that as the planet warms, atmospheric circulation becomes more sluggish, so that even though there is more water in the system, the rate at which it condenses into rain does not increase much (Vecchi et al., 2006). Atmospheric circulation weakens with the increase of atmospheric CO₂ levels in our simulations, in particular, the mean annual stream function zonally averaged is divided by 5 between the 20 PAL and the 400 PAL runs (figure not shown).

It has often been thought that infrared radiative cooling places an *a priori* constraint on precipitation, given that the dominant energy balance in the atmosphere is between latent heat release and radiative cooling (see Allen and Ingram, 2002 for a recent review of such arguments). However, Pierrehumbert (Pierrehumbert, 1999) showed that the argument is fallacious, or at least incomplete, since in principle the surface conditions can adjust to allow whatever column-integrated infrared cooling is necessary to accommodate the evaporation required to sustain a given precipitation level. In order to understand constraints on precipitation, a detailed consideration of the surface energy budget is unavoidable. In particular, it was argued
that the evaporation from the surface is in effect limited by the supply of solar radiation absorbed at the surface, so that if one warms the climate by increasing CO₂ without changing the solar radiation, there comes a threshold above which all solar radiation is used to sustain evaporation, whereafter precipitation no longer increases with temperature. This argument was extended and made more quantitative in Pierrehumbert (2002), where it was pointed out that the evaporation can somewhat exceed the solar absorption, but that the excess tends to be limited by the properties of turbulent fluxes in a stable boundary layer. This limitation in the increase in precipitation rates with increased pCO₂ has important implications for silicate weathering, and it is these implications that we explore in the Snowball recovery simulations here.

Table 1 and Fig. 3b summarize key boundary layer and surface budget properties pertinent to the hydrological cycle. As CO₂ is increased from 10 PAL to 400 PAL, the low level saturation mixing ratio increases by a factor of 4.6, and the actual boundary layer mixing ratio increases by a similar factor. However, the global precipitation rate only increases by 60%, with a proportional increase over the continents. To understand this relationship, it is necessary to examine the components of the surface energy budget. The following equation describes the global energy budget of the Earth surface:

$$S_{\text{abs}} = L + F_{\text{sens}} + IR. \qquad (1)$$

where Sₐbs is the global solar radiation absorbed by the surface, L is the latent heat flux (evaporation), Fₐsens is the sensible heat flux and IR the net infrared cooling of the surface. The mass flux of water due to evaporation is proportional to L, which hence determines the global precipitation rate. The net infrared cooling rate of the surface is the upward flux σT⁴ₐsurface minus the infrared back-radiation from the atmosphere. To the extent that Fₐsens and IR balance some of the incoming solar absorption, they capture energy that would otherwise be used to sustain evaporation. Both Fₐsens and IR increase as the surface is made increasingly warmer than the overlying atmosphere; they tend to be small when the air and surface temperatures are close, and moreover as the temperature increases and the boundary layer becomes more moist the net infrared cooling falls to zero because the atmosphere radiates to the ground at the same temperature as the ground radiates to the atmosphere. Looking at the calculated energy budget in the Fig. 3b, we see that the solar absorption remains nearly constant as CO₂ is increased, but that the sensible and infrared cooling falls precipitously as the climate is made warmer. The slack is taken up by evaporation, but the total “extra” energy in the 10 PAL case is only 86 W/m². Indeed, we see that by the time we reach the 400 PAL case, evaporation balances 82% of the surface absorbed solar radiation.

Thus, the limited increase in runoff in these simulations is seen as a robust consequence of the behavior of the surface energy budget. Nevertheless, it should be reminded that the runoff depends also on how available water vapor is redistributed between land and ocean. In our particular case where most of the continents are located in the tropical regions, deep atmospheric convection as well as sea surface temperature distribution should be two major factors controlling the water transport between the ocean and the continent. While we are confident in the way FOAM simulates the former one, the latter one may suffer from uncertainties owing to the use of the mixed-layer version of the FOAM model. Nevertheless, although we acknowledge that the mixed-layer model may bias our results, we suspect that the global energetic arguments developed in the paper remain the first order process controlling the absolute value of the global runoff. The use of a fully coupled climate model is warranted in order to provide a definitive answer. Potential shortcomings of our results may also come from the cloud parameterizations and its behaviour in a very different climate than the modern one. One of the primary determinants of cloud forcing is the cloud water content. The FOAM cloud water parameterization and its limitations were described in Pierrehumbert

(2004). In particular, while the scale height of cloud water increases in a warming climate, allowing high clouds to become more optically thick in a warmer world, the base-level cloud water parameter is independent on temperature. Insofar as there is more supply of water vapor in a hotter world, it is not unreasonable to expect that clouds at all levels might have higher water vapor content. In that case, the present simulations may underestimate the cloud water content. To test the implications of this possibility, we performed another run with 400× present CO₂, but with the cloud water parameter increased by a factor of 2. As a test in this run, both the cloud albedo effect and the cloud greenhouse effect roughly doubled compared to the control case, but the net cloud forcing remained nearly the same, with a slight shift toward more net cooling in the subtropics and more warming at the poles. The higher albedo did, however, reduce the surface absorbed solar radiation from 167 W/m² in the 400 PAL control case to 150 W/m² in the thick-cloud case. In accordance, this also reduced the global mean precipitation rate from 141 mm/mo to 114 mm/mo. Thus, higher cloud water content further reduces the precipitation. Nevertheless, the factors governing cloud fraction are poorly understood. In particular, there is observational evidence that for warm rain processes, cloud fraction goes down with temperature, owing to increased precipitation efficiency (Rapp et al., 2005). This feedback is not present in the cloud fraction scheme we are using. This is a fruitful area for future inquiry.

## 3. Continental weathering during the supergreenhouse climate

We now consider the weathering pattern as a function of the simulated climate. As shown in Section 2, the simulated post-glacial period is warm and moderately humid, with maximum continental runoff limited at 33 cm/yr above 100 PAL CO₂. Depending on this climate regime, the silicate weathering acts as a restoring force of the carbon cycle that reduces atmospheric CO₂. Continental weathering is a function of climate, increasing with temperature and continental runoff as demonstrated by numerous field studies (Brady, 1991; White and Blum, 1995; Dessert et al., 2001; Oliva et al., 2003). Based on the compilation of the main cation fluxes exported by small basaltic and granitic catchments, Dessert et al. (2001) and Oliva et al. (2003) have proposed parametric relationships linking the climatic parameters (mean annual temperature and runoff) to CO₂ consumption by the dissolution of silicate minerals. These laws are strictly valid only for the present day environment and do not include a direct dependence of rock weathering on the acidity of rainwater (Oliva et al., 2003). Exporting these laws into the post-snowball environment is tricky for three reasons: (1) The climate is exotic and is beyond the validity field of the parametric laws. (2) The Earth surface might have been largely modified by the global glaciation. That is, dynamic glaciation persisting for several millions of years might have produced fine grained materials (i.e. rock flour) with very high reactive surface areas, promoting the dissolution of minerals and the atmospheric CO₂ consumption. (3) The very high levels of atmospheric CO₂ in the snowball Earth aftermath (100–400 PAL) would have acidified the rainfall. These CO₂ levels, which are potentially above the common soil CO₂ levels encountered today, have promoted the mineral dissolution rates. Because the use of parametric laws to describe continental weathering cannot account for these three major changes, it is necessary to investigate the runoff-temperature-weathering feedback with a mechanistic model of silicate rock weathering in natural environments.

### 3.1. Description of the continental weathering model

Here we utilize the WITCH weathering model (Goddéris et al., 2006), which was initially developed to simulate weathering processes in soil profiles and in river catchments above silicate bedrock and calibrated on a granitic watershed. It successfully
describes the seasonal evolution of the chemical composition of soil solutions as a function of water-rock-biosphere interactions, as well as the mean annual chemical composition of the main stream. For present day studies, WITCH is coupled to a biospheric model that calculates the water circulation in the soil, the partial pressure of CO₂ below ground level as a result of root respiration and organic carbon decay, and the uptake and release of cations by the living and dead biomass. WITCH calculates the dissolution of minerals as a function of the soil solution composition at each time step through kinetic laws derived from transition state theory concept (Eyring, 1935). The kinetic and thermodynamic parameters have been experimentally determined. Four dissolving species are accounted for (H⁺, OH⁻, H₂O and organic ligands, but hereafter the organic ligands dissolution will be neglected due to absence of terrestrial vegetation). The dissolution rate of a given mineral R_s is thus equals to:

$$R_s = A \cdot S_w \cdot \left[ \sum_i k_i \exp\left(\frac{-E_a^i}{RT}\right) \cdot a_i^{n_i} \cdot f_{\text{inh}} \right] (1 - \Omega^s) \quad (2)$$

where A stands for the mineral reactive surface (m²/m³ of soil), S_w is the soil moisture saturation (function of the water volumetric content, Sverdrup and Warfinge, 1995), k_i and E_a^i are the dissolution rate constant and the activation energy for the i species-promoted dissolution, a_i and n_i stand for the activity of species i and the order of reaction with respect to the dissolving species i, respectively. f_inh describes the inhibiting effect of aqueous species on mineral dissolution (Oelkers et al., 1994; Schott and Oelkers, 1995; Devidal et al., 1997). Finally, the effect of departure from equilibrium on mineral dissolution/precipitation rates is accounted for by the factor (1 - Ω^s), where Ω^s stands for the solution saturation index with respect to the particular solid and s is a stoichiometric number.

The WITCH model includes 32 silicate minerals. The number of soil layers that are modelled can be adapted to any specific problem. The mineral reactive surface is calculated according to a simple parametric law assumed to be an approximation of the solid BET surface (Brunauer et al., 1938; Sverdrup and Warfinge, 1995):

$$A = (8.0x_{\text{clay}} + 2.2x_{\text{silt}} + 0.3x_{\text{sand}} + 0x_{\text{coarse}}) \cdot \rho \quad (3)$$

where x\_clay, x\_silt, x\_sand and x\_coarse represent the textural fractions of the clay, silt, sand and coarse materials in the soil, such that x\_clay + x\_silt + x\_sand + x\_coarse = 1. ρ is the soil density in g/m³. We also include continental calcite and dolomite dissolution in the WITCH model. Because carbonate dissolution is several orders of magnitude faster than silicate dissolution, we assume that continental waters are at equilibrium respectively with the carbonate phases when their abundance exceeds 2% by volume of the bedrock (thus when carbonate rock is the dominant phase, geochemically speaking), and with the CO₂ partial pressure below ground level. The solubility product of dolomite is calculated as a function of the air temperature T (Drever, 1997):

$$K_{\text{dol}} = 10^{-17.09} \exp \left[ - \frac{39480}{R} \cdot \left( \frac{1}{298.15} - \frac{1}{T} \right) \right] \quad (4)$$

where the standard enthalpy of reaction is in J/mol and R is the perfect gas constant.

For calcite, the solubility product is given by Drever (1997):

$$K_{\text{cal}} = 10^{-8.48} \exp \left[ - \frac{9610}{R} \cdot \left( \frac{1}{298.15} - \frac{1}{T} \right) \right] \quad (5)$$

## 3.2. Design of the WITCH simulations

A number of parameters have to be fixed for using the WITCH model in the Neoproterozoic context. (1) The dissolution of minerals

by organic ligands, and exchanges of elements with the living and dead biomass are set to zero because of the absence of continental vegetation at this time. (2) Since we aim at constraining the maximum CO₂ consumption in the post snowball environment, we assume granitic lithology covering all continental surface except where the occurrence of basaltic outcrops can be reasonably inferred from geological data (see compilation by Goddéris et al., 2003 and Fig. 2). Mineralogy of the granite and basalt used in WITCH are shown on Table 2. Carbonate weathering does not consume CO₂ on a 10⁵ to 10⁶ year timescale. Its role during the melting phase will be discussed in Section 5 (3) The volumetric content of water in the weathering profile has to be fixed. This parameter is important as it defines the fraction of the total reactive mineral surface that is wet and weatherable. In a first step, we fix it to a standard value of 0.2 m³ water/m³ soil. In a second step, we use the soil water content simulated by FOAM. (4) Given the assumed absence of soil biomass, CO₂ level is fixed in pore spaces to the atmospheric concentration. The WITCH model does not include a carbon budget of the below ground solutions so far. This means that the model is providing unlimited CO₂ to the water solutions, resulting in a possible slight overestimation of the weathering rates. (5) Assumptions about the thickness of the weathering profile and its texture are required in the model. We have chosen a shallow two layer structure with a thickness of 25 cm for the surface layer and of 100 cm for the deepest one. These relatively shallow thicknesses have been chosen based on the influence of persistent mechanical weathering during glaciation (Donnadieu et al., 2003), which would not allow the survival of any pre-glacial thick soil profiles. The soil system is assumed to have been largely reset by the glacial event. A sensitivity test to soil thickness is performed below. (6) With regards to the grain size, three tests have been run in which the reactive surface is changed according to Eq. (3) and intended to span the range of materials observed on the present earth surface. (7) temperature of the weathering profiles is taken from the GCM output field for the various atmospheric CO₂, and the GCM runoff is assumed to be equal to the vertical drainage of the weathering profile for each continental grid element. This assumption is valid because of the large size of the continental pixels in this study: all the water flowing through the weathering profiles can be assumed to end up in a river draining the pixel.

## 3.3. Weathering rates in a supergreenhouse environment

Under 400 PAL, calculated CO₂ consumption by silicate weathering (Fig. 4a and b) ranges from 2.4 to 28 times the present-day flux, as estimated by Gaillardet et al. (1999), with the lowest value calculated assuming a coarse grain size (sand) material with reactive surface of 0.5 × 10⁶ m²/m³ (LR run, Table 3). While such a coarse-grained surface layer on the continents is unlikely due to grinding by glaciers, this test provides a useful absolutely minimum level of silicate weathering to be expected. To test the sensitivity of weathering rates as a function of the mineral reactive surface, a finer grain size is used (MR run, Table 3), which yields a 6-fold increase over present day CO₂ consumption (Fig. 5a), showing that the mineral reactive surface is a crucial factor governing the post glacial weathering. In order to represent the Earth surface resulting from a global ice cover more precisely, we apply the reactive surface value measured on loesses derived from glacial erosion during the Last Glacial Maximum (Kahle et al., 2002) (13 × 10⁶ m²/m³, HR run, Table 3) to the full 1.25 m-thick soil profile. In this run, which yields what we regard as a maximum value for weathering rates, the CO₂ consumption rate is 10 times

<table><caption>**Table 2**
Standard mineralogy for the silicated lithologies (in percent)</caption><tr><td>Granite</td><td>Quartz 31%</td><td>Orthose 26.5%</td><td>Andesine 42%</td><td>Apatite 0.3%</td></tr><tr><td>Basalt</td><td>Labradorite 45%</td><td>Basaltic glass 17.7%</td><td>Diopside 34%</td><td>Forsterite 2%</td><td>Apatite 1.3%</td></tr></table>
greater than the present-day value (Fig. 4a). A second run on a 6 m thick loess layer over the whole continental surface (Fig. 4a, run HR/thick) results in an increase in the calculated CO₂ consumption by a factor of only 2 to 3, which shows that predicted weathering rates are a non-linear function of loess thickness. This relationship is due to the limitation in weathering rates by the increase in the saturation state of pore water with respect to the primary silicate minerals as it percolates through the soil horizon. The result is significant because

continental runoff is not commensurate with increased CO₂ levels, meaning soil waters rapidly reach oversaturation with respect to primary minerals.

Predicted CO₂ consumption by silicate weathering at low CO₂ levels (10 PAL) are all at least an order of magnitude below the present value of 11.7×10¹² mol CO₂/yr (Gaillardet et al., 1999). These low predicted values are mainly driven by (1) the low below ground CO₂ levels (fixed at 10 PAL, which can be regarded as a lower hand value for

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

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

<figure><img src="image_49.png" /><figcaption>**Fig. 4.** a) Total CO₂ consumption by continental silicate weathering at various pCO₂ levels normalized to the present day value, as calculated by the WITCH model forced with the climatic output of the FOAM GCM. HR: high reactivity experiment, MR: middle reactivity experiment, LR: low reactivity experiment, HR/thick: high reactivity experiment with a 6 m-thick loess, and MR/FOAM: middle reactivity experiment with the water volumetric content of the loess calculated by FOAM. b) Calculated pattern of the CO₂ consumption by continental silicate weathering, at 400 PAL by the WITCH model(MR simulation, in 10¹² mol/yr). Basaltic outcrop are identified by a higher flux.</figcaption></figure>
<table><caption>**Table 3**
WITCH simulations</caption><thead><tr><th>WITCH simulation set</th><th>Temperature and runoff</th><th>Loess reactive surface (upper layer)</th><th>Water volumetric content of the loess</th></tr></thead><tbody><tr><td>HR</td><td>FOAM</td><td>13 × 10<sup>6</sup> m<sup>2</sup>/m<sup>3</sup></td><td>20%</td></tr><tr><td>MR</td><td>FOAM</td><td>6.2 × 10<sup>6</sup> m<sup>2</sup>/m<sup>3</sup></td><td>20%</td></tr><tr><td>LR</td><td>FOAM</td><td>0.51 × 10<sup>6</sup> m<sup>2</sup>/m<sup>3</sup></td><td>20%</td></tr><tr><td>MR/FOAM</td><td>FOAM</td><td>6.2 × 10<sup>6</sup> m<sup>2</sup>/m<sup>3</sup></td><td>FOAM</td></tr></tbody></table>

All sets of run are using the FOAM climatology. Reactive surface of the upper layer (0.25 m thick) for the MR run (LR run) is calculated according to Eq. (3), assuming a grain size distribution with 25% clay and 75% silt materials (50% sandy and 50% coarse materials). Reactive surface of the upper layer for the HR run is prescribed at $13 \times 10^6$ m$^2$/m$^3$ based on Ref. The lower layer (1 m thick) contains sandy (50%) and coarse (50%) materials in all runs (yielding a reactive surface of $0.51 \times 10^6$ m$^2$/m$^3$), except for the HR run where we assumed a mineral reactive surface of the lower layer of $1 \times 10^6$ m$^2$/m$^3$.

present day soils), (2) the absence of organic acid produced by the biotic activity, (3) the relative aridity of the continental surfaces, (4) the smaller total continental surfaces ($100 \times 10^6$ km$^2$ in the Neoproterozoic, compared to the $150 \times 10^6$ km$^2$ today).

## 4. Duration of the post-snowball supergreenhouse climate

Because the increase in CO$_2$ consumption by silicate weathering does not scale linearly with higher CO$_2$ levels, the post-snowball supergreenhouse climate should last much longer than the few hundred thousands years envisaged in the original snowball Earth model (Higgins and Schrag, 2003; Hoffman et al., 1998). Indeed, long term atmospheric CO$_2$ evolution is controlled by the balance between the solid Earth degassing and the total silicate weathering (Walker et al., 1981; Berner, 2004). No clear consensus exists concerning changes in the CO$_2$ outgassing flux in the geological past (Gaffin, 1987; Rowley, 2002; Cogné and Humler, 2004); thus we conservatively assume that the total solid earth degassing rate during the late Neoproterozoic was close to its present day value. Modelling the evolution of atmospheric CO$_2$ through geological times during a timespan where silicate weathering largely exceeds solid Earth degassing requires at least the calculation of the alkalinity and carbon budget of the surficial reservoirs (ocean + atmosphere), as well as the carbonate speciation in seawater. To perform these calculations, we use a simplified version of the COMBINE model (Goddéris and Joachimski, 2004), a box model with two oceanic reservoirs (surface and deep) and one for the atmosphere that has previously been applied to simulating the geochemical evolution of the snowball ocean (Le Hir et al., 2008a,b). CO$_2$ consumption by silicate weathering is taken from the FOAM/WITCH simulations, as a function of atmospheric CO$_2$. Since the WITCH model only gives cation and CO$_2$ consumption fluxes for the 8 solved climate scenarios (at 400, 200, 100, 50, 25, 20, 15 and 10 PAL CO$_2$), a linear interpolation for each continental grid cell (Fig. 2) is performed to define a continuous trend for cation fluxes as a function of atmospheric CO$_2$ between 400 and 10 PAL. Through this interpolation, we start from the CO$_2$ consumption fluxes calculated for each continental grid cell by the WITCH model, depending on the FOAM temperature and runoff, and end up with a simple trend describing the global CO$_2$ consumption through silicate rock weathering as a function of atmospheric CO$_2$. This interpolated flux is then incorporated at each timestep into the simplified COMBINE numerical model (Le Hir et al., 2008a,b). Since no GCM simulation was performed above 400 PAL our CO$_2$ history begins at 0.12 bar.

Fig. 5 shows the atmospheric CO$_2$ evolution in the aftermath of snowball Earth. Four shapes of global CO$_2$ consumption have been used to decipher the role of weathering parameters in controlling how atmospheric CO$_2$ would have decreased in the snowball aftermath. Using the FOAM water volumetric content (MR/FOAM run), the continental weathering becomes less efficient due to enhanced dryness of the weathering profiles. Indeed, the hothouse climate

plays an inhibiting role on weathering by dropping the soil water volumetric content below 20% which was the fixed value of the other experiments. Nevertheless, changes linked to the use of the FOAM water volumetric content remain minor. This is the result of the first order control played by the runoff on the weathering, bringing waters in contact with the minerals and transporting the weathering products. The soil water content only plays a secondary role on how long the runoff interacts with the minerals (i.e. residence time). Hence, the simulated CO$_2$ reductions are of the same order for the MR/FOAM run and the MR run. Larger effects are visible when changing the reactive surface (HR run) and even more when also changing the soil thickness (HR/thick run). After 1 million years, atmospheric CO$_2$ decreases to 100 PAL, 70 PAL and 20 PAL for MR, HR and HR/thick runs, respectively. Recovery timescales are thus dependent on the weathering profile parameters, though the HR run is probably the one better representing soil conditions that may have occurred after several million years of global glaciation. Hence, we can conclude that the recovery timescale from a snowball Earth is 1) on the order of a few million years in any case and 2) mostly dependent on the soil thickness, the reactive surface, and the runoff pattern.

## 5. The deglaciation and the cap dolostone deposition, a model-data discussion

Is it possible to explore the weathering rates during the melting phase itself? As discussed above, modelling the melting of a snowball Earth at very high CO$_2$ levels using a coupled ice-sheet – climate model remains a technical challenge for a while because of computational limitations. Our GCM cannot at present handle rapidly melting ice sheets and sea glaciers in the initial stages of the high CO$_2$ world. Nevertheless, guesses can be assumed regarding the continental runoff and temperature resulting from the melting phase. A

<figure><img src="image_50.png" /><figcaption>**Fig. 5.** Post-snowball atmospheric CO$_2$ evolution calculated by the simplified COMBINE geochemical module, using the weathering rates calculated by the WITCH mechanism weathering model in the four cases described above.</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>**Fig. 6.** a) Numerical model used to calculate the deposition rate of the cap dolostones. b) Cap dolostones thickness accumulated during the deglaciation phase (lasting 10 kyr), as a function of the level below the pre-glacial sealevel (assumed to be located at 0 m). 0.5 km represents the sea-level drop linked to Neoproterozic global glaciations. *L* denotes the thickness of the water layer intercepting the seafloor and from which cap dolostone are precipitated.</figcaption></figure>

previous modelling study coupling a climate and an ice-sheet models calculates a land-ice volume reaching ~200×10⁶ km³ for a global glaciation (Donnadieu et al., 2003), a value in close agreement with a glacio-eustatic sea level fall of 500 m or more deduced from the stratigraphic record of end-Cryogenian glaciation in Namibia (Hoffman et al., 2007). If the melting spans 10 kyrs (Hoffman et al., 2007) and if the icecap melts linearly with time, this would correspond to a runoff flux of 20 cm/yr during 10 kyr. The contribution of rainfall to

this runoff cannot be estimated accurately without a coupled GCM-ice sheet simulation of the melting of the snowball. This is a fruitful area for future inquiry. Indeed, during deglaciation, it is possible that precipitation was focused on low latitude continents owing to newly formed ice-free oceanic areas and to resulting steep meridional temperature gradient (Peyser and Poulsen, 2008), but this process cannot be explored here. Regarding temperature, we assume that the melt waters responsible for the dissolution of the minerals are around
1 °C. At such a low temperature, the solubility of carbonate minerals increases, but silicate weathering is strongly inhibited. Atmospheric CO₂ is fixed at 400 PAL. Under such conditions (Deglac run), the predicted silicate weathering flux obtained by WITCH only consumes 15.11 × 10¹² mol CO₂/yr despite very high CO₂ pressure acidifying the percolating solutions. This value is only slightly above the present day estimate (Gaillardet et al., 1999). The net atmospheric pCO₂ drawdown linked to silicate weathering itself is insignificant owing to the melting phase duration (10 kyr) and the high atmospheric CO₂ content. This result is in agreement with the observation that the maximum silicate weathering is reached after cap dolostone deposition, i.e. after the melting and transgressive phase, as shown by Sr and Ca isotope data (Kaseman et al., 2005; Nédélec et al., 2007; Nogueira et al., 2007).

While silicate weathering rates during the deglaciation would have been relatively low, carbonate weathering should have been much higher because the dissolution of carbonate minerals is enhanced by the low temperature. Higgins and Schrag (2003) argued that the carbonate weathering flux might be responsible for the accumulation of the 18 m-thick (on average) post-glacial cap dolostones (Hoffman et al., 2007) during the melting phase, mainly resulting from the weathering of emerged carbonate platforms and intra-cratonic carbonate outcrops. To test this hypothesis, an additional carbonate lithology is added to the model, randomly along the coast line so that the carbonate weathering can be investigated with WITCH. Indeed, given the low sea-level prior to deglaciation, large carbonate (calcite and dolomite) platforms would have been exposed to subaerial weathering. The exact location of the carbonate outcrops is not important in this numerical experiment, since global climate is assigned at a fixed value for each continental grid element (runoff=20 cm/yr, temperature=1 °C), whatever its location. Assuming that these carbonate platforms accumulate prior to the snowball Earth beneath a 100 m-thick epicontinental oceanic layer, their maximum areal extent cannot exceed 20 × 10⁶ km², assuming a present day oceanic hypsometry. In addition to these emerged carbonate shelves, ancient carbonates in the continental interior would also have been subjected to weathering. Today, intra-continental carbonate outcrops cover about 20 × 10⁶ km² of the land surface (Amiotte-Suchet et al., 2003). We use the same value for the Neoproterozoic continental surface. Thus, the assumed maximum carbonate surface area in the direct aftermath of the snowball (as long as sealevel is below the 100 m depth level) was 40 × 10⁶ km². In our model, 50% of these outcrops are assumed to be calcite and 50% dolomite and the water draining these areas are assumed to be saturated with respect to the relevant carbonate mineral. The WITCH simulation can be performed at very high CO₂ levels. Indeed we prescribe the climate for each continental grid element (constant runoff and temperature). The FOAM GCM output are not used here and the 400 PAL threshold is not a limitation for the WITCH model. Because it has been suggested that 400 PAL might be too low (Pierrehumbert, 2004) to melt the snowball Earth, we perform our simulations assuming that the CO₂ level was around 1000 PAL, thus further maximizing the continental weathering rates and consequently the cap-carbonates deposition rate.

Since the WITCH model calculates the Mg²⁺ and Ca²⁺ flux to the ocean, we can convert these fluxes in thickness of cap dolostones being deposited, assuming that chemical conditions for cap dolostones to precipitate are satisfied (Font et al., 2006), a deposition surface and a mean density of 2600 kg/m³ for the sedimentary carbonates. We deliberately choose this mass balance method since the kinetics of cap dolostone precipitation is not known. To be as realistic as possible, we build up a numerical model calculating the deposition surface as a function of the sealevel during the melting phase (the so-called "hypsometric model", Fig. 1). The depositional surface is estimated assuming the present day oceanic hypsometry, and that cap dolostones where deposited within a water layer of thickness L (Fig. 6a). This water layer rises from the slope along the shelf, as the oceanic volume rises linearly in response to the icecap melting, forcing sealevel to rise. The interception between the seafloor surface and this water layer defines

the area of cap dolostones accumulation. All Mg²⁺ originating from continental weathering is assumed to precipitate instantaneously on this area, assuming consequently that the chemical conditions required to precipitate dolostones are satisfied within the considered water layer. We also account for the decrease in the surface of emerged carbonate platforms once the sealevel reaches 100 m below its pre-perturbation value (see Fig. 6a).

With L=20 m, we found that, if the melting phase lasts 10 kyr, a maximum of 50 cm of cap dolostones could be accumulated by using all the Mg²⁺ from coeval continental weathering (Fig. 6b). The calculated thickness of cap-carbonates is not heavily dependent on the assumed thickness L of the water layer in which cap dolostones are deposited (Fig. 7), nor on the duration of the melting phase (since a longer duration implies a lower runoff, and hence a lower weathering). The calculated thickness of the cap dolostones is about 40 times lower than the mean compiled value of 18.5 m (Hoffman et al., 2007), although this is the mean arithmetic thickness, without weighting by the original depositional area. To reconcile the model output with the data would require that the surface of dolomitic outcrop being weathered was 40 times greater than the assumed 20×10⁶ km², which is clearly unrealistic. Conversely, if we allow that the area over which cap dolostones were precipitated was lower by a factor of 10, then cap dolostones on average 5 meters thick would have been deposited. While the base-Ediacarian cap dolostones in Namibia and Death Valley, for example, are very thick, the equivalent cap dolostones in Svalbard, East Greenland, northern Norway, Australia, and elsewhere are on average much thinner. If these thinner cap dolostones are in fact more typical and their areal extent restricted, then the carbonate weathering model for cap dolostone precipitation may in fact be valid. This also suggests that the thickness of cap dolostones was heterogeneous over continental shelves and that high average thicknesses of preserved cap dolostones are a preservative bias. In addition, as noted above, we do not take into account the runoff component resulting from the precipitation occurring during the deglaciation itself but the one resulting from the continental ice cover melting. Although an accurate estimate of the former one depends on coupled GCM ice-sheet simulations in this peculiar context, we can note that the runoff calculated by FOAM at 400 PAL saturates around 30–33 cm/yr. By adding both components (the one from the ice-sheet melting and the one from the ocean-atmosphere system), we obtain a mean thickness of 82 cm, a value which remains far for the one suggested by Hoffman et al. (Hoffman et al., 2007). Hence, even with the shortcomings of our model, our results still support the possibility that thinner cap dolostones are in fact more typical during the aftermath of the snowball Earth. This conclusion remains to be confirmed in the future.

Alternatively, the timescale for cap dolostone precipitation may be grossly underestimated. In fact, the timescale of cap dolostone precipitation is a subject of significant debate (Hoffman et al., 2007). Importantly, the stratigraphic context of the cap dolostone unambiguously constrains its timing of deposition to the period spanning the post-glacial eustatic rise in sea level. Rapid deposition of the cap dolostone (2000–10000 years; Hoffman et al., 2007) has thus been inferred based on modelled high rates of ice sheet decay at the end of the snowball (Hyde et al., 2000) and by analogy with Quaternary deglaciation (Bard et al., 1990). However, magnetostratigraphic data from multiple successions imply that cap dolostones span multiple magnetic reversal events (Trindade et al., 2003; Kilner et al., 2005), suggesting that the cap dolostones accumulated over several hundreds of thousands of years. Whether this timescale may be reconciled with the climate models is still an open question, pending on a careful modeling of the melting phase.

## 6. Conclusion

In this study, we force a numerical model predicting continental weathering with 3D-climate model output to investigate the climate
and weathering response during the supergreenhouse climate following the snowball Earth event. The main conclusion of this study is that post-glacial silicate weathering was likely not as vigorous as commonly assumed (Higgins and Schrag, 2003; Hoffman et al., 2007). This result is largely due to the non-linear relationship between average continental temperatures and runoff rates under very high CO₂ levels. This relative dryness of the supergreenhouse climate is linked to physical constraints governing the surface energy budget. Indeed the energy released by the evaporation cannot overcome the net solar energy absorbed in surface. Since the incoming solar energy is not enhanced by increasing CO₂, the evaporation rates, and thus the intensity of the hydrological cycle, level off at high CO₂ levels. For example, even at 400 times the present CO₂ pressure, continental runoff is only 22% higher than the calculated present day values. As a consequence, under a 400 PAL CO₂ climate, our simulated discharge of dissolved elements from continental weathering into the ocean shows at best a 10-fold increase compared to the present-day estimate. As a consequence, the time required to restore CO₂ concentrations to roughly pre-glacial levels is on the order of a few million years, rather than a few thousand years, as has been previously assumed (Higgins and Schrag, 2003). Hence, an unusually warm climate may have persisted for many million years, thus extending the duration of the environmental filter on the biospheric evolution. A more protracted timescale for the recovery from the post-snowball supergreenhouse should be testable using isotopic proxies (e.g. ⁸⁷Sr/⁸⁶Sr, ⁶²Mg) of the continental weathering.

A final important result of our linked numerical models is that, assuming a timescale of deglaciation of 10 kyr, total continental weathering (carbonate+silicate) supplies only enough Mg²⁺ to accumulate a maximum of 50 cm of cap dolostone over the continental shelves, which is over an order of magnitude lower than the calculate average global thickness of base-Ediacaran cap dolostones (Hoffman et al., 2007). Our history of continental rock dissolution in the aftermath of the snowball glaciation is also in agreement with the isotopic proxies for continental weathering that suggest that the major increase in silicate weathering occurred directly after cap dolostone deposition (Kaseman et al., 2005).

## Acknowledgments

This study was funded by the French National Agency for Research program ANR-06-BLAN-0347 and by the CNRS/INSU through the ECLIPSE program. We thank Praksys for the maintenance of the cluster of PC at the LMTG (www.praksys.org). We thank Lee Kump and an anonymous referee for their insightful reviews, and Peter deMenocal for the quality of the editorial procedure.

## References

Allen, M.R., Ingram, J., 2002. Constraints on future changes in climate and the hydrologic cycle. Nature 419, 224-232.

Amiotte-Suchet, P., Probst, J.L., Ludwig, W., 2003. World wide distribution of continental rock lithology: implications for atmospheric/soil CO₂ uptake by continental weathering and alkalinity river transport to the oceans. Glob. Biogeochem. Cycles 17, doi:10.1029/2002GB001891.

Bao, H.M., Lyons, J.R., Zhou, C., 2008. Triple oxygen isotope evidence for elevated CO₂ levels after a Neoproterozoic glaciation. Nature 453, 504-506.

Bard, E., Hamelin, B., Fairbanks, R.G., Zindler, A., 1990. Calibration of the ¹⁴C timescale over the past 30,000 years using mass spectrometric U–Th ages from Barbados corals. Nature 345, 405–410.

Berner, R.A., 2004. The Phanerozoic Carbon Cycle. Oxford University Press, New York. 150 pp.

Berner, E.K., Berner, R.A., 1987. The global water cycle: geochemistry and environment. Prentice Hall. 397 pp.

Brady, P.V., 1991. The effect of silicate weathering on global temperature and atmospheric CO₂. J. Geophys. Res. 96, 18101–18106.

Brunauer, S., Emmett, P.H., Teller, E., 1938. Adsorption of gases in multimolecular layers. J. Am. Chem. Soc. 60, 309–319.

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

Cogné, J.P., Humler, E., 2004. Temporal variation of oceanic spreading and crustal production rates during the last 180 My. Earth Planet. Sci. Lett. 227, 427–439.

Condon, D., Zhu, M.Y., Bowring, S., Wang, W., Yang, A.H., Jin, Y.G., 2005. U–Pb ages from the Neoproterozoic Doushantuo Formation, China. Science 308, 95–98.

Dessert, C., Dupré, B., François, L.M., Schott, J., Gaillardet, J., Chakrapani, G.J., Bajpai, S., 2001. Erosion of Deccan Traps determined by river geochemistry: impact on the global climate and the ⁸⁷Sr/⁸⁶Sr ratio of seawater. Earth Planet. Sci. Lett. 188, 459–474.

Devidal, J.L., Schott, J., Dandurand, J.L., 1997. An experimental study of kaolinite dissolution and precipitation kinetics as a function of chemical affinity and solution composition at 150 °C, 40 bars, and pH 2, 6.8, and 7.8. Geochim. Cosmochim. Acta 61, 5165–5186.

Donnadieu, Y., Fluteau, F., Ramstein, G., Ritz, C., Besse, J., 2003. Is there a conflict between the Neoproterozoic glacial deposits and the snowball Earth interpretation: an improved understanding with numerical modeling. Earth Planet. Sci. Lett. 208, 101–112.

Donnadieu, Y., Pierrehumbert, R.T., Fluteau, F., Jacob, R., 2006. Modelling the primary control of paleogeography on Cretaceous climate. Earth Planet. Sci. Lett. 248, 426–437.

Drever, J.L., 1997. The geochemistry of natural waters. Prentice Hall, Upper Saddle River, New Jersey 07458. 436 pp.

Eyring, H., 1935. The activated complex in chemical reactions. J. Chem. Phys. 3, 107–115.

Font, E., Nédélec, A., Trindade, R., Macouin, M., Charrière, A., 2006. Chemostratigraphy of the Neoproterozoic Mirassol d'Oeste cap dolostones (Mato Grosso, Brazil): a model for Marinoan cap dolostone formation. Earth Planet. Sci. Lett. 250, 89–103.

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

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

Goddéris, Y., Joachimski, M.M., 2004. Global change in the late Devonian: modelling the Frasian-Famennian short-term carbon isotope excursions. Palaeogeogr. Palaeoclimatol. Palaeoecol. 202, 309–329.

Goddéris, Y., Nédélec, A., Donnadieu, Y., Dupré, B., Dessert, C., François, L.M., Grard, A., Ramstein, G., 2003. The Sturtian glaciation: Fire and ice. Earth Planet. Sci. Lett. 211, 1–12.

Goddéris, Y., François, L.M., Probst, A., Schott, J., Moncoulon, D., Labat, D., Viville, D., 2006. Modelling weathering processes at the catchment scale with the WITCH numerical model. Geochim. Cosmochim. Acta 70, 1128–1147.

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

Held, I.M., Soden, B.J., 2006. Robust responses of the hydrological cycle to global warming. J. Clim. 19, 5686–5699.

Higgins, J.A., Schrag, D.P., 2003. Aftermath of a snowball Earth. Geochem. Geophys. Geosyst. 4, doi:10.1029/2002GC000403.

Hoffman, P.F., Kaufman, A.J., Halverson, G.P., Schrag, D.P., 1998. A Neoproterozoic Snowball Earth. Science 281, 1342–1346.

Hoffman, P.F., Halverson, G.P., Schrag, D.P., 2004. On the consistent sequence of mineralogy and bedforms in post-glacial cap carbonates of Marinoan (ca 635 Ma) age on different paleomargins: steady-state ocean stratification during a post-snowball Earth marine transgression? Annual Geological Society of America Annual meeting, Denver.

Hoffman, P.F., Halverson, G.P., Domack, E.W., Husson, J.M., Higgins, J.A., Schrag, D.P., 2007. Are basal Ediacaran (635 Ma) post-glacial "cap-dolostones" diachronous? Earth Planet. Sci. Lett. 258, 114–131.

Holland, H.D., 1978. The chemistry of Oceans and Atmosphere. Wiley, New York. 351 pp.

Hyde, W.T., Crowley, T.J., Baum, S.K., Peltier, W.R., 2000. Neoproterozoic "snowball Earth" simulations with a coupled climate/ice sheet model. Nature 405, 425–429.

Jacob, R., 1997. Low frequency variability in a simulated atmosphere ocean system. Univ. Wisconsin-Madison.

Kahle, M., Kieber, M., Reinhold, J., 2002. Carbon storage in loess derived surface soils from Central Germany: influence of mineral phase variables. J. Plant. Nutr. Soil Sci. 165, 141–149.

Kaseman, S.A., Hawkesworth, C.J., Prave, A.R., Fallick, A.E., Pearson, P.N., 2005. Boron and calcium isotope composition in Neoproterozoic carbonate rocks from Namibia: evidence for extreme environmental change. Earth Planet. Sci. Lett. 231, 73–86.

Kilner, B., MacNiocail, C., Brasier, M., 2005. Low-latitude glaciation in the Neoproterozoic of Oman. Geology 33, 413–416.

Le Hir, G., Goddéris, Y., Donnadieu, Y., Ramstein, G., 2008a. A geochemical modelling study of the evolution of the chemical composition of seawater linked to a snowball glaciation. Biogeosciences 5, 253–267.

Le Hir, G., Goddéris, Y., Ramstein, G., Donnadieu, Y., 2008b. A scenario for the evolution of the atmospheric pCO₂ during a snowball Earth. Geology 36, 47–50.

Meert, J.G., Torsvik, T.H., 2003. The making and unmaking of a supercontinent: Rodinia revisited. Tectonophysics 375, 261–288.

Nédélec, A., Affaton, P., France-Lanord, C., Charrière, A., Alvaro, J.J., 2007. Sedimentology and chemostratigraphy of the Bwipe Neoproterozoic cap carbonates (Ghana, Volta Basin): a record of microbial activity in a peritidal environment. C. R. Geosci. 339, 516–518.

Nogueira, A.C.R., Riccomini, C., Sial, A.N., Moura, C.A.V., Trindade, R., Fairchild, T.R., 2007. Carbon and strontium isotope fluctuations and paleoceanographic changes in the late Neoproterozoic Araras carbonate platform, southern Amazon craton, Brazil. Chem. Geol. 237, 168–190.

Oelkers, E.H., Schott, J., Devidal, J.L., 1994. The effect of aluminum, pH and chemical affinity on the rates of aluminosilicate dissolution reactions. Geochim. Cosmochim. Acta 58, 2011–2024.

Oliva, P., Viers, J., Dupré, B., 2003. Chemical weathering in granitic crystalline environments. Chem. Geol. 202, 225–256.
Peyser, C.E., Poulsen, C.J., 2008. Controls on Permo-Carboniferous precipitation over tropical Pangaea: a GCM sensitivity study. Palaeogeogr. Palaeoclimatol. Palaeoecol. 268, 181–192.

Pierrehumbert, R.T., 1999. Subtropical water vapor as a mediator of rapid global climate change. In: Clarks, P.U., Webb, R.S., Keigwin, L.D. (Eds.), Mechanisms of global change at the Millenial time scales. Geophys. Monogr. Ser., vol. 112. American Geophysical Union, Washington D.C., p. 394.

Pierrehumbert, R.T., 2002. The hydrologic cycle in deep-time climate problems. Nature 419, 191–198.

Pierrehumbert, R.T., 2004. High levels of atmospheric carbon dioxyde necessary for the termination of global glaciation. Nature 429, 646–649.

Poulsen, C.J., 2003. Absence of a runaway ice-albedo feedback in the Neoproterozoic. Geology 31, 115–118.

Poulsen, C.J., Jacob, R.L., 2004. Factors that inhibit snowball Earth simulation. Paleooceanography 19, doi:10.1029/2004PA001056.

Rapp, A.D., Kummerow, C., Berg, W., Griffith, B., 2005. An evaluation of the proposed mechanism of the adaptative infrared iris hypothesis using TRMM VIRS and PR measurements. J. Clim. 18, 4185–4194.

Rowley, D.B., 2002. Rate of plate creation and destruction: 180 Ma to present. Geol. Soc. Am. Bull. 114, 927–933.

Schott, J., Oelkers, E.H., 1995. Dissolution and crystallization rates of silicate minerals as a function of chemical affinity. Appl. Chem. 67, 903–910.

Sverdrup, H., Warfinge, P., 1995. Estimating field weathering rates using laboratory kinetics. Rev. Mineral. 31, 485–541.

Trindade, R.I.F., Font, E., D'Agrella-Filho, M.S., Nogueira, A.C.R., Riccomini, C., 2003. Low-latitude and multiple geomagnetic reversals in the Neoproterozoic Puga cap carbonate, Amazon craton. Terra Nova 15, 441–446.

Vecchi, G.A., Soden, B.J., Wittenberg, A.T., Held, I.M., Leetma, A., Harrison, M.J., 2006. Weakening of tropical Pacific atmospheric circulation due to anthropogenic forcing. Nature 441, 73–76.

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

White, A.F., Blum, A.E., 1995. Effects of climate on chemical weathering in watersheds. Geochim. Cosmochim. Acta 59, 1729–1747.
<figure><img src="image_52.png" /></figure>

# The climate change caused by the land plant invasion in the Devonian

Guillaume Le Hir<sup>a,*</sup>, Yannick Donnadieu<sup>b</sup>, Yves Goddéris<sup>c</sup>, Brigitte Meyer-Berthaud<sup>d,e</sup>, Gilles Ramstein<sup>b</sup>, Ronald C. Blakey<sup>f</sup>

## ARTICLE INFO

**Article history:**

Received 2 March 2011

Received in revised form 15 July 2011

Accepted 30 August 2011

Available online 29 September 2011

Editor: P. DeMenocal

**Keywords:**

Devonian

land plants

trees

climate modeling

Palaeozoic

carbon cycle

## ABSTRACT

Land plants invaded continents during the Mid-Paleozoic. Their spreading and diversification have been compared to the Cambrian explosion in terms of intensity and impact on the diversification of life on Earth. Whereas prior studies were focused on the evolution of the root system and its weathering contribution, here we used a coupled climate/carbon/vegetation model to investigate the biophysical impacts of plant colonization on the surface climate through changes in continental albedo, roughness, thermal properties, and potential evaporation. From the Early to the Late Devonian, our model simulates a significant atmospheric CO₂ drop from 6300 to 2100 ppmv that is due to an increase in the consumption of CO₂ though continental silicate weathering. The continental drift and the climatic changes promoted by land plants explain this trend. The simulated CO₂ drawdown is paradoxically associated with unchanged temperatures. We show here that the CO₂ drop is counteracted by a large warming resulting from the surface albedo reduction caused by the appearance of an extended plant-cover. If CO₂ is consensually assumed as the main driver of the Phanerozoic climate, this paper demonstrates that, during land-plant invasion, the modifications of soil properties could have played in the opposite direction of the carbon dioxide fall, hence maintaining warm temperatures during part of the Devonian.

© 2011 Elsevier B.V. All rights reserved.

## 1. Introduction

Earliest bryophyte-like plants evolved in the mid-Ordovician (~470 Ma). Land plants diversified but remained reduced in size following the origin of the first vascular organisms approximately 450 million years ago, then gradually increased in height and complexity during the Early Devonian (~415 Ma), achieved the tree habit in the Middle Devonian (~390 Ma), and evolved large leaves and subsequently seeds, in the Late Devonian (~365 Ma) (Algeo and Scheckler, 1998; Beerling et al., 2001; Kenrick and Crane, 1997; Stein et al., 2007). The first stages of plant colonization were confined to wetlands (Greb et al., 2006). Prior to the Early Devonian, the interplay between plants and their immediate abiotic environment is considered as insignificant due to their small size, shallow root system, and lack of leaves (Driese and Mora, 2001). Land plants then progressively colonized all habitats in lowlands. Seed plants ultimately inhabited drier habitats, possibly uplands, in the late Late Devonian

(Algeo and Scheckler, 1998; Fairon-Demaret and Hartkopf-Fröder, 2004). The rise of large vascular plants with an extended leaf-cover and their expansion into extrabasal settings led to substantial changes in the terrestrial landscapes, and therefore in the global environment.

It is frequently assumed that the increase in plant size and a wider geographic distribution promoted chemical weathering processes (Algeo and Scheckler, 1998; Beerling and Berner, 2005; Berner, 2004). Studies on the fluxes of dissolved cations for water draining forested and bare surfaces have suggested an enhancement of the weathering by factors ranging from 3 to 10 (Berner, 2004; Bormann et al., 1998; Moulton and Berner, 1998). Based on these field results, carbon cycle models estimate that the paleo-pCO₂ should probably be 2 to 5 times higher if vascular plants had never evolved (Berner, 1987, 1991, 2001, 2006). These studies were focused on the acceleration of the continental weathering rate owing to deeply rooted soils. Conversely and to our knowledge, none modeling studies have investigated how changes in the surface properties may have affected the water cycle, the continental weathering rate and hence the climate in response to the land plant invasion during the Paleozoic. Indeed, apart from directly impacting the weathering rates as described above, vascular plants also promote transpiration, reduce the surface

<footnote><sup>a</sup> Institut de Physique du Globe de Paris, Université Paris7-Denis Diderot, 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote><sup>b</sup> LSCE, CNRS-CEA-UVSQ, Gif-sur-Yvette, France</footnote>

<footnote><sup>c</sup> Geosciences-Environnement-Toulouse, CNRS-Université de Toulouse, Observatoire Midi-Pyrénées, Toulouse, France</footnote>

<footnote><sup>d</sup> Université Montpellier 2, UMR AMAP, Montpellier, F-34000, France</footnote>

<footnote><sup>e</sup> CNRS, UMR AMAP, Montpellier, F-34000, France</footnote>

<footnote><sup>f</sup> Northern Arizona University, Flagstaff, USA</footnote>

<footnote>* Corresponding author at: Institut de Physique du Globe de Paris, 1, rue Jussieu, 75238 Paris cedex 05, France. Tel.: + 33 1 83 95 74 97.</footnote>

<footnote>E-mail address: lehir@ipgp.fr (G. Le Hir).</footnote>
albedo, and increase the soil roughness (Algeo et al., 2001; Berner, 2001; Donnadieu et al., 2009). Because continental temperature and runoff directly influence the chemical weathering, quantifying such changes appears essential to better constrain the climatic evolution in relation with the plant continental invasion and the carbon cycle perturbation. To explore the climate-vegetation interplay, we have used a fully coupled 3D carbon-climate-vegetation model (Fig. 1) to simulate the long-term climate evolution during the advent of the land plants. Finally, by using a set of sensitivity experiments, we have determined how each characteristic of the Devonian vegetation has influenced the climate.

## 2. Methods and tools

### 2.1. Climate-carbon model (GEOCLIM)

All experiments conducted in this study are performed with the GEOCLIM model (Donnadieu et al., 2004, 2006; Le Hir et al., 2009). This model includes a General Circulation Model (GCM), here FOAM (for Fast Ocean Atmosphere Model) (Jacob, 1997) where the ocean is represented by a 50 m depth mixed-layer (slab ocean), coupled to a box model of the geological carbon-alkalinity cycles (COMBINE (Godderis and Joachimski, 2004)) (Fig. 1). The spatial resolution of GEOCLIM (7.5° long × 4.5° lat) allows us to capture the evolution of the main features of Devonian geographies and their impacts on the carbon cycle. According to a given set of boundary conditions (paleo-geographic configuration, solar constant, and solid Earth degassing), GEOCLIM calculates the corresponding steady-state atmospheric CO₂ and climate pattern. The pCO₂ calculation relies on interactive processes loop between the FOAM and COMBINE models. Fixing the solid Earth CO₂ degassing to a constant value, the GEOCLIM model estimates the pCO₂ required to balance the degassed carbon flux by the flux of carbon consumed through continental silicate weathering (Walker et al., 1981). The outgassed rate of carbon is kept constant at its supposed present day value: 6.8 × 10¹² mol CO₂ yr⁻¹ required to balance the present day CO₂ consumption by silicate weathering, (Gaillardet et al., 1999), which appears to be a conservative assumption, due to the lack of consensus about the Earth degassing rate for the last 200 million years (see discussion in Donnadieu et al. (2006)). The carbon consumption flux through silicate weathering is calculated by COMBINE using phenomenological laws (Dessert et al., 2003; Oliva et al., 2004). The model accounts for the weathering

of basaltic and granitic outcrops (all grid cells have rock types in modern proportions: 4% of basalt/96% of granite).

$$F_{sw} = k \sum_{j} A_{j} Run_{j} \exp \left[ - \frac{E_{a}}{R} \left( \frac{1}{T_{j}} - \frac{1}{T_{0}} \right) \right] \quad (1)$$

where $T_j$ is the mean annual temperature and $Run_j$ the mean annual runoff of the continental grid box j. $A_j$ is the area of the grid box. $E_a$ is an apparent activation energy set at 48,200 J/mol for granites (Oliva et al., 2004) and 42,300 J/mol for basalts (Dessert et al., 2001). k is a constant calibrated for each lithology, so that total continental silicate weathering consumes a total amount of 6.8 × 10¹² mol/yr of CO₂, 30% being consumed by basalt weathering (Dessert et al., 2003).

The climatic variables $T_j$ and $Run_j$ are computed by FOAM forced by a large range of pCO₂ (30, 25, 20, 15, 10, 8, 6 and 4PAL, one PAL (Preindustrial Atmospheric Level) being equal to 280 ppmv). Climate reaches equilibrium with the prescribed boundary conditions after 30 simulated years; the analyses are performed on the average of the next ten simulated years. Starting from an arbitrary atmospheric CO₂ value, the GEOCLIM model calculates the temperature and runoff of each grid cell through a linear interpolation of the temperature and runoff predicted by each climatic simulation and for each continental grid cell. A new CO₂ consumption through silicate weathering is then estimated for each grid box with Eq. (1). This flux is summed up over the whole continental surface and total silicate weathering is estimated, and consequently a new atmospheric CO₂. This procedure is repeated until a steady-state is reached that corresponds to a stable CO₂ and climate. For more details, an exhaustive GEOCLIM description can be found in Donnadieu et al. (2006).

### 2.2. Vegetation model description (SLAVE-Jr) and coupling procedure

Similar to Donnadieu et al. (2009), we improve the GEOCLIM model by incorporating a dynamic model for the vegetation to account for the vegetation impact on the climate. Indeed, the influence of vegetation on climate occurs through the fluxes of energy, water vapor and momentum between land and atmosphere. Vegetation modifies the albedo of the continental surfaces and absorbs part of the heat which modifies the sensible heat flux. Plants also transpire and change the latent heat flux. Furthermore, they hinder air movements which transfer momentum to the atmosphere. The calculated vegetation distribution directly impacts on the following physical parameters in FOAM: (1) the albedo of each continental grid element, (2) its thermal properties, (3) its roughness (height of plants governing the aerodynamic resistance of wind speed), (4) and its potential evaporation (incorporating stomata effect).

Consequently, to explore vegetation/climate interplays, the main improvement of GEOCLIM is the coupling with the vegetation model SLAVE-Junior (hereafter SLAVE-Jr) so that vegetation distribution can be coevally calculated with the climatic reconstruction. SLAVE-Jr is a slightly modified version of SLAVE (Freidlingstein, 1995; Friedlingstein et al., 1992) adapted for deep-time purpose. It originally represented 9 vegetation types, but for the Devonian, only three vegetation types (basal embryophyte communities, cladoxylopid forests comprising tree fern-like trees, and Archeopteris forests consisting of woody trees) are considered as appropriate representation of main steps in the land-plant evolution and given our limited knowledge of the biosphere characteristics in the Devonian.

The spatial distribution of each vegetation type is based on two climatic variables: (1) the annual mean continental temperature, and (2) water availability for plants (or degree of saturation of the soil) for each grid point. The degree of saturation (hereafter DOS) is defined as the water volume (soil moisture variable) divided by the maximum water content in the soil (15 cm, bucket model

<figure><img src="image_53.png" /><figcaption>Fig. 1. Models architecture. The GEOCLIM model comes from the coupling of the FOAM model and of the COMBINE model. DOS is degree of saturation of soils, Air temp is the annual mean continental air temperature at the surface level, Plant Functional Types (PFT) are fully described in Section 2.3 and Table 1.</figcaption></figure>
assumption). This water content is directly governed by the water availability to the pore space and ranges from 0 (dry) to 1 (saturated). In contrast with modern plants, the shallow root systems of the first land plants were probably limited by the available soil water. This parameterization tends to better account for this effect. Concerning the temperature sensitivity, we consider that a grid box with an annual continental temperature below 0 °C is a factor prohibiting the growth of the tree ferns that are not adapted to cold seasons. Explicit conditions to install vegetation types for Plant Functional Types (hereafter PFT1, PFT2 and PFT3) are given in Table 1b.

To couple FOAM with SLAVE-Jr, we have written a *python* program (*python* is a dynamic object-oriented programming language) that changes the climatic parameters in the FOAM boundary files according to the vegetation distribution given by SLAVE-Jr. Starting from the climate and atmospheric information as inputs, SLAVE-Jr dynamically computes spatially explicit transient vegetation composition in terms of plant functional groups, and their associated water budgets. In the coupling procedure, FOAM is run for five years, the annual terrestrial temperature and precipitations of the last year are saved to force SLAVE-Jr which calculates the PFT pattern. This process is repeated throughout the first 20 model years to take into account the vegetation-climate feedbacks, a time to ensure that the vegetation reaches the equilibrium coevally with the climatic conditions. The climates calculated by the coupled FOAM-SLAVE-Jr model are then used in COMBINE to estimate past atmospheric CO₂ and steady state climate following the procedure described above. We eventually reach a climate at steady-state associated with a biome distribution for each simulated geological period in equilibrium with the calculated pCO₂.

## 2.3. Devonian land plants

Previous conclusions on early land plant evolution distinguished three successive phases in the terrestrialization of plants (Kenrick and Crane, 1997). (1) The “Eoembryophytic phase” when most land plants are at a bryophytic grade of organization is followed by (2) an “Eotracheophytic phase” corresponding to the rise of the earliest vascular plants in the Silurian and Early Devonian. (3) The “Eutracheophytic phase”, which is the main focus of this paper, initiates in the Early Devonian and corresponds to an extensive diversification of the morphologies and reproductive strategies of land plants. During this phase, the evolution of trees that occurred separately in different clades represents a major breakthrough with a putative

impact on climate changes (Algeo et al., 2001; Beerling and Berner, 2005). Interestingly, subdividing the “Eotracheophytic phase” based on the presence versus the absence of trees and on the successive dominant types of trees of the Devonian results in recognizing the Early, Middle and Late Devonian respectively as three significant phases in the evolution of land plants (see below).

An accurate appraisal of the composition, structure, ecological requirements and geographical distribution of the palaeovegetations for each phase is currently difficult. Our assessments are based on extrapolations from the limited number of palaeoecological studies that have been realized so far, most from Euramerican localities (Cressler et al., 2010; Driese et al., 1997; Hotton et al., 2001; Mintz et al., 2010; Scheckler, 1986; Steemans et al., 2010).

Finally, assessing the most appropriate extant analogs for the Devonian vegetations in order to parameterize the land surface features of past landscapes is also a challenge. Our approach consists of adopting the same criteria to define the main vegetation types of the Devonian as those used by Matthews (1984) to classify the extant types of vegetation. These criteria include the height, presumed density and main physiological traits of the dominant lifeforms. Because we use ecological rather than systematic criteria, and despite the fact that they comprise no closely related plants, we assume that the extant types of vegetation we have selected are the closest available analogs for the Devonian forms (hereafter called PFTs for Plant Functional Types). We then derive the land-surface properties for the Devonian PFTs from the values of the corresponding extant vegetation types provided by Matthews (1984) and Hack et al. (1993). Considering these assumptions, the proposed vegetation types and extant analogs used to parameterize land-surface boundary conditions for the Devonian (up to the mid-Famennian) are the following. They are summarized in Table 1.

(1) Vegetated Early and Middle Devonian landscapes were dominated by PFT1 which comprises basal land plants. These include the eotracheophytic rhyniophytes and rhyniophytoids, and the eutracheophytic zosterophylls, early lycopsids and basal euphylophytes of low to moderate stature. These plants are shallowly rooted and characterized by a homosporous reproduction necessitating free water for fertilization. Interpreted as obligate dwellers of wet substrates (Greb et al., 2006), some may have withstood seasonal freezing or aridity (Wellman and Gray, 2000; Wellman et al., 2000). Many formed

<table><caption>**Table 1**
Devonian vegetation (upper part) and boundary conditions (lower part).</caption><thead><tr><th>Simulations</th><th>Periods</th><th>Main plant type</th><th>Height</th><th>Density</th><th>Phenology</th><th>Morphology of photosynthetic structures</th><th>Closest extant vegetation analog (Matthews, 1984)</th></tr></thead><tbody><tr><td>LP1</td><td>Early Devonian</td><td>PFT1, basal embryophytes</td><td>0.5–1 m up to 2 m</td><td>–</td><td>–</td><td>Narrow, dissected</td><td>Grassland/forb formation</td></tr><tr><td>LP2</td><td>Mid-Devonian (or Late Early Devonian)</td><td>PFT1, basal embryophytes<br>PFT2, Tree-ferns (Cladoxylopsid)</td><td>0.5–1 m up to 2 m potentially >5 m<br>Potentially >5 m</td><td>–</td><td>–</td><td>Narrow, dissected</td><td>Grassland/forb formation</td></tr><tr><td>LP3</td><td>Late Devonian (mid Frasnian to mid-Famennian)</td><td>PFT3, Woody trees (Archeopteris)</td><td>Potentially >5 m</td><td>Open stands<br>Closed forest</td><td>Evergreen</td><td>Narrow, dissected</td><td>Evergreen needleleaved</td></tr><tr><td>Surface properties (PFT)</td><td>Albedo (%) (total/vis/near IR)</td><td>Roughness (m)</td><td>Evapotranspiration coefficient</td><td></td><td>Thermal property</td><td>Broad-leaved</td><td>Rainforest</td></tr><tr><td>Rocky desert</td><td>32/24/40</td><td>0.05</td><td>0.1</td><td>7</td><td></td><td></td><td></td></tr><tr><td>Basal embryophytes (PFT1)</td><td>20.5/9/32</td><td>0.1</td><td>0.1</td><td>6</td><td></td><td></td><td>0.5&lt;DOS&lt;0.8</td></tr><tr><td>Tree-ferns (PFT2)</td><td>19.5/14/25</td><td>1</td><td>0.25</td><td>5</td><td></td><td></td><td>0.8&lt;DOS and annual temp &gt; 0 °C</td></tr><tr><td>Woody trees (PFT3)</td><td>12/6/18</td><td>2</td><td>0.25</td><td>4</td><td></td><td></td><td>0.3&lt;DOS</td></tr></tbody></table>

To define Devonian plant properties, we used values of their respective analog extracted from the Matthews classification (1984) (see Section 2.3 for an exhaustive discussion concerning the choice of modern analogs). Then, land surface boundary conditions are specified in FOAM using vegetation-related surface features. The data base used to prescribe the plant albedo is the NASA Tech. Memo (Matthews, 1984). For the roughness, evapo-transpiration coefficient, and thermal property, values come from the NCAR Tech. Note (Hack et al., 1993). The surface extension conditions are determined as a function of assumptions based on paleo-botanical features (see Section 2.3).
monotaxic patches of clonal plants forming dense stands that excluded other species (Hotton et al., 2001). The closest analogs found in Matthews (1984) are “tundra formations” (referred to as # 16 in Matthews' Table 1) and “grasslands” (referred to as # 19 in Matthews' Table 1). The former consists of dwarf-scrub communities, the latter to either actual grasslands or to forb formations, i.e. formations of non graminoid herbaceous plants such as ferns that include plants exceeding 50 cm in height. “Tundra formations” are more appropriate analogs for the eotracheophytic vegetations of late Silurian and earliest Devonian age. However, in order to take into account the increase in size of the plants accompanying the diversification of the early eutracheophytes prior to the advent of the tree habit, we have selected the “grasslands/forb formations” as the closest analogs for the Early and Middle Devonian vegetation types devoid of trees. This vegetation type better corresponds to the communities reported from the Emsian Battery Point Formation from Gaspé (Québec) for example (Hotton et al., 2001). Selecting “tundra formations” as the closest analog for PFT1 would have probably underestimated the impact of this PFT for both the late Early Devonian and the Middle Devonian. Our choice of “grasslands/forb formations” probably overestimates the impact of this PFT for the Lochkovian and Pragian stages of the Devonian. The ecological requirements mentioned above for PFT1 supposes a small range of acceptance concerning the DOS (from 0.5 to 0.8) and absence of limitation for the temperature.

(2) The Middle Devonian saw the rise of PFT2 which corresponds to the earliest tree communities. It comprised tree fern-like cladoxylopsids that could exceed 5 m in height (Stein et al., 2007). By analogy with extant trees, the limited extent of their root system and narrowness of their crown lacking large laminate leaves suggest that these cladoxylopsids may have lived in open stands. Photosynthesis was performed by distal branches that were green and may have superficially resembled those of extant Cupressaceae. Stumps rooted in paleosols indicate that they occupied wetlands (Driese et al., 1997; Mintz et al., 2010). For these reasons, we choose Matthews (1984) “evergreen needle-leaved woodlands” (referred to as # 10 in Matthews' Table 1), a composite vegetation type comprised of 50% of evergreen woodlands, 25% of “grasslands” and 25% of barren surface, as the closest analog for the cladoxylopsid forests with the requirements for an elevated DOS (here fixed over 0.8) of their habitat. Because of the large amount of living cells characterizing their trunk, cladoxylopsid trees probably avoided freezing climates (Meyer-Berthaud et al., 2010) which exclude lands experiencing temperatures below zero and long cold seasons. Indeed, no indication of seasonality such as growth rings or differences in the size of the branch scars along the trunk has been reported for these trees, which suggests their preference for equable environmental conditions.

(3) In the Late Devonian spread a second type of forests, PFT3, comprised of the progymnosperm genus *Archeopteris*. *Archeopteris* had a worldwide distribution that ranged from low to high paleolatitudes (Anderson et al., 1995; Cressler et al., 2010; Retallack, 1997). These trees, which may have reached 40 m in height, possessed large leaves and an extended root system. Their trunks, branches, and roots were essentially made of wood. They formed dense shady forests, a strategy that overcame the necessity of moisture linked to their heterosporous reproduction (Beck, 1964; Driese et al., 1997). Because of these traits — depth of the root system, wood instead of living cells as the essential component of the biomass, optimal environment for reproduction independent from external conditions and provided by the forest itself — we assume that the Late Devonian PFT3, could survive with a DOS below 30%. In

addition, the structural resemblance of *Archeopteris* with that of extant conifers suggests that these trees may have supported freezing conditions and suppose an absence of thermal restriction for their colonization. They inhabited diverse types of habitats from poorly to well-drained lowlands (Driese et al., 1997; Scheckler, 2001), and showed a wide range of ecological strategies including opportunistic growth (Mintz et al., 2010). They had the ability to stop their growth, produce growth rings (Arnold, 1930; Lemoigne et al., 1983), and generate epicormic shoots (Meyer-Berthaud et al., 2000), a set of traits that may have favored their settlements in habitats subjected to environmental stresses or even disturbances, and under a wider range of climatic conditions than any of the other trees of the Late Devonian (Cressler, 2006; Cressler et al., 2010; Fairon-Demaret et al., 2001). *Archeopteris* shed leafy branches instead of leaves. But the deciduous/persistent character of the leafy branches is equivocal (Beck, 1971). Many woody remains of *Archeopteris* from equatorial and paleotropical deposits show either subtle (Lemoigne et al., 1983; Meyer-Berthaud et al., 1997; Snigirevskaya and Snigirevsky, 2001) or no growth rings (Arnold, 1930; Galtier et al., 1996; Meyer-Berthaud et al., 2004), suggesting that leafy branches remained attached to the trees for several years (Falcon-Lang, 2000). Added to the height and density of the stands, the persistency of the leafy units resulted in forests that were presumably closed and evergreen which legitimizes our choice for “rainforests” as the closest extant analogs for most of the PFT3 represented in the fossil record. Dense broad-leaved forests thrived at high paleolatitudes (85°S) during the Permian (Taylor et al., 1992) and shed their leaves during winter, a strategy possibly adopted by their Late Devonian counterparts. However, although archaeopterid spores have been recorded in the Givetian (385–392 Ma) of Antarctica, and compressions of leafy branches in the Famennian (374–359 Ma) of South Africa, no woody remains were yielded from high paleolatitude outcrops of Devonian age, and the phenology of archaeopterid forests in these zones is unknown. Treating all Late Devonian forests as a single PFT analogous to extant rainforests is another limitation of this study which may overestimate the impact of the actual vegetation at this time.

Evidence for the occurrence of early seed plants before the late Famennian is poor. Published analyses of the distribution of plants in the late Famennian localities of the Catskill Delta show that early seed plants were ruderal in disturbed habitats within lowlands (Cressler et al., 2010; Scheckler, 1986). This indicates the pivotal role of such plants in ecological successions but does not infer that these plants covered large surfaces in the Devonian landscapes even if it is suspected that they have been more significant than suggested by their record at this time (Fairon-Demaret and Hartkopf-Fröder, 2004). Because of these uncertainties about the actual distribution of this plant group in Late Devonian landscapes, especially in uplands, and the spatial resolution of the model, we choose to follow a conservative approach and did not take into account the seed plants in this work.

## 2.4. Boundary conditions and simulation design

The only three forcing functions are the solid Earth degassing rate (here fixed at its present-day value due to the lack of constraints), the continental configuration, and the slowly increasing solar constant. In the present study, we use the palaeogeographies from Ron Blakey (2003), where the Devonian is mainly characterized by a southward migration of Gondwanaland while Laurentia, Baltica and Siberia craton remained located at low latitudes. To represent this continental drift, two maps have been used for the Devonian. The first one is
assumed to be representative of the 400–370 Ma interval (Early and Middle Devonian), and the second one stands for the 370–360 Ma interval (Late Devonian). Concerning the solar constant, we have supposed that the solar luminosity increases from 1323 W/m² (3.3% reduction relative to present day value) in the Early Devonian to 1326 W/m² in the Late Devonian, according to stellar evolution models (Gough, 1981).

To model the land plant invasion and investigate climate–carbon-vegetation interplays, five GEOCLIM runs are presented in this contribution. We first quantify the effect of the continental drifting between the Early and the Late Devonian without plants (no land-plant simulations, Table 2). For these runs, we prescribe a continental surface covered by a rocky regolith (Tables 1 and 2). Then, the land-plant invasion is modeled by three sets of simulations (land-plant runs or LP, see Tables 1–2). LP1 simulations represent the Earliest Devonian case characterized by the existence of primitive plants (PFT1). LP2 runs simulate the Middle Devonian (or Late Early Devonian) where PFT1 and PFT2 coexist. Both LP1 and LP2 are performed assuming an Early Devonian continental configuration. Finally LP3 simulation describes the Late Devonian vegetation dominated by woody trees (PFT3). The Late Devonian continental configuration is used for LP3. For all the LP (for land plant) simulations, GEOCLIM is coupled to the terrestrial biosphere model SLAVE-Jr. In addition, sensitivity runs have been performed to test the robustness of our results according to the boundary conditions that we have initially fixed (see Sections 3.4 and 3.5).

### 3. Surface climate and vegetation interplays during the land-plant invasion

The impact of plant colonization on the surface climate is investigated using a series of general circulation model simulations, coupled to a carbon and vegetation models. Since the modeling outputs are vegetation maps, CO₂ and climate, those three aspects will be successively described.

#### 3.1. Calculated vegetation maps during Devonian

In vegetation maps calculated by the LP1 run, basal embryophytes (Fig. 2A) are mainly restricted to the coastal environments of north Gondwana, Laurussia and north Siberia. This limited expansion is the result of (1) an arid climate persisting during the Early Devonian due to the continental configuration, and (2) the water requirement

<figure><img src="image_54.png" /><figcaption>**Fig. 2.** Devonian potential vegetation distribution. Maps show the maximal potential distribution pattern for three timeslices, in equilibrium with the calculated pCO₂ (respectively 5144, 4843 and 2125 ppmv for LP1, LP2 and LP3 runs).</figcaption></figure>

<table><caption>**Table 2**
Main modeling results.</caption><thead><tr><th rowspan="2">Simulations</th><th colspan="3">Early Devonian</th><th colspan="2">Late Devonian</th></tr><tr><th>No land plants</th><th>LP1</th><th>LP2</th><th>No land plants</th><th>LP3</th></tr></thead><tbody><tr><td>Surface (10<sup>6</sup> km<sup>2</sup>)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rocky regolith</td><td>122</td><td>95</td><td>72</td><td>104</td><td>27</td></tr><tr><td>Primitive plants</td><td>-</td><td>27</td><td>30</td><td>-</td><td>-</td></tr><tr><td>Giant ferns</td><td>-</td><td>-</td><td>20</td><td>-</td><td>-</td></tr><tr><td>Tree-like</td><td>-</td><td>-</td><td>-</td><td>-</td><td>77</td></tr><tr><td>Cont albedo (%)</td><td>32.0</td><td>29.9</td><td>27.0</td><td>32.0</td><td>17.2</td></tr><tr><td>Planetary albedo (%)</td><td>17.8</td><td>17.5</td><td>17.2</td><td>17.4</td><td>15.8</td></tr><tr><td>pCO<sub>2</sub> (ppmv)</td><td>6296</td><td>5144</td><td>4843</td><td>3954</td><td>2125</td></tr><tr><td>Tcont (°C)</td><td>16.2</td><td>15.5</td><td>15.5</td><td>16.9</td><td>16.5</td></tr><tr><td>SST (°C)</td><td>26.5</td><td>25.9</td><td>26.0</td><td>24.6</td><td>23.7</td></tr><tr><td>T<sub>global</sub> (°C)</td><td>24.1</td><td>23.4</td><td>23.5</td><td>23.1</td><td>22.2</td></tr><tr><td>Runoff (cm/yr)</td><td>36.6</td><td>36.4</td><td>36.3</td><td>41.3</td><td>40.0</td></tr><tr><td>Cont. precip (cm/yr)</td><td>103.1</td><td>103.8</td><td>107.6</td><td>112.3</td><td>129.5</td></tr><tr><td>Cont. evap (cm/yr)</td><td>66.4</td><td>67.3</td><td>71.3</td><td>70.9</td><td>89.4</td></tr><tr><td>Cont. latent heat flux (W/m<sup>2</sup>)</td><td>53.8</td><td>54.6</td><td>57.8</td><td>57.5</td><td>72.6</td></tr></tbody></table>

The continental albedo does not include the snow-cover. Planetary albedo computed in clear sky.

of basal embryophytes (Table 1b). Looking at LP2 simulations, the distribution of basal embryophytes remains quasi-unchanged, but the new PFT2, consisting of “tree ferns”, spreads into low latitudes (<15° latitudes) settling in the wet environments (P-E<sub>low latitudes</sub> > 100 cm/yr) (Fig. 2B). Mid-latitudes (40 to 50°) also take advantage of soil remaining humid throughout the year due to the sea-side proximity (P-E<sub>mid latitudes</sub> > 40 cm/yr). Between 15° and 30° of latitudes, the subsiding branch of the Hadley cell promotes the development of subtropical arid deserts and explains why these regions remain free of vegetation (P-E < 20 cm/yr). Regarding LP3 simulations, 75% of continents are potentially colonized by land plants (PFT3), the remaining 25% being too dry (deserts) to be invaded (but this ratio is slightly depending on topography assumptions). This expansion is the result of less restrictive living conditions (Table 1b). A second cause is related to the tectonic setting. The land distribution, coupled to the sea-level rise, reduces the continental aridity. These two features facilitate the spreading of plants during the Late Devonian, even at high latitudes (Fig. 2C). It should be noted that Fig. 2 displays
maps of potential vegetation, without representing the plant abundance, or the associated productivity.

## 3.2. Devonian simulations predict an important pCO₂ reduction

Regarding runs performed without vegetation (No-LP runs), the pCO₂ drops from 6300 to 3950 ppmv throughout the Devonian (Table 2). This long-term reduction is mainly due to modifications of the global runoff simulated between the Upper and Lower Devonian, respectively 36.6 and 41.3 cm/yr (Table 2). In the Lower Devonian case, in order to consume the 6.8·10¹² mol of carbon released by volcanoes, the CO₂ has to rise (6296 ppmv) until the Earth surface is warm enough ($T_{global} = 24.1$ °C) to counteract the inhibition of weathering by low runoff and balance the solid Earth prescribed degassing. Conversely, the Upper Devonian configuration predicts a wetter climate (+11.4% of runoff, and +8.2% for precipitations, in annual mean), and thus requires a lower pCO₂ (3954 ppmv) to reach the same continental silicate weathering rate.

Fig. 3a shows the silicate weathering fluxes for both continental configurations at the same atmospheric CO₂ concentration (4200 ppmv). Two main differences arise, the first one over the tropical latitudes and the second one over the northern mid-latitudes. They are due to the general northward drift of continents from the Lower Devonian to the Upper Devonian that brings the Laurussia plate into the intertropical convergence zone and the Siberia plate from the tropical divergence zone to the wetter middle latitudes of the Northern hemisphere. These patterns of the plates drifting result in an increase of the silicate weathering through an increase of the continental runoff and promote a more efficient CO₂ consumption.

Fig. 4 shows the simulated atmospheric CO₂ concentrations during the land-plant colonization. The first step of land-plant appearance (i.e. the basal embryophytes, PFT 1) induces a decrease of pCO₂ of 1100 ppmv (from 6296 to 5144 ppmv) compared to the run performed without land plants. An additional 300 ppmv CO₂ decrease is simulated when including the arborescent cladoxylopsids vegetation (PFT 2). Finally, during the Late Devonian, the emergence of the PFT3 associated to the continental drift causes the largest effect on atmospheric CO₂ level with a decrease of more than 2500 ppmv (from 4843 to 2125 ppmv, where 1800 ppmv of reduction is the contribution accountable to land plants, see Fig. 4 and Table 2). These successive steps, in atmospheric CO₂ levels, are the results of an increase in the CO₂ consumption over low latitudes for all simulations and over the southern mid-latitudes for the LP3 runs. Regarding zonal weathering, runoff and temperature profiles at a fixed pCO₂ (Fig. 3b–d), we show that the three steps of the colonization of the land plants have induced a general increase of the continental temperature at all latitudes and a more complex modification of the runoff pattern with positive and negative anomalies. Nevertheless, when comparing runoff modifications and the changes in the weathering fluxes, one can easily attribute most of the increase in the CO₂ consumption to the runoff modifications promoted by the rise of land plants (Fig. 3d) and by the continental drift.

## 3.3. Climate evolution: a temperature decoupled from the pCO₂

While our model simulates a huge CO₂ decrease, the world-wide surface temperature remains rather stable during the Devonian, respectively 23.4, 23.5 and 22.2 °C for LP1, LP2 and LP3 simulations (Table 2). For instance, despite a significant drop in CO₂ (4843 to 2125 ppmv) from LP2 (Middle Devonian) to LP3 (Late Devonian), the temperature remains paradoxically quasi-unchanged (from 23.5 to 22.2 °C), whereas runoff is significantly boosted (36.3 to 40 cm/yr).

To explain this unusual behavior linked to the land-plant expansion, we compare the Earth energy budget evolution from 400 Ma (No-LP runs) to 370 Ma (LP3) to estimate the impact of the vegetated surface rise. At the global scale and in clear sky diagnostic, the net

solar flux absorbed at the TOA (Top of Atmosphere) increases by 6.4 W/m² over the course of the Devonian. This result is explained by the significant reduction of the planetary Earth albedo (from 17.8 to 15.8%) caused by the continental albedo decrease (32.0 to 17.2%, see Table 2). However this albedo reduction does not provide a net warming. Indeed, as our model includes carbon-climate feedbacks, the transient warming caused by the albedo reduction is counteracted by the increase of the weathering rate which drives a pCO₂ drawdown. In consequence, the continental heating provided by the land-plant spreading is balanced by the cooling caused by the CO₂ drawdown. According to the relationship linking the pCO₂ and its radiative forcing, the CO₂ dropping from 6296 to 2125 ppmv creates a negative forcing estimated to –8.9 W/m² (computation based on the standard radiative model of Kiehl and Dickinson (1987)). Hence our model predicts a pCO₂ decrease over the course of the Devonian paradoxically associated to a mean global temperature remaining constant.

<figure><img src="image_55.png" /><figcaption>Fig. 3. Simulated silicate weathering rates, runoff and surface temperatures anomalies at a prescribed pCO₂ (15PAL or 4200 ppmv) for various timeslices and run conditions.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>**Fig. 4.** Simulated pCO₂ obtained with GEOCLIM/SLAVE-JR (green curve) through the Devonian. The blue curve represents pCO₂ estimations based on GEOCARBIII (modified from Van Geldern et al., 2006). Light brown diamonds represent pCO₂ obtained with GEOCLIM without land plants (no-LP runs). The ice-age threshold has been estimated close to 4PAL (1120 ppmv) using the Late Devonian conditions.</figcaption></figure>

### 3.4. Impact of the land plants on the water cycle and weathering rates

In Section 3.2, we have demonstrated that the increase of the silicate weathering was due to an increase in the runoff. Coevally with the colonization of the continental surfaces by land plants, three physical parameters of the climatic model are impacted: the potential evaporation, the roughness and the albedo of the surface. Each of these parameters may modify the continental runoff and hence weathering rates. Regarding Early and Mid-Devonian simulations (LP1, LP2), the enhancement of the runoff (Fig. 3d) and associated weathering rates remain moderate (Fig. 3c). Indeed, (1) basal embryophytes have surface characteristics similar to the rocky regolith, and (2) the arborescent cladoxylopsids vegetation remains too confined (20·10⁶ km²) to significantly affect the global climate via an albedo reduction or via the evapotranspiration process (Table 2). The largest effect occurs during the Late Devonian where runoff and temperature patterns show large positive anomalies (Fig. 3c and d). Another diagnostic obtained with a roughly constant continental temperature but a changing pCO₂ (Table 2) shows that the runoff reaches 40 cm/yr in LP3, compared to 36 cm/yr in LP1 and LP2 runs.

To identify the mechanisms that modify the runoff and consequently the weathering, a supplementary set of sensitivity simulations have been performed, wherein one element of the reference run (albedo, roughness or potential evaporation coefficient) is replaced by the corresponding value typical of a rocky regolith. To inhibit the carbon-climate feedback, the pCO₂ is fixed at 4200 ppmv (15PAL). Regarding the anomalies obtained for weathering rates (Fig. 5), a large difference occurs when the vegetation albedo (0.12) is replaced by its desert value (0.32) (Fig. 5, black curve). This result suggests that a higher continental albedo decreases the runoff, and hence the weathering efficiency. A second order factor is the roughness affecting the latent and sensible heat fluxes. A surprising result is the unchanged pattern when the potential evaporation factor is set to its desert value, which supposes a minor effect caused by the evapo-transpiration process. Consequently the runoff, which drives weathering changes, is mainly controlled by the decrease in surface albedo promoted by the spreading of dense forests.

### 3.5. Importance of the plant albedo

The results of this study suggest that the land plant albedo largely drives an enhancement of the silicate weathering and a decrease in CO₂. As the plant albedo also governs the global hydrological cycle at the first order, it must be considered with care in order to account

<figure><img src="image_57.png" /><figcaption>**Fig. 5.** Sensitivity test to plants surface properties. In each simulation one element of the reference run is replaced by the corresponding desert value (albedo, roughness and evapotranspiration coefficient). To avoid climate-carbon feedbacks and isolate the climate-vegetation interplay, the pCO₂ is fixed at 15PAL.</figcaption></figure>

for the long term response of land plants. To investigate the impact of our parameterization concerning the plant albedo, two alternative albedo scenarios have been tested. In the first one, Late Devonian plants have the tree-fern albedo (19.5%) instead of the woody trees forests (12%). In the second case, we have replaced the initial PFT3 albedo by the basal embryophytes value (20.5%). In all cases the plant pattern distribution and other surface properties are kept constant, and correspond to LP3 run values.

As expected, our results indicate that tree-fern and basal embryophyte albedo lead to the calculation of higher pCO₂, respectively 3161 and 2979 ppmv, in the Late Devonian (compared to 2125 ppmv in the original LP3 run). We have calculated the climate sensitivity to Land Cover Change (LCC, see method in Appendix A), which corresponds to the warming caused by the land albedo radiative forcing. According to our model, when the LCC produces a forcing of 1 W/m² on the continent, the continental temperature increases by ~0.2 °C. The same calculation as a function of the continental albedo suggests a ratio of 1 °C for each 0.036 albedo variation. Hence, assuming that the pre-land plants' surface is correctly represented by a rocky desert, the progressive change of the continental albedo to the Late Devonian, from 32.0 to 17.2% (Table 2) indicates that land plants have induced a global warming of at least 4 °C. It should be noted that this warming is not observed in simulations due to the parallel reduction of the pCO₂ (see Section 3.2).

## 4. Discussion

In such deep time periods, it remains difficult to compare our predicted plant distributions with a coarse resolution to the fossil record (frequently punctual in time and space and allochthonous). However several similarities could be distinguished. For Mid-Devonian (LP2), the modeled phytogeography matches reasonably the fossils outcrops based on spore and megafossil assemblages (Edwards and Wellman, 2001; Scotese et al., 1999), except for South Laurussia. Indeed, according to our model, this region located between 15 and 30°S latitude, is too arid to bear any primitive land-plant implantation. This pattern makes sense because, beyond 20° of latitude, the subsiding branch of the Hadley cell leads to the development of subtropical arid zones. Nevertheless, we also emphasize that the use of the ocean slab model commonly induces slightly extended subsiding regions. Hence, the equator belt, with a very humid and warm climate, tends to be too narrow.
The Middle and Upper Devonian are considered as periods characterized by a weak latitudinal temperature gradient due to generally warm climatic conditions (Joachimski et al., 2009; Van Geldern et al., 2006; Veizer et al., 1999). Indeed miospores typical of mid-to-high latitude environment disappeared in the Mid-Devonian and tropical forms spread into subpolar regions during the Late Devonian. Similarly the ⁸⁷Sr/⁸⁶Sr records in carbonates tend to show an intensified weathering during the Late Devonian (Van Geldern et al., 2006), although the interpretation of the ⁸⁷Sr/⁸⁶Sr as indicative of weathering intensity is highly disputable (Godderis and Francois, 1995; Godderis et al., 2008; Nardin et al., 2010). Finally the oxygen isotopic composition of the surface seawater recorded in sedimentary calcite suggests a global warming during the Late Devonian. However, a global warming episode in the Middle–Late Devonian would be anti-correlated with the CO₂ evolution simulated by numerical models such as GEOCARBII or GEOCARBSULF which predicts a drastic decrease from approximately 3400 to 1100 ppmv (Van Geldern et al., 2006). To reconcile proxies recording a global warming and the modeled evolution of the global carbon cycle during the Devonian, a likely solution would be that the negative climatic forcing driven by the pCO₂ reduction would be transiently overcome by a positive forcing linked to decrease continental albedo associated with a rapid land-plant expansion during the Middle to Late Devonian (Joachimski et al., 2009; Van Geldern et al., 2006). The present study yields tropical sea surface temperatures of 28 °C, 24.7 °C and 26.5 °C for the No-LP, LP2 and LP3 simulations, respectively. Using the Joachimski et al. temperature trend, we get tropical SSTs estimates of 28–32 °C at the Silurian–Devonian boundary (=No-LP), 21–25 °C in the Middle Devonian (=LP2), and 28–32 °C in the Late Devonian (=LP3). Trend in temperature from the Early to the Middle Devonian is in fairly good agreement between our study and the one of Joachimski et al. (2009), although our cooling trend is a bit weaker than the one estimated with data. However, the Late Devonian estimates from data are clearly warmer than our estimates. The model-data disagreement suggests either that (1) pCO₂ did not decline as much during this interval, or (2) the positive radiative forcing associated with Earth-surface albedo changes was larger than that calculated in the present-study. For example, we estimate that the land plant albedo changes associated to their massive extension could provide an important warming (+4.6 °C) if the pCO₂ is not influenced by the inorganic carbon feedback. In addition, moderate differences in term of CO₂ may come from other processes not included in our carbon cycle model (roots, organic burial).

In addition, our results have to be interpreted with caution. During the Latest Devonian, it seems that other factors may have supplanted the effect of the tectonic and of the land plant invasion. Indeed, a cold event recorded by the ⁶¹O in brachiopod and conodont tests is found at the Famennian (365 Ma) (Algeo and Scheckler, 1998; Algeo et al., 2001). Regarding the timing of major changes occurring during the latest Devonian and in agreement with the “Devonian Plant Hypothesis” (Algeo et al., 2001), it could be suggested that the spread of seed plants into upland areas may have been the factor that reduced pCO₂ to levels triggering a continental glaciation event (Fig. 4).<fnref n="1" />

Previous modeling studies suggest that during the Devonian, atmospheric CO₂ was reduced by a factor of 2 (Berner, 2004; Berner and Kothavala, 2001) due to vast quantities of carbon previously located in the atmosphere and now stored in the sedimentary layer. The most accepted mechanism to explain this downward trend was the enhancement of chemical weathering during land plant colonization, through root respiration, mechanical effect of roots and organic

acid release (Berner, 1994; Berner and Kothavala, 2001; Ekart et al., 1999; Mora et al., 1996; Yapp and Poths, 1992). This study shed light on other processes, notably the changes in surface properties, which may increase the silicate weathering enough to decrease atmospheric CO₂ by a factor of 1.55 ± 0.3 (see footnote) according to the vegetation albedo value assumed (from 0.12 to 0.2). Therefore this decrease in CO₂ is of the same order of magnitude as prior estimations which are based on other processes. Hence surface property changes are one of the major actors involved in the CO₂ reduction during Mid-Paleozoic.

The processes we pointed out in this paper are robust, nevertheless their impact on precipitation, silicate weathering and finally, atmospheric CO₂ is certainly model dependent essentially because of our limited knowledge of Devonian vegetation. For example, in an atmospheric general circulation model (here FOAM), the runoff pattern is significantly affected by topography, a factor very difficult to constrain in deep times notably due to the absence of proxies. For example the Eastern/Western Laurussia remains not vegetated in the Late Devonian (LP3) because humid air masses tend to be locked on the central orogenic belt of the Laurussia. A second limitation of our study concerns obviously the physical properties used in SLAVE-Jr to mimic the Devonian vegetation. We are aware that representing the complexity of the vegetation by one or two PFT, for a flora which was already much diversified during the Devonian, is a crude simplification. However we motivate this choice by the necessity to distinguish fundamental steps in the land-plant invasion, with a time resolution coherent with the other boundary conditions we used, notably the geography.

## 5. Conclusion

In the present contribution, we explore the impact of the land plant evolution during the Devonian on the evolution of the global carbon cycle and climate. The numerical model used (GEOCLIM) couples a 3D-climate model to a model of the global biogeochemical cycles and to a simple model for vegetation dynamics adapted to Devonian conditions. More specifically, we quantify the impact of the changing continental conditions during terrestrialization (albedo, roughness, thermal properties and potential evaporation) on the climate and hence on the CO₂ consumption through continental silicate weathering; we also account for the continental drift.

We found that, in response to the continental drift, the atmospheric CO₂ decreases from 6300 ppmv to 3950 ppmv from the Early to the Middle Devonian. This decrease is triggered by the drifting of Laurussia and Siberia into more humid latitudinal belts over the course of the Devonian. The Late Devonian pCO₂ is further reduced to 1800 ppmv by the spreading of the land plants on continental surfaces. This further decrease is mainly triggered by decreasing continental albedo, which warms the continental surfaces and hence promotes CO₂ consumption by silicate rock weathering.

A salient result of our simulations is that the evolutionary development and the spread of land plants during the Devonian have been an important controlling factor of the Paleozoic climatic evolution. We show that continental surface changes induced by land plants and tectonic drift have produced a large CO₂ drawdown without being associated to a global cooling, because the cooling trend is counteracted by a warming trend resulting from the surface albedo reduction. Future works should now integrate our results to the previously modeled effects of the land plant expansion during the Devonian, as the impact of plant rooting on the weathering rates, and the burial of terrestrial organic carbon.

## Acknowledgements

This work has been funded by the French National Agency for Research, project ANR-06-BLAN-0347-01. Y.D. thanks the ANR project

<footnote n="1">¹ 1.85 is reached for a pCO₂ decrease from 3954 to 2125 ppm (Table 2), and only 1.25 for a drop from 3954 to 3161 ppmv of CO₂. These two end-members represent the CO₂ reduction predicted for an albedo reduction from 0.32 to respectively 0.12 and 0.20.</footnote>
COLORS (no. ANR-09-JCJC-0105). We thank Pierre Friedlingstein for providing the original SLAVE model. AMAP (Botany and Computational Plant Architecture) is a joint research unit which associates CIRAD (UMR51), CNRS (UMR5120), INRA (UMR931), IRD (R123), and Montpellier 2 University (UM27); http://amap.cirad.fr/.

We gratefully acknowledge Thomas Algeo and an anonymous reviewer for providing constructive comments on this manuscript, and Peter DeMenocal for the quality of the editorial procedure. This is IPGP contribution no. 3218.

## Appendix A

The Land Cover Change (LCC) is determined as the ΔTemp/ΔSNFSC where LP3 is the simulation of reference. The LCC being computed with an unchanged pCO₂, all results are obtained with a pCO₂ fixed at 10PAL (2800 ppmv). All simulations are done with LP3 boundary conditions, except the vegetation albedo. SNFSC represents the solar net flux absorbed at the surface in clear sky diagnostic.

<table><thead><tr><th>Simulations</th><th>LP3</th><th>LP3 - scenario 1</th><th>LP3 - scenario 2</th></tr></thead><tbody><tr><td>Vegetal albedo (vis/IR)</td><td>0.06/0.18</td><td>0.14/0.25</td><td>0.09/0.32</td></tr><tr><td>Continental albedo (%)</td><td>17.2</td><td>22.7</td><td>23.5</td></tr><tr><td>SNFSC (W/m²)</td><td>199.5</td><td>190.4</td><td>190.0</td></tr><tr><td>Continentale temp (°C)</td><td>18.3</td><td>16.8</td><td>16.6</td></tr><tr><td>LCC (°C/Wm<sup>-2</sup>)</td><td>-</td><td>0.17</td><td>0.18</td></tr></tbody></table>

## References

Algeo, T.J., Scheckler, S.E., 1998. Terrestrial-marine teleconnections in the Devonian: links between the evolution of land plants, weathering processes, and marine an-oxic events. Philos. Trans. R. Soc. B-Biol. Sci. 353 (1365), 113–128.

Algeo, T.J., Scheckler, S.E., Maynard, B., 2001. Effects of the Middle to Late Devonian spread of vascular land plants on weathering regimes, marine biotas, and global climate. In: Gensel, P.G., Edwards, D. (Eds.), Plants Invade the Land: Evolutionary and Environmental Perspectives, 213–236. Columbia University Press, New York.

Anderson, H.M., Hiller, N., Gess, R.W., 1995. Archeopteris (progymnospermopsida) from the Devonian of Southern Africa. Bot. J. Linn. Soc. 117 (4), 305–320.

Arnold, C.A., 1930. The genus *Callixylon* from the Upper Devonian of central and western New York. T. Pap. Mic. Acad. Sci. 11, 1–50.

Beck, C.B., 1964. Predominance of *Archeopteris* in Upper Devonian flora of western Catskills and adjacent Pennsylvania. Bot. Gaz. 125 (2), 126–128.

Beck, C.B., 1971. On the anatomy and morphology of lateral branch systems of *Archeopteris*. Am. J. Bot. 58 (8), 758–784.

Beerling, D.J., Berner, R.A., 2005. Feedbacks and the coevolution of plants and atmospheric CO₂. Proc. Natl. Acad. Sci. U. S. A. 102 (5), 1302–1305.

Beerling, D.J., Osborne, C.P., Chaloner, W.G., 2001. Evolution of leaf-form in land plants linked to atmospheric CO₂ decline in the Late Palaeozoic era. Nature 410 (6826), 352–354.

Berner, R.A., 1987. Models for carbon and sulfur cycles and atmospheric oxygen – application to Palaeozoic geologic history. Am. J. Sci. 287 (3), 177–196.

Berner, R.A., 1991. A model for atmospheric CO₂ over Phanerozoic time. Am. J. Sci. 291 (4), 339–376.

Berner, R.A., 1994. GEOCARBII – a revised model of atmospheric CO₂ over Phanerozoic time. Am. J. Sci. 294 (1), 56–91.

Berner, R.A., 2001. The effect of the rise of land plants on atmospheric CO₂ during the Paleozoic. In: Gensel, P.G., Edwards, D. (Eds.), Plants Invade the Land: Evolutionary and Environmental Perspectives, 173–178. Columbia University Press, New York.

Berner, R.A., 2004. The Phanerozoic Carbon Cycle: CO₂ and O₂. Oxford University Press, p. 158.

Berner, R.A., 2006. GEOCARBSURF: a combined model for Phanerozoic atmospheric O₂ and CO₂. Geochim. Cosmochim. Acta 70, 5653–5664.

Berner, R.A., Kothavala, Z., 2001. GEOCARB III: a revised model of atmospheric CO₂ over Phanerozoic time. Am. J. Sci. 301 (2), 182–204.

Blakey, R.C., 2003. Carboniferous–Permian paleogeography of the assembly of Pangea. In: Wong, T.E. (Ed.), XVth International Congress on Carboniferous and Permian Stratigraphy. Royal Netherlands Academy of Arts and Sciences, Utrecht, the Netherlands, pp. 443–465.

Bormann, B.T., Wang, D., Bormann, F.H., Benoit, G., April, R., Snyder, M.C., 1998. Rapid, plant-induced weathering in an aggrading experimental ecosystem. Biogeochemistry 43 (2), 129–155.

Cressler III, W.L., 2006. Plant paleoecology of the Late Devonian Red Hill locality, north- central Pennsylvania, an *Archeopteris*-dominated wetland plant community and early tetrapot site. In: Greb, S.F., DiMichele, W.A. (Eds.), Geological Society of America, Boulder, CO, Special Paper, Wetlands through Time (399), pp. 79–102.

Cressler III, W.L., Daeschler, E.S., Slingerland, R., Peterson, D.A., 2010. Terrestrialization in the Late Devonian: a palaeoecological overview of the Red Hill site, Pennsylvanian, USA. In: Vecoli, M., Clément, G., Meyer-Berthaud, B. (Eds.), Geological Society, London, Special Publications, The Terrestrialization Process: Modelling Complex Interactions at the Biosphere–Geosphere Interface (339), pp. 111–128.

Dessert, C., Dupre, B., Francois, L.M., Schott, J., Gaillardet, J., Chakrapani, G., Bajpai, S., 2001. Erosion of Deccan Traps determined by river geochemistry: impact on the global climate and the Sr-87/Sr-86 ratio of seawater. Earth Planet. Sci. Lett. 188 (3–4), 459–474.

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

Donnadieu, Y., Godderis, Y., Ramstein, G., Nedelec, A., Meert, J., 2004. A 'snowball Earth' climate triggered by continental break-up through changes in runoff. Nature 428 (6980), 303–306.

Donnadieu, Y., Godderis, Y., Pierrehumbert, R., Dromart, G., Fluteau, F., Jacob, R., 2006. A GEOCLIM simulation of climatic and biogeochemical consequences of Pangea breakup. Geochem. Geophys. Geosyst. 7, Q10109.

Donnadieu, Y., Godderis, Y., Bouttes, N., 2009. Exploring the climatic impact of the continental vegetation on the Mezoso atmospheric CO₂ and climate history. Clim. Past 5 (1), 85–96.

Driese, S.G., Mora, C.I., 2001. Diversification of Siluro-Devonian plant traces in paleosols and influence on estimates of paleoatmospheric CO₂ levels. In: Gensel, P.G., Edwards, D. (Eds.), Plants Invade the Land: Evolutionary and Environmental Perspectives. Columbia University Press, New York, pp. 237–253.

Driese, S.G., Mora, C.I., Elick, J.M., 1997. Morphology and taphonomy of root and stump casts of the earliest trees (Middle to Late Devonian, Pennsylvania and New York, USA. Palaios 12 (6), 524–537.

Edwards, D., Wellman, C.H., 2001. Embryophytes on land: the Ordocvician to Lochkovian (Lower Devonian) record. In: Gensel, P.G., Edwards, D. (Eds.), Plants Invade the Land: Evolutionary and Environmental Perspectives, Columbia, pp. 3–28.

Ekart, D.D., Cerling, T.E., Montanez, I.P., Tabor, N.J., 1999. A 400 million year carbon isotope record of pedogenic carbonate: implications for paleoatmospheric carbon dioxide. Am. J. Sci. 299 (10), 805–827.

Fairon-Demaret, M., Hartkopf-Fröder, C., 2004. Late Famennian plant mesofossils from the Refrath 1 Borehole (Bergisch Gladbach-Paffrath Syncline; Ardennes-Rhenish Massif, Germany). Cour. Forsch.-Inst. Senckenberg 251, 89–121.

Fairon-Demaret, M., Leponce, I., Streel, M., 2001. *Archeopteris* from the Upper Famennian of Belgium: heterospority, nomenclature, and palaeobiogeography. Rev. Palaeobot. Palynol. 115 (1–2), 79–97.

Falcon-Lang, H.J., 2000. The relationship between leaf longevity and growth ring markedness in modern conifer woods and its implications for palaeoclimatic studies. Palaeogeogr. Palaeoclimatol. Palaeocoeol. 160 (3–4), 317–328.

Freidlingstein, P., 1995. Modélisation du cycle du carbone biosphérique et étude du couplage biosphère-atmosphère, de Brussel.

Friedlingstein, P., Delire, C., Muller, J.F., Gerard, J.C., 1992. The climate induced variation of the continental biosphere – a model simulation of the last glacial maximum. Geophys. Res. Lett. 19 (9), 897–900.

Gaillardet, J., Dupre, B., Louvat, P., Allegre, C.J., 1999. Global silicate weathering and CO₂ consumption rates deduced from the chemistry of large rivers. Chem. Geol. 159 (1–4), 3–30.

Galtier, J., Paris, F., ElAouadDebbaj, Z., 1996. The occurrence of Callixylon in the upper Devonian of Morocco and its palaeogeographical significance. C. R. Acad. Sci. Ser. Ila: Fascicule Sci. Terre Planets 322 (10), 893–900.

Godderis, Y., Francois, L.M., 1995. The Cenozoic evolution of the strontium and carbon cycles: relative importance of continental erosion and mantle exchanges. Chem. Geol. 126 (2), 169–190.

Godderis, Y., Joachimski, M.M., 2004. Global change in the Late Devonian: modelling the Frasnian–Famennian short-term carbon isotope excursions. Palaeogeogr. Palaeoclimatol. Palaeocoeol. 202 (3–4), 309–329.

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

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

Greb, S.F., Dimichele, W.A., Gastaldo, R.A., 2006. Evolution and importance of wetlands in earth history. Geol. Soc. Am. 399, 1–40 (special paper).

Hack, J.J., Boville, B.A., Briegleb, B.P., Kiehl, J.T., Rasch, P.J., Williamson, D.L., 1993. Description of the NCAR Community Climate Model (CCM2). NCAR Technical Note NCAR/TN-382+STR: 108.

Hotton, C., Hueber, F., Griffing, D., Bridge, J., 2001. Early terrestrial plant environments: an example from the Emsian of Gaspé, Canada. In: Gensel, P.G., Edwards, D. (Eds.), Plants Invade the Land: Evolutionary and Environmental Perspectives, 179–212. Columbia University Press, New York, pp. 179–212.

Jacob, R., 1997. Low frequency variability in a simulated atmosphere ocean system. Joachimski, M.M., Breisig, S., Buggisch, W., Talent, J.A., Mawson, R., Gereke, M., Morrow, J.R., Day, J., Weddige, K., 2009. Devonian climate and reef evolution: insights from oxygen isotopes in apatites. Earth Planet. Sci. Lett. 284 (3–4), 599–609.

Kenrick, P., Crane, P.R., 1997. The origin and early evolution of plants on land. Nature 389 (6646), 33–39.

Kiehl, J.T., Dickinson, R.E., 1987. A study of the radiative effects of enhanced atmospheric CO₂ and CH₄ on early Earth surface temperatures. J. Geophys. Res. Atmos. 92 (D3), 2991–2998.
Le Hir, G., Donnadieu, Y., Godderis, Y., Pierrehumbert, R.T., Halverson, G.R., Macouin, M., Nedelec, A., Ramstein, G., 2009. The snowball Earth aftermath: exploring the limits of continental weathering processes. Earth Planet. Sci. Lett. 277 (3–4), 453–463.

Lemoigne, Y., Iurina, A., Snigirevskaya, N., 1983. Révision du genre *Callixylon* Zalessky 1911 (*Archaeopteris*) du Dévonien. Palaeontographica 186B, 81–120.

Matthews, E., 1984. Prescription of land-surface boundary-conditions in GISS GCM II: a simple method based on fine-resolution data bases. NASA Tech. Memo. 86096, 1–21 June.

Meyer-Berthaud, B., Wendt, J., Galtier, J., 1997. First record of a large Callixylon trunk from the late Devonian of Gondwana. Geol. Mag. 134 (6), 847–853.

Meyer-Berthaud, B., Scheckler, S.E., Bousquet, J.L., 2000. The development of *Archaeopteris*: new evolutionary characters from the structural analysis of an early Famenian trunk from southeast Morocco. Am. J. Bot. 87 (4), 456–468.

Meyer-Berthaud, B., Rücklin, M., Soria, A., Belka, Z., Lardeux, H., 2004. Frasnian plants from the Dra Valley, southern Anti-Atlas, Morocco. Geol. Mag. 141 (6), 675–686.

Meyer-Berthaud, B., Soria, A., Decombeix, A.-L., 2010. The land plant cover in the Devonian: a reassessment of the evolution of the tree habit. In: Vecoli, M., Clément, G., Meyer-Berthaud, B. (Eds.), The Terrestrialization Process: Modelling Complex Interactions at the Biosphere-geosphere Interface. : Special Publications, 339. Geological Society, London, pp. 59–70.

Mintz, J.S., Driese, S.G., White, J.D., 2010. Environmental and ecological variability of middle Devonian (Givetian) forests in Appalachian basin paleosols, New York, United States. Palaios 25 (1–2), 85–96.

Mora, C.I., Driese, S.G., Colarusso, L.A., 1996. Middle to late Paleozoic atmospheric CO₂ levels from soil carbonate and organic matter. Science 271 (5252), 1105–1107.

Moulton, K.L., Berner, R.A., 1998. Quantification of the effect of plants on weathering: studies in Iceland. Geology 26 (10), 895–898.

Nardin, E., Goddéris, Y., Donnadieu, Y., Le Hir, G., Blakey, R.C., Pucéat, E., Aretz, M., 2010. Modeling the early Paleozoic long-term climatic trend. Geol. Soc. Am. Bull. 123 (5–6), 1181–1192.

Oliva, P., Dupre, B., Martin, F., Viers, J., 2004. The role of trace minerals in chemical weathering in a high-elevation granitic watershed (Estibere, France): chemical and mineralogical evidence. Geochim. Cosmochim. Acta 68 (10), 2223–2243.

Retallack, G.J., 1997. Early forest soils and their role in Devonian global change. Science 276 (5312), 583–585.

Scheckler, S.E., 1986. Ancestral character states of Archeopteridales (progymnospermopsida). Am. J. Bot. 73 (5), 704–705.

Scheckler, S.E., 2001. Afforestation — the first forests. In: Briggs, D.E.G., Crowther, P. (Eds.), Palaeobiology, II. Blackwell Science Ltd, Oxford, pp. 67–71.

Scotese, C.R., Boucot, A.J., McKerrow, W.S., 1999. Gondwanan palaeogeography and palaeoclimatology. J. Afr. Earth Sci. 28 (1), 99–114.

Snigirevskaya, N.S., Snigirevsky, S.M., 2001. New locality of *Callixylon* (*Archaeopteridae*) in the Late Devonian of Andoma Mountain (Vologda Region, north-west Russia) and its importance for the reconstruction of archaeopterid distribution. Acta Palaeobotanica 41, 97–105.

Steemans, P., Wellman, C.H., Gerrienne, P., 2010. Palaeogeographic and palaeoclimatic considerations based on Ordovician to Lochkovian vegetation. In: Vecoli, M., Clément, G., Meyer-Berthaud, B. (Eds.), The Terrestrialization Process: Modelling Complex Interactions at the Biosphere-geosphere Interface. : Geological Society, London, Special Publications. The Geological Society of London, pp. 47–56.

Stein, W.E., Mannolini, F., Hernick, L.V., Landing, E., Berry, C.M., 2007. Giant cladoxylop- sid trees resolve the enigma of the Earth's earliest forest stumps at Gilboa. Nature 446 (7138), 904–907.

Taylor, E.L., Taylor, T.N., Cuneo, N.R., 1992. The present is not the key to the past – a polar forest from the Permian of Antarctica. Science 257 (5077), 1675–1677.

Van Geldern, R., Joachimski, M.M., Day, J., Jansen, U., Alvarez, F., Yolkin, E.A., Ma, X.P., 2006. Carbon, oxygen and strontium isotope records of Devonian brachiopod shell calcite. Palaeogeogr. Palaeoclimatol. Palaeoecol. 240 (1–2), 47–67.

Veizer, J., Ala, D., Azmy, K., Bruckschen, P., Buhl, D., Bruhn, F., Carden, G.A.F., Diener, A., Ebneth, S., Godderis, Y., Jasper, T., Korte, C., Pawellek, F., Podlaha, O.G., Strauss, H., 1999. Sr<sup>87</sup>Sr/<sup>86</sup>Sr, δ<sup>13</sup>C and δ<sup>18</sup>O evolution of Phanerozoic seawater. Chem. Geol. 161, 59–88.

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

Wellman, C.H., Gray, J., 2000. The microfossil record of early land plants. Philos. Trans. R. Soc. London, Ser. B-Biol. Sci. 355 (1398), 717–731.

Wellman, C.H., Habgood, K., Jenkins, G., Richardson, J.B., 2000. A new plant assemblage (microfossil and megafossil) from the Lower Old Red Sandstone of the Anglo-Welsh Basin: its implications for the palaeoecology of early terrestrial ecosystems. Rev. Palaeobot. Palynol. 109 (3–4), 161–196.

Yapp, C.J., Poths, H., 1992. Ancient atmospheric CO₂ pressures inferred from natural goethites. Nature 355 (6358), 342–344.
<figure><img src="image_58.png" /></figure>

# The faint young Sun problem revisited with a 3-D climate–carbon model – Part 1

G. Le Hir¹, Y. Teitler¹,*, F. Fluteau¹, Y. Donnadieu², and P. Philippot¹

¹IPGP, Institut de Physique du Globe de Paris, Sorbonne Paris cité, Université Paris Diderot, UMR7154, 75005, France

²LSCE, CNRS-CEA-UVSQ, 91191 Gif-sur-Yvette, France

*now at: Centre for Exploration Targeting, University of Western Australia, Crawley, Australia

Correspondence to: G. Le Hir (lehir@ipgp.fr)

Received: 6 February 2013 – Published in Clim. Past Discuss.: 22 March 2013

Revised: 17 February 2014 – Accepted: 18 February 2014 – Published: 3 April 2014

**Abstract.** During the Archaean, the Sun's luminosity was 18 to 25% lower than the present day. One-dimensional radiative convective models (RCM) generally infer that high concentrations of greenhouse gases (CO₂, CH₄) are required to prevent the early Earth's surface temperature from dropping below the freezing point of liquid water and satisfying the faint young Sun paradox (FYSP, an Earth temperature at least as warm as today). Using a one-dimensional (1-D) model, it was proposed in 2010 that the association of a reduced albedo and less reflective clouds may have been responsible for the maintenance of a warm climate during the Archaean without requiring high concentrations of atmospheric CO₂ (*pCO₂*). More recently, 3-D climate simulations have been performed using atmospheric general circulation models (AGCM) and Earth system models of intermediate complexity (EMIC). These studies were able to solve the FYSP through a large range of carbon dioxide concentrations, from 0.6 bar with an EMIC to several millibars with AGCMs. To better understand this wide range in *pCO₂*, we investigated the early Earth climate using an atmospheric GCM coupled to a slab ocean. Our simulations include the ice-albedo feedback and specific Archaean climatic factors such as a faster Earth rotation rate, high atmospheric concentrations of CO₂ and/or CH₄, a reduced continental surface, a saltier ocean, and different cloudiness. We estimated full glaciation thresholds for the early Archaean and quantified positive radiative forcing required to solve the FYSP. We also demonstrated why RCM and EMIC tend to overestimate greenhouse gas concentrations required to avoid full glaciations or solve the FYSP. Carbon cycle–climate interplays and conditions for sustaining *pCO₂* will be discussed in a companion paper.

## 1 Introduction

Sagan and Mullen (1972) argued from solar models that, in the early Archaean, the young Sun's luminosity was 25% lower compared to now. According to climate models, the weakness of the young Sun's brightness would have implied that Earth's surface temperature was well below the freezing point of water during the first two billion years of its history. A completely frozen Earth, however, conflicts with evidence of the presence of liquid water at the surface and life development during the Archaean. More intriguingly, geological evidence suggests that the Archaean climate was similar or even warmer than the present-day climate. Greenhouse gases, and more specifically carbon dioxide, were soon put forward as good candidates for solving this issue known as the “faint young Sun paradox” (FYSP; see Feulner, 2012 for a detailed review). Using a vertical one-dimensional radiative convective climate model (RCM), Kasting (1984, 2010) calculated that an atmospheric carbon dioxide partial pressure (*pCO₂*) of 0.3 bar or ~1000 times the preindustrial atmospheric level (PAL) (i.e. 0.28 mbar or 280 ppm) was needed to maintain a mean global surface temperature of 15 °C for a Sun 75% dimmer than today.

In parallel with the modelling approach, atmospheric carbon dioxide concentration was semi-quantitatively estimated from Archaean palaeosols. Considering the absence of siderite in 2.8 Ga-old palaeosols, Rye et al. (1995) suggested from thermodynamical modelling that the *pCO₂* had not exceeded ~0.04 bar (or ~130 PAL) during the Late Archaean. Questioning the relevance of the Rye et al. (1995) approach, Sheldon (2006) conducted mass-balance calculations of weathering in several Archaean to Palaeoproterozoic
palaeosols and proposed that the *pCO₂* fluctuated from 3 mbar to 26 mbar (i.e. 11 to 93 PAL) from 2.2 to 1.0 Ga. More recently, Driese et al. (2011), using palaeoweather-ing indicators, suggested moderate CO₂ partial pressures of about 10 mbar at 2.7 Ga ago. Given the composition of iron-rich minerals found in Archaean banded iron forma-tions (BIFs), Rosing et al. (2010) argued that the coexis-tence of magnetite (i.e. Fe₃O₄, a moderately oxidised min-eral) and siderite (i.e. FeCO₃, a reduced mineral) in BIF was dependent on the atmospheric level in di-hydrogen (pH₂). Di-hydrogen (H₂) is a nutrient for methanogenic bacteria, and biological constraints impose a minimum pH₂. Using this constraint, Rosing et al. (2010) defined an upper limit for *pCO₂* of about 0.9 mbar (900 ppmv or 3.2 PAL). Such a low *pCO₂* has been the core of vivid debate among spe-cialists, notably because BIFs are formed in marine envi-ronments far from thermodynamic equilibrium with the at-mosphere and therefore are poor proxies of the early Earth atmosphere composition (Dauphas and Kasting, 2011; Rein-hard and Planavsky, 2011).

Given the difficulties in estimating the partial pressure of atmospheric greenhouse gases on the primitive Earth, other mechanisms influencing climate such as changes in cloud properties have been explored. The mean forcing of clouds results from two opposite effects (Pierrehumbert, 2010). Clouds are highly reflective in the spectrum of visible solar radiation (shortwave radiation) and highly absorbent in the spectrum of thermal terrestrial radiation (longwave radiation). This leads to two well-known opposite effects; the shortwave reflection effect contributing to the planetary albedo and the longwave absorption effect contributing to the greenhouse. These mechanisms are related to the cloud type, the greenhouse effect being stronger than the albedo effect in the case of cirrus clouds, which are quite trans-parent at visible wavelengths but good absorbers of ther-mal terrestrial radiation. Using a one-dimensional radiative model, Goldblatt and Zahnle (2011b) estimated that in or-der to compensate for a weak Sun, Earth must have been en-tirely cloud-covered with poorly reflective (i.e. absence of stratus) but 3.5 times thicker clouds than today. In the ab-sence of a mechanism able to change cloud properties in such proportions, the authors concluded that this scenario was highly unlikely. Goldblatt and Zahnle (2011b) also chal-lenged Rosing et al.'s (2010) hypothesis that thinner clouds, caused by sparse cloud condensation nuclei (CCN), could have acted as an alternative to atmospheric greenhouse gases to warm up the early Earth's surface. In the present-day atmo-sphere, abundant CCN permit the formation of small liquid water droplets, leading to optically dense and highly reflec-tive clouds because of their direct scattering effect (Charlson et al., 1987). During the Archaean, the concentration of CCN would have been about 10 to 20% of its natural present-day level due to lower organic and inorganic sources. As a consequence, fewer droplets, but of larger size, would have formed optically thin clouds characterised by a low scattering

of shortwave radiations and shorter lifetime. This hypothesis had been quantified with RCM simulations (Rosing et al., 2010). However, using a more realistic treatment of clouds, Goldblatt and Zahnle (2011b) demonstrated that the radia-tive forcing due to such changes in cloud properties has been overestimated.

To go one step further in the understanding of the FYSP, Kienert et al. (2012) conducted the first 3-D climate mod-elling of the Archaean climate using the EMIC (Earth sys-tem models of intermediate complexity) CLIMBER, which consists in an ocean general circulation model coupled with a simplified, low-resolution atmosphere. Using Archaean boundary conditions including a reduced incoming solar ra-diation, a faster Earth rotation and a reduced continental sur-face, CLIMBER simulates a present-day climate for a car-bon dioxide partial pressure fixed at 0.6 bar. The authors ex-plained their results by the ice-albedo feedback, which is en-hanced during the Archaean due to a faster Earth rotation rate. This result contrasts with conclusions obtained by Jenk-ins (1993b) and Jenkins et al. (1993) with CCMAO (CC-MAO was the version of the NCAR model at the begin-ning of the 90s, an atmospheric GCM coupled to a swamp ocean without heat capacity), where a faster rotation rate would have induced a warmer climate. The recent mod-elling of Archaean climate by Charney et al. (2013) using LMDz (an atmospheric GCM coupled to a slab ocean de-veloped by Laboratoire de Météorologie Dynamique) chal-lenged the CLIMBER results (Kienert et al., 2012, 2013) and revealed a very limited impact of the faster Earth ro-tation. However, LMDz results suggest that, for 3.8 Ga-old conditions, 100 mbar of carbon dioxide and 2 mbar of CH₄ are required to maintain an Earth climate as warm as today. Similar results are predicted for the Late Archaean (2.7 Ga), with the atmospheric GCM coupled to a slab ocean CAM3 (Wolf and Toon, 2013), where a combination of 20 mbar of CO₂ and 1 mbar of CH₄ is required to yield a global mean surface temperature near that of the present day.

To understand why AGCMs (Charney et al., 2013; Wolf and Toon, 2013; Kunze et al., 2014), EMIC (Kienert et al., 2012, 2013) and RCMs (Rosing et al., 2010; Kasting, 1984) differ in simulating the early Earth climate, we investigated the Precambrian climate using an atmospheric general circu-lation model coupled to a slab ocean FOAM. The paper is organised as follows. First, we estimate the radiative deficit required for the onset of a pan-glaciation. We run experi-ments ranging from 3.5 to 1 Ga to explore the impacts of Archaean boundary conditions (solar constant, palaeogeogra-phy, amount of emerged continental surface, and cloud prop-erties). Second, we investigate the faint young Sun problem by performing additional sets of simulations where the ef-fects of greenhouse gases (CO₂ and CH₄), saltier oceans, faster Earth rotation or sea ice albedo are tested. Finally, we compare our results with the above-mentioned studies to explore the causes of observed differences among RCMs,
<figure><img src="image_59.png" /><figcaption>Fig. 1. Net long-wave fluxes computed for CO₂ and CH₄ with the FOAM radiative code (NCAR ccm3) (the total atmospheric is held constant, 1 bar). Thick arrows indicate the radiative forcing corresponding to 0.9 mbar of CO₂ and CH₄. Using the present-day atmosphere as a reference, the change in net flux at the tropopause (ΔF_trop) is 19 W m⁻², which reduces to 6 and 13 W m⁻² for 0.9 mbar of CO₂ and CH₄, respectively.</figcaption></figure>

EMIC and AGCMs, and discuss how these models differ in representing the ice-albedo feedback mechanism.

# 2 Model description and simulation specifications

## 2.1 Model description

Climate simulations were run with the atmospheric general circulation climate model (AGCM) FOAM 1.5. The atmospheric component of FOAM is a parallelised version of NCAR's Community Climate Model 2 (CCM2) with the upgraded radiative and hydrologic physics incorporated in CCM3 v. 3.2 (Jacob, 1997). All simulations were performed with an R15 spectral resolution (4.5° × 7.5°) and 18 vertical levels. The radiative code (extracted from NCAR CCM3) has already shown its validity in atmospheres highly enriched in CO₂ (Pierrehumbert, 2004; Le Hir et al., 2010). The infrared absorption has been computed (Fig. 1) to replace CO₂-CH₄ conversion proposed by Rosing et al. (2010). Under an enriched CH₄ atmosphere, no significant differences are observed compared to prior estimates (Kiehl and Dickinson, 1987; Halevy et al., 2009). FOAM integrates the CCM3 version of cloud water parameterisation. As in the CCM3, cloud optical properties and liquid water droplet clouds used the Slingo (1989) parameterisation, where the liquid water droplet clouds (Cloud Water Path or CWP) is determined by

integrating the liquid water concentration ($q_{clw}$):

$$CWP = \int q_{clw} \cdot dz \quad \text{where} \quad q_{clw} = q_{clw0} e^{(-z/H_{clw})} \quad (1)$$

$z$ being the altitude, $q_{clw0}$ an empirical constant equals to 0.21 g m⁻³ and $H_{clw}$ the scale height representing the total precipitable water column, which is dependent of the precipitable water ($q_{tot}$). This last equation implies an exponential decay for in-cloud water concentration,

$$H_{clw} = 700 \ln(1 + q_{tot}) \quad (2)$$

with $q_{tot}$ being computed by integrating the specific humidity $Q$ (in kg kg⁻¹) over $dp$, $p$ being the pressure and $g$ the gravitational acceleration (9.8 m s⁻²):

$$q_{\text{tot}} = 1/g * \int Q \, dp. \quad (3)$$

FOAM is used in mixed-layer mode, meaning the atmospheric component is linked to a 50 m mixed-layer ocean with heat transport parameterised through diffusion. The diffusive rate of the slab ocean ($Q_{flux}$) is derived from the net surface fluxes and fixed to represent sea ice-covered seas based on estimates of the present-day sea ice mass budget. The sea ice module uses the thermodynamics of NCARs to compute ice fraction (growth and melting), snow cover, and penetrating radiation, and includes a Semtner model to calculate the sea ice temperature. The sea ice module does not treat sea ice dynamics, a lack which may slow the sea ice flow and underestimate the sea ice cover (Voigt and Abbot, 2012). Sea ice albedo is fixed at 0.7 (visible spectrum) and 0.5 (near-infrared spectrum).

## 2.2 Boundary conditions and experimental design

In this section, we describe boundary conditions used to explore the effects of solar constants, palaeogeography (amount of continental surface), and cloud properties from 3.5 to 1 Ga ago. Sets of simulations are summarised in Table 1.

### 2.2.1 Solar irradiance and palaeogeographies

Following Gough's formula (1981), the solar constant has been increased by 19%, from 3.5 Ga (1053 W m⁻²) to 1 Ga (1258 W m⁻²) (Fig. 2a). Orbital parameters were set at their present-day values (orbital eccentricity: 0.016724, obliquity: 23.4463, vernal equinox latitude of perihelion: 77.9610) in all experiments. We varied the continental surface in response to the crustal growth and emergence, a process which affected the mean surface albedo ($\alpha_{ocean} = 0.06$, $\alpha_{baresoil} = 0.32$). Different scenarios ranging from a progressive formation to a quasi-instantaneous formation of all continental crust have been proposed (see Lowe and Tice, 2007 for a detailed review). Here we considered a gradual continental growth (Fig. 2b), consistent with the changes of $^{87}$Sr/$^{86}$Sr