of the model simulations (at least several tens of million years), we assume this budget to be close to steady state:

$$ \frac{r_{\text{MOR}} - r_{\text{OC}}}{9.43 + r_{\text{MOR}}} F_{\text{MOR}}^{\text{Sr}} + \frac{r_{\text{CW}} - r_{\text{OC}}}{9.43 + r_{\text{CW}}} F_{\text{CW}}^{\text{Sr}} + \frac{r_{\text{swy}} - r_{\text{OC}}}{9.43 + r_{\text{swy}}} F_{\text{swy}}^{\text{Sr}} \\ + \frac{r_{\text{swo}} - r_{\text{OC}}}{9.43 + r_{\text{swo}}} F_{\text{swo}}^{\text{Sr}} \\ = 0 \qquad (A7) $$

where $r_{MOR}$, $r_{CW}$, $r_{swy}$, and $r_{swo}$ are the $^{87}$Sr/$^{86}$Sr ratio of the mantle, carbonates, young volcanics and old cratonic rocks. These are set to 0.703, 0.708, 0.703 and 0.718 respectively. $r_{OC}$ is the measured $^{87}$Sr/$^{86}$Sr ratio of seawater recorded in sedimentary carbonates (Veizer et al., 1999). The superscript Sr to the flux F means that we are considering the strontium fluxes associated to the carbon and alkalinity fluxes.

(5) Lastly, three additional equations describe the dependence of silicate ($F_{swy}$ and $F_{swo}$) and carbonate weathering ($F_{CW}$) on atmospheric CO$_2$, through their dependencies on runoff and air temperature. In GEOCARB-style models, only global runoff and mean average global temperature are calculated as a function of CO$_2$ (Eqs. (1) and (2)), and then used to calculate mean global weathering fluxes. In GEOCLIM, a physically-based 3D climate model is used to calculate the spatial distribution of runoff and temperature as a function of CO$_2$, which are subsequently used to reconstruct spatially-distributed weathering fluxes.

If the solid Earth degassing rate is known ($F_{vol}$ and $F_{MOR}$), the 9 Eqs. ((A1), (A2), (A4), (A5), (A6), (A7) + the three climatic dependences) can be solved for $F_{CW}$, $F_{OW}$, $F_{swy}$, $F_{swo}$, $\sigma_y$, $F_{cd}$, $F_{od}$, and atmospheric CO$_2$ (9 unknowns).

## References

Barron, E.J., Washington, W.M., 1982. Cretaceous climate: a comparison of atmospheric simulations with the geologic record. Palaeogeogr. Palaeoclimatol. Palaeoecol. 40, 103–133.

Barron, E.J., Hay, W.W., Thompson, S., 1989. The hydrologic cycle: a major variable during Earth history. Global Planet. Chang. 75, 157–174.

Barron, E.J., Fawcett, P.J., Pollard, D., Thompson, S., 1993. Model simulations of Cretaceous climates – the role of geography and carbon-dioxide. Philos. Trans. R. Soc. Lond. B Biol. Sci. 341, 307–315.

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

Berner, R.A., 1991. A model for atmospheric CO$_2$ over Phanerozoic time. Am. J. Sci. 291, 339–376.

Berner, R.A., 1994. 3GEOCARB II: a revised model of atmospheric CO$_2$ over Phanerozoic time. Am. J. Sci. 294, 56–91.

Berner, R.A., 2004. The Phanerozoic Carbon Cycle. Oxford University Press (150 pp.).

Berner, R.A., 2006. GEOCARBSULF: a combined model for Phanerozoic atmospheric O$_2$ and CO$_2$. Cosmochim. Acta 70, 5653–5664.

Berner, R.A., Caldeira, K., 1997. The need for mass balance and feedback in the geochemical carbon cycle. Geology 25, 955–956.

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

Bluth, G.J.S., Kump, L.R., 1991. Phanerozoic paleogeology. Am. J. Sci. 291, 284–308.

Craggs, H.J., Valdes, P.J., Widdowson, M., 2012. Climate model predictions for the latest Cretaceous: an evaluation using climatically sensitive sediments as proxy indicators. Palaeogeogr. Palaeoclimatol. Palaeoecol. 315–316, 12–23.

Crowley, T.J., Baum, S.K., 1991. Modelling late Paleozoic glaciation. Geology 20, 507–510.

Crowley, T.J., Baum, S.K., 1995. Reconciling Late Ordovician (440 Ma) glaciation with very high (14×) CO$_2$ levels. J. Geophys. Res. 100, 1093–1101.

Dessert, C., Dupré, B., Gaillardet, J., François, L.M., Allègre, C.J., 2003. Basalt weathering laws and the impact of basalt weathering on the global carbon cycle. Chem. Geol. 202, 267–273.

Donnadieu, Y., Goddéris, Y., Ramstein, G., Nédélec, A., Meert, J., 2004. A 'snowball Earth' climate triggered by continental break-up through changes in runoff. Nature 428, 303–306.

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

Donnadieu, Y., Goddéris, Y., Bouttes, N., 2009. Exploring the climatic impact of the continental vegetation on the Mesozoic atmospheric CO$_2$ and climate history. Clim. Past. 5, 85–96.

Dunne, T., 1978. Rate of chemical denudation of silicate rocks in tropical catchments. Nature 274, 244–246.

François, L.M., Goddéris, Y., 1998. Isotopic constraints on the Cenozoic evolution of the carbon cycle. Chem. Geol. 145, 177–212.

François, L.M., Walker, J.C.G., 1992. Modelling the Phanerozoic carbon cycle and climate: constraints from the $^{87}$Sr/$^{86}$Sr isotopic ratio of seawater. Am. J. Sci. 292, 81–135.

Gibbs, M.T., Kump, L.R., 1994. Global chemical erosion during the last glacial maximum and the present: sensitivity to changes in lithology and hydrology. Paleogeography 9, 529–543.

Gibbs, M.T., Barron, E.J., Kump, L.R., 1997. An atmospheric pCO$_2$ threshold for glaciation in the Late Ordovician. Geology 25, 447–450.

Gibbs, M.T., Bluth, G.J.S., Fawcett, P.J., Kump, L.R., 1999. Global chemical erosion over the last 250 My: variations due to changes in paleogeography, paleoclimate, and paleoecology. Am. J. Sci. 299, 611–651.

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

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

Goddéris, Y., Donnadieu, Y., de Vargas, C., Pierrehumbert, R.T., Dromart, G., van de Schootbrugge, B., 2008a. Causal or casual link between the rise of nannoplankton calcification and a tectonically-driven massive decrease in Late Triassic atmospheric CO$_2$? Earth Planet. Sci. Lett. 267, 247–255.

Goddéris, Y., Donnadieu, Y., Tombozafy, M., Dessert, C., 2008b. Shield effect on continental weathering: implication for climatic evolution of the Earth at the geological time-scale. Geoderma 145, 439–448.

Gyllenhaal, E.D., Engberts, C.J., Markwick, P.J., Smith, L.H., Patzkowsky, M.E., 1991. The Fujita–Ziegler model: a new semi-quantitative technique for estimating paleoclimate from paleogeographic maps. Palaeogeogr. Palaeoclimatol. Palaeoecol. 86, 41–66.

Hay, W.W., Barron, E.J., Thompson, S.L., 1990. Global atmospheric circulation experiments on an Earth with polar and tropical continents. J. Geol. Soc. London 147, 749–757.

Hay, W.W., DeConto, R.M., Wold, C.N., Wilson, K.M., Voight, S., Schulz, M., Wold-Rossby, A., Dullo, W.-C., Ronov, A.B., Balukhovsky, A.N., Soding, E., 1999. Alternative global Cretaceous paleogeography. In: Barrera, E., Johnson, C.C. (Eds.), Evolution of the Cretaceous ocean–climate system. Boulder, Colorado, Geological Society of America Special Paper, 332, pp. 1–48.

Hayes, J.M., Strauss, H., Kaufman, A.J., 1999. The abundance of $^{13}$C in marine organic matter and isotopic fractionation in the global biogeochemical cycle of carbon during the past 800 Ma. Chem. Geol. 161, 103–125.

Herrmann, A.D., Patzkowsky, M.E., Pollard, D., 2004. The impact of paleogeography, pCO$_2$, poleward ocean heat transport and sea level change on global cooling during the Late Ordovician: a climate model analysis. Palaeogeogr. Palaeoclimatol. Palaeoecol. 206, 59–74.

Horton, D.E., Poulsen, C.J., Pollard, D., 2010. Influence of high-latitude vegetation feedbacks on late Palaeozoic glacial cycles. Nat. Geosci. 3, 572–577.

Horton, D.E., Poulsen, C.J., Montañez, I.P., DiMichele, W.A., 2012. Eccentricity-paced late Paleozoic climate change. Palaeogeogr. Palaeoclimatol. Palaeoecol. 331–332, 150–161.

Hyde, W.T., Crowley, T.J., Tarasov, L., Peltier, W.R., 1999. The Pangean ice age: studies with a coupled climate–ice sheet model. Clim. Dyn. 15, 619–629.

Jacob, R., 1997. Low Frequency Variability in a Simulated Atmosphere Ocean System. Ph.D. Thesis University of Wisconsin, Madison, USA.

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 apatite. Earth Planet. Sci. Lett. 284, 599–609.

Kennett, J.P., 1977. Cenozoic evolution of Antarctic glaciation, the circum-Antarctic Ocean, and their impact on global paleoceanography. J. Geophys. Res. 82 (27), 3843–3860.

Kent, D.V., Muttoni, G., 2013. Modulation of Late Cretaceous and Cenozoic climate by variable drawdown of atmospheric pCO$_2$ from weathering of basaltic provinces on continents drifting through the equatorial humid belt. Clim. Past. 9, 525–546.

Köppen, W., Wegener, A., 1924. Die Klimate der geologischen Vorzeit. Gebrüder Bornträger, Berlin.

Korte, C., Kozur, H.W., Veizer, J., 2005. $^{δ}$13C and $^{δ}$18O values of Triassic brachiopods and carbonate rocks as proxies for coeval seawater and palaeotemperature. Palaeogeogr. Palaeoclimatol. Palaeoecol. 226, 287–306.

Kothavala, Z., Oglesby, R.J., Saltzman, B., 1999. Sensitivity of equilibrium surface temperature of CCM3 to systematic changes in atmospheric CO$_2$. Geophys. Res. Lett. 26, 209–212.

Kutzbach, J.E., Gallimore, R.G., 1989. Pangaean climates: megamonsoons of the megacontinent. J. Geophys. Res. 94, 3341–3357.

Labat, D., Goddéris, Y., Probst, J.-L., Guyot, J.-L., 2004. Evidence for global runoff increase related to climate warming. Adv. Water Resour. 27, 631–642.

Lagabrielle, Y., Goddéris, Y., Donnadieu, Y., Malavieille, J., Suarez, M., 2009. The tectonic history of Drake Passage and its possible impacts on global climate. Earth Planet. Sci. Lett. 279, 197–211.

Le Hir, G., Donnadieu, Y., Goddéris, Y., Meyer-Berthaud, B., Ramstein, G., Blakey, R.C., 2011. The climate change caused by the land plants invasion in the Devonian. Earth Planet. Sci. Lett. 310, 203–212.

Lefebvre, V., Donnadieu, Y., Sepulchre, P., Swingedouw, D., Zhang, Z.-S., 2012. Deciphering the role of southern gateways and carbon dioxide on the onset of the Antarctic Circumpolar Current. Paleogeography 27. http://dx.doi.org/10.1029/2012PA002345.

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

Longinelli, A., 2013. Comment on work by Pucéat et al. (2010) on a revised phosphate-water fractionation equation. Earth Planet. Sci. Lett. 377–378, 378–379.

Lunt, D.J., Haywood, A.M., Schmidt, G.A., Salzmann, U., Valdes, P.J., Dowsett, H.J., 2010. Earth system sensitivity inferred from Pliocene modelling and data. Nat. Geosci. 3, 60–64.
Lyell, C., 1830. Principles of Geology, Being an Attempt to Explain the Former Changes of the Earth's Surface by Reference to Causes Now in Operation. , Vol. 1. John Murray, London.

Marshall, H.G., Walker, J.C.G., Khun, W.R., 1988. Long-term climate change and the geochemical cycle of carbon. *J. Geophys. Res.* 93, 791–801.

Nardin, E., Goddéris, Y., Donnadieu, Y., Le Hir, G., Blakey, R.C., Pucéat, E., Aretz, M., 2011. Modelling the early Paleozoic long-term climatic trend. *Geol. Soc. Am. Bull.* 123, 1181–1192.

Oglesby, R., Park, J., 1989. The effect of precessional isolation changes on Cretaceous climate and cyclic sedimentation. *J. Geophys. Res.* 94, 14793–14816.

Oliva, P., Viers, J., Dupré, B., 2003. Chemical weathering in granitic environments. *Chem. Geol.* 202, 255–256.

Otto-Bliesner, B.L., 1995. Continental drift, runoff, and weathering feedbacks: implications from climate model experiments. *J. Geophys. Res.* 100, 11537–11548.

Park, J., Royer, D.L., 2011. Geologic constraints on the glacial amplification of Phanerozoic climate sensitivity. *Am. J. Sci.* 311, 1–26.

Peters, N.E., 1984. Evaluation of Environmental Factors Affecting Yields of Major Dissolved Ions of Streams in the United States. , 2228. U.S. Geological Survey Water Supply Paper (39 pp.).

Prokoph, A., Shields, G.A., Veizer, J., 2008. Compilation and time-series analysis of a marine carbonate δ¹⁸O, δ¹³C, ⁸⁷Sr/⁸⁶Sr and δ³⁴S database through Earth history. *Earth Sci. Rev.* 87, 113–133.

Pucéat, E., Joachimski, M.M., Bouilloux, A., Monna, F., Bonin, A., Motreuil, S., Morinière, P., Hénard, S., Mourin, J., Dera, G., Quesne, D., 2013. Reply on comment by Longinelli (2013) on a revised phosphate-water fractionation equation. *Earth Planet. Sci. Lett.* 377–378, 380–382.

Raymo, M.E., 1991. Geochemical evidence supporting T.C. Chamberlin's theory of glaciation. *Geology* 19, 344–347.

Raymo, M.E., Ruddiman, W.F., Froelich, P.N., 1988. Influence of late Cenozoic mountain building on ocean geochemical cycles. *Geology* 16, 649–653.

Sitch, S., Smith, B., Prentice, I.C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J.O., Levis, S., Lucht, W., Sykes, M.T., Thonicke, K., Venevsky, S., 2003. Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model. *Glob. Chang. Biol.* 9, 161–185.

Spicer, R.A., Ahlberg, A., Herman, A.B., Hofmann, C.-C., Raikevich, M., Valdes, P.J., Markwick, P.J., 2008. The Late Cretaceous continental interior of Siberia: a challenge for climate models. *Earth Planet. Sci. Lett.* 267, 228–235.

Sun, Y., Joachimski, M.M., Wignall, P.B., Yan, C., Chen, Y., Jiang, H., Wang, L., Lai, X., 2012. Lethally hot temperatures during the Early Triassic greenhouse. *Science* 338, 366–370.

Tajika, E., 2007. Long-term stability of climate and global glaciations throughout the evolution of the Earth. *Earth Planets Space* 59, 293–299.

Taylor, L.L., Banwart, S.A., Leake, J.R., Beerling, D.J., 2011. Modelling the evolutionary rise of the ectomycorrhiza on subsurface weathering environments and the geochemical carbon cycle. *Am. J. Sci.* 311, 369–403.

Taylor, L.L., Banwart, S.A., Valdes, P.J., Leake, J.R., Beerling, D.J., 2012. Evaluating the effects of terrestrial ecosystems, climate and carbon dioxide on weathering over geological time: a global-scale process-based approach. *Phil. Trans. R. Soc. B* 367, 565–582.

Trotter, J.A., Williams, I.S., Barnes, C.R., Lécuyer, C., Nicoll, R.S., 2008. Did cooling oceans trigger the Ordovician biodiversification? Evidence from conodont thermometry. *Science* 321, 550–554.

Veizer, J., Ala, D., Azmy, K., Bruckschen, P., Buhl, D., Bruhn, F., Carden, G.A.F., Diener, A., Ebneth, S., Goddéris, Y., Jasper, T., Korte, C., Pawellek, F., Podlaha, O.G., Strauss, H., 1999. ⁸⁷Sr/⁸⁶Sr, δ¹³C and δ¹⁸O evolution of Phanerozoic seawater. *Chem. Geol.* 161, 59–88.

Veizer, J., Goddéris, Y., François, L.M., 2000. Evidence for decoupling of atmospheric CO₂ and global climate during the Phanerozoic eon. *Nature* 408, 698–701.

Von der Heydt, A., Dijkstra, H.A., 2006. Effect of ocean gateways on the global ocean circulation in the late Oligocene and early Miocene. *Paleoceanography* 21. http://dx.doi.org/10.1029/2005PA001149.

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.

Wallmann, K., 2001. Controls on the Cretaceous and Cenozoic evolution of seawater composition, atmospheric CO₂ and climate. *Geochem. Cosmochim. Acta* 65, 3005–3025.

Wallmann, K., 2004. Impact of atmospheric CO₂ and galactic cosmic radiation on Phanerozoic climate change and the marine δ¹⁸O record. *Geochem. Geophys. Geosyst.* 5, Q06004.

Wegener, A., 1912. Die Entstehung der Kontinente. *Geol. Rundsch.* 3, 276–292.

Worsley, T.R., Kidder, D.L., 1991. First-order coupling of paleogeography and CO₂, with global surface temperature and its latitudinal contrast. *Geology* 19, 1161–1164.

Zhang, X., Prange, M., Steph, S., Butzin, M., Krebs, U., Lunt, D.J., Nisancioglu, K.H., Park, W., Schmitter, A., Schneider, B., Schulz, M., 2012. Changes in equatorial Pacific thermocline depth in response to Panamanian seaway closure: insights from a multi-model proxy. *Earth Planet. Sci. Lett.* 317–318, 76–84.
# Orbitally forced ice sheet fluctuations during the Marinoan Snowball Earth glaciation

Douglas I. Benn<fnref n="1" />,²*, Guillaume Le Hir³, Huiming Bao⁴, Yannick Donnadieu⁵, Christophe Dumas⁵, Edward J. Fleming<fnref n="1" />,⁶†, Michael J. Hambrey<fnref n="7" />, Emily A. McMillan⁶, Michael S. Petronis⁸, Gilles Ramstein⁵, Carl T. E. Stevenson⁶, Peter M. Wynn⁹ and Ian J. Fairchild⁶

Two global glaciations occurred during the Neoproterozoic. Snowball Earth theory posits that these were terminated after millions of years of frigidity when initial warming from rising atmospheric CO₂ concentrations was amplified by the reduction of ice cover and hence a reduction in planetary albedo<fnref n="1" />,². This scenario implies that most of the geological record of ice cover was deposited in a brief period of melt-back³. However, deposits in low palaeo-latitudes show evidence of glacial-interglacial cycles⁴⁻⁶. Here we analyse the sedimentology and oxygen and sulphur isotopic signatures of Marinoan Snowball glaciation deposits from Svalbard, in the Norwegian High Arctic. The deposits preserve a record of oscillations in glacier extent and hydrologic conditions under uniformly high atmospheric CO₂ concentrations. We use simulations from a coupled three-dimensional ice sheet and atmospheric general circulation model to show that such oscillations can be explained by orbital forcing in the late stages of a Snowball glaciation. The simulations suggest that while atmospheric CO₂ concentrations were rising, but not yet at the threshold required for complete melt-back, the ice sheets would have been sensitive to orbital forcing. We conclude that a similar dynamic can potentially explain the complex successions observed at other localities.

The Wilsonbreen Formation in northeast Svalbard contains a detailed record of environmental change during the Marinoan, the second of the major Cryogenian glaciations (650–635 Ma; refs 7,8). At this time, Svalbard was located in the Tropics on the eastern side of Rodinia⁹,¹⁰. The <180 m thick Wilsonbreen Formation was deposited within a long-lived intracratonic sedimentary basin¹¹. It is subdivided into three members (W1, W2 and W3) based on the relative abundance of diamictite and carbonate beds<fnref n="7" />,⁸ (Fig. 1 and Supplementary Figs 1 and 2). The occurrence throughout the succession of lacustrine sediments containing both precipitated carbonate and ice-rafted detritus, and intermittent evaporative carbonates and fluvial deposits, indicates that the basin remained isolated from the sea, consistent with eustatic sea-level fall of several hundred metres and limited local isostatic depression (Supplementary Information; ref. 12). This makes it ideal for investigating environmental change within a Neoproterozoic panglaciation, as it provides direct evidence of subaerial environments and climatic conditions.

We made detailed sedimentary logs at ten known and new localities extending over 60 km of strike (Fig. 1 and Supplementary Fig. 1; see Methods). Seven sediment facies associations were identified, recording distinct depositional environments that varied in spatial extent through time (Supplementary Fig. 3 and Supplementary Information). These are: FA1: subglacial, recording direct presence of glacier ice; FA2: fluvial channels; FA3: dolomitic floodplain, recording episodic flooding, evaporation and microbial communities; FA4: carbonate lake margin, including evidence of wave action; FA5: carbonate lacustrine, including annual rhythmites and intermittent ice-rafted debris; FA6: glacilacustrine, consisting of ice-proximal grounding-line fans (FA6-G) and ice-distal rainout deposits (FA6-D); and FA7: periglacial, recording cold, non-glacial conditions. Further descriptions are provided in the Supplementary Information. The vertical and horizontal distribution of these facies associations (Fig. 1) allows the sequence of environmental changes to be reconstructed in detail:

(1) The base of the Wilsonbreen Formation is a well-marked periglacial weathered horizon with thin wind-blown sands (Supplementary Fig. 4a,b). This surface records very limited sediment cycling in cold, arid conditions.

(2) At all localities, the weathering horizon is overlain by fluvial channel facies (FA2) and mudstones, marking the appearance of flowing water in the basin and implying positive air temperatures for at least part of the time (Supplementary Fig. 5a).

(3) Glacilacustrine deposits (FA6-D) record flooding of the basin and delivery of sediment by ice rafting (Supplementary Fig. 4c,d). Far-travelled clasts are common, indicating transport by a large, continental ice sheet.

(4) Warm-based, active ice advanced into the basin, indicated by traction tills and glacitectonic shearing (FA1; Supplementary Fig. 4e-g). (1-4 make up Member W1).

(5) Ice retreat is recorded by a second periglacial weathering surface (FA7). This is overlain by fluvial channel, floodplain, lake-margin and carbonate lacustrine sediments of W2 (FA2-5; Supplementary Fig. 5), recording a shifting mosaic of playa lakes and ephemeral streams. Lakes and river channels supported microbial communities. Millimetre-scale carbonate-siliciclastic rhythmites indicate seasonal cycles of

<footnote n="1">¹Department of Geology, The University Centre in Svalbard (UNIS), N-9171 Longyearbyen, Norway. ²School of Geography and Geosciences, University of St Andrews, St Andrews KY16 8YA, UK. ³Institut de Physique du Globe de Paris, 75238 Paris, France. ⁴Department of Geology and Geophysics, E235 Howe-Russell Complex, Louisiana State University, Baton Rouge, Louisiana 70803, USA. ⁵Laboratoire des Sciences du Climat et de l'Environnement, CNRS-CEA, 91190 Gif-sur-Yvette, France. ⁶School of Geography, Earth and Environmental Sciences, University of Birmingham, Birmingham B15 2TT, UK.</footnote>

<footnote n="7">⁷Institute of Geography and Earth Sciences, Aberystwyth University, Aberystwyth SY23 3DB, UK. ⁸Natural Resource Management, Environmental Geology, New Mexico Highlands University, Las Vegas, New Mexico 87701, USA. ⁹Lancaster Environment Centre, University of Lancaster, Lancaster LA1 4YQ, UK. ¹Present address: CASP, West Building, 181A Huntingdon Road, Cambridge CB3 0DH, UK. *e-mail: Doug.Benn@unis.no</footnote>
<figure><img src="image_98.png" /></figure>

**Figure 1 | Sedimentary architecture and palaeoenvironments of the Wilsonbreen Formation.** Regional correlation of facies associations and members W1, W2 and W3 across northeast Svalbard. From north to south, study locations are: DRA, Dracoisen; DIT, Ditlovtoppen; AND, East Andromedafjellet; REIN, Reinsryggen (informal name); KLO, Klofjellet; McD, MacDonaldryggen; BAC, Backlundtoppen-Kvitfjellet ridge; PIN, Pinnsvinryggen (informal name); SLA, Slangen and ORM, Ormen.

photosynthesis. The environment seems to have been closely similar to that of the present-day McMurdo Dry Valleys in Antarctica, although with less extreme seasonality owing to its low latitude¹³.

(6) Water levels and glacier extent underwent a series of oscillations, recorded by switches between glacilacustrine diamicite (FA6-D) and fluvial, lacustrine and lake-margin sediments (FA2-5) in W2. Sedimentation rates inferred from annual rhythmites in W2 suggest that each retreat phase may have lasted ~10⁴ years.

(7) A second major ice advance marks the base of W3, with widespread deposition of subglacial tills and glaciectonism of underlying sediments. Basal tills are absent from the northernmost locality, but close proximity of glacier ice is recorded by grounding-line fans (FA6-G; Supplementary Fig. 4h,i).

(8) Ice retreated while the basin remained flooded and glaciogenic sediment continued to be delivered to the lake by ice rafting. Thin laminated carbonates (FA5) in W3 indicate periods of reduced glaciogenic sedimentation, indicative of minor climatic fluctuations over timescales of ~10³ years (Supplementary Fig. 5g).

(9) A sharp contact with overlying laminated 'cap' carbonate (Supplementary Fig. 2) records the transition to post-glacial conditions. At some localities, basal conglomerates provide evidence of subaerial exposure, followed by marine transgression. The cap carbonate closely resembles basal Ediacaran carbonates elsewhere, and marks global deglaciation, eustatic sea-level rise and connection of the basin to the sea¹,¹²,¹⁴.

Environmental and atmospheric conditions during deposition of W2 and W3 can be further elucidated by isotopic data from carbonate-associated sulphate in lacustrine limestones (Fig. 2 and Supplementary Fig. 6). These exhibit negative to extremely negative Δ¹⁷O values, with consistent linear co-variation with δ³⁴S, indicating mixing of pre-glacial sulphate and isotopically light sulphate formed in a CO₂-enriched atmosphere¹⁵,¹⁶. The observed

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

**Figure 2 | Co-variation of Δ¹⁷O and δ³⁴S from carbonate-associated sulphate in W2 and W3.** 'Existing data' (ref. 16) and new data define a mixing line between pre-glacial sulphate (top) and an isotopically light sulphate formed by oxidation of pyrite including incorporation of a light-Δ¹⁷O signature from a CO₂-enriched atmosphere. Data from W2 and W3 lie on closely similar trend lines, indicating no detectable change in $p_{CO_2}$ between deposition of the two members.

values could reflect non-unique combinations of $p_{CO_2}$, $p_{O_2}$, $O_2$ residence time and other factors, but a box model¹⁷ indicates $p_{CO_2}$ was most likely ~10 to 100 mbar (1 mbar = 1,000 ppmv).

These values are far too high to permit formation of low-latitude ice sheets in the Neoproterozoic, but they are consistent with a late-stage Snowball Earth. For an ice-free Neoproterozoic Earth, model studies indicate mean terrestrial temperatures in the range 30–50 °C for $p_{CO_2}$ = 10 to 100 mbar (ref. 18). Formation of low-latitude ice sheets requires much lower $p_{CO_2}$, on the order of 0.1–1 mbar (refs 2,19,20). Once formed, however, ice sheets can persist despite rising CO₂ from volcanic outgassing, as a result of a high planetary albedo. This hysteresis in the relationship between $p_{CO_2}$ and planetary temperature is a key element of Snowball Earth theory. It implies that W2 and W3 were deposited relatively late in the Marinoan, after volcanic outgassing had raised $p_{CO_2}$ from 0.1 or 1 mbar to 10 or 100 mbar. Modelled silicate weathering and volcanic outgassing rates indicate that this would require 10⁶ to 10⁷ years²¹.

The consistent co-variation of Δ¹⁷O and δ³⁴S in lacustrine limestones in both W2 and W3 suggests no detectable rise in atmospheric $p_{CO_2}$, as this would alter the slope of the mixing line (Fig. 2). This implies that the glacier oscillations recorded in W2 and W3 occurred during a relatively short time interval (<10⁵ years²¹) towards the end of the Marinoan. In turn, this implies that the remainder of the Wilsonbreen Formation (including the basal weathering horizon) represents many millions of years, during which $p_{CO_2}$ built up from the low values necessary for inception of low-latitude glaciation to those indicated by the geochemical evidence. The weathering horizon provides direct evidence of cold, arid conditions during this interval, before the appearance of fluvial and glacilacustrine sediments in the basin.

The evidence for ice sheet advance/retreat cycles at low latitudes in a CO₂-enriched atmosphere motivated a series of numerical simulations to test the hypothesis that these cycles were linked to Milankovitch orbital variations. We employed asynchronous coupling of a three-dimensional ice sheet model and an atmospheric general circulation model using the continental configuration of ref. 22. We first ran simulations with a modern orbital configuration to examine ice sheet behaviour through a large range of $p_{CO_2}$ values from 0.1 to 100 mbar (ref. 23; Supplementary Figs 7–10). Consistently with previous results²,²⁰, at low $p_{CO_2}$ (0.1 mbar), global ice volume reaches 170 × 10⁶ km³, but substantial tropical land
<figure><img src="image_100.png" /></figure>

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

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

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

**Figure 3** | Modelled ice sheet oscillations in response to orbital forcing. **a,b,** Shaded contours show land ice thickness obtained with 20 mbar of carbon dioxide in response to changes of orbital forcing (WSO (**a**)) and CSO (**b**), warm/cold summer orbit for the northern hemisphere) over the course of two precession cycles (40 kyr of simulation). In the light brown continental areas without ice, the white line is used to represent the old ice sheet extension (WSO case). The Svalbard area is indicated by a red circle. **c,** Ice thickness variation in 10 kyr (WSO case after 20 kyr minus CSO case after 30 kyr of simulation). In **a–c** the continental outline is shown by the 0.5 m elevation contour. **d,** Surface of hemisphere covered by ice (m²) through time ([WSO] and [CSO] indicate which orbital configuration is used).

areas remain ice free as a result of sublimation exceeding snowfall (Supplementary Fig. 10a). Ice volume remains relatively constant for $p_{CO_2}$ = 0.1 to 20 mbar (Supplementary Fig. 10b), owing to an increase in accumulation that compensates for higher ablation rates (Supplementary Fig. 13). In contrast, above 20 mbar, ice extent in the eastern Tropics significantly decreases (Supplementary Fig. 10c). At $p_{CO_2}$ = 100 mbar, most of the continental ice cover disappears, except for remnants over mountain ranges (Supplementary Fig. 10d).

To test the sensitivity of the tropical ice sheets to Milankovitch forcing, experiments with changing orbital parameters were initialized using the steady-state ice sheets for $p_{CO_2}$ = 20 mbar. Although obliquity has been invoked as a possible cause of Neoproterozoic glaciations²⁴, this mechanism remains problematical and cannot account for significant climatic oscillations at low latitudes²⁵,²⁶. We therefore focused on precession as a possible driver, and used two opposite orbital configurations, favouring cold and warm summers, respectively, over the northern tropics (CSO: cold summer orbit and WSO: warm summer orbit; Supplementary Fig. 14). Switching between these configurations causes tropical ice sheets to advance/retreat over several hundred kilometres in 10 kyr (Supplementary Movie 1), with strong asymmetry between hemispheres (Fig. 3). Shifting from WSO to CSO causes ice retreat in the Southern Hemisphere and ice sheet expansion in the Northern Hemisphere (Supplementary Fig. 14c,d). Significant ice volume changes occur between 30°N and S, but are less apparent at higher latitudes. This reflects higher ablation rates in the warmer low latitudes (Supplementary Fig. 14e,h), and higher ice sheet sensitivity to shifting patterns of melt. Larger greenhouse forcing at the end of the Snowball event implies increasing ice sheet sensitivity to subtle insolation changes. Given a strong diurnal cycle²³, our simulations also predict a significant number of days above 0°C in the tropics (Supplementary Fig. 15), consistent with geological evidence for ice rafting, liquid water in lakes and rivers, and photosynthetic microbial communities.

Our results show that geological evidence for glacial–interglacial cycles⁵⁻⁷ is consistent with an enriched Snowball Earth theory. Termination of the Marinoan panglaciation was not a simple switch

from icehouse to greenhouse states, but was characterized by a climate transition during which glacial cycles could be forced by Milankovitch orbital variations. The geochemical evidence presented here implies that at least the upper 60–70% of the Wilsonbreen Formation was deposited in ~10⁵ years, on the assumption that a trend in $p_{CO_2}$ would be evident over longer timescales²¹. Rates of $CO_2$ build-up, however, may have slowed in the later stages of Snowball Earth owing to silicate weathering of exposed land surfaces, so it is possible that the oscillatory phase was more prolonged.

Initiation of low-latitude glaciation in the Neoproterozoic requires low $p_{CO_2}$ (0.1–1 mbar; refs 2,19,20), implying that the oscillatory phase was preceded by a prolonged colder period (~10⁶ to 10⁷ years) during which $p_{CO_2}$ increased gradually as a result of volcanic outgassing²¹. This timescale is in agreement with recent dates indicating the Marinoan lasted ~15 million years²⁷. The basal weathering horizon is consistent with a period of low temperatures and limited hydrologic cycle before the oscillatory phase²,¹⁹.

Further work is needed to refine the upper and lower limits of $p_{CO_2}$ conducive to climate and ice sheet oscillations in Snowball Earth. Factors not included in the present model, such as supraglacial dust or areas of ice-free tropical ocean²⁸⁻³⁰, can be expected to make the Earth system more sensitive to orbital forcing. While many details remain to be investigated, our overall conclusions remain robust.

The Neoproterozoic Snowball Earth was nuanced, varied and rich. We anticipate that detailed studies of the rock record in other parts of the world, in conjunction with numerical modelling studies, will continue to yield insight into the temporal and regional diversity of this pivotal period in Earth history.

## Methods

Methods and any associated references are available in the online version of the paper.

Received 17 February 2015; accepted 8 July 2015; published online 24 August 2015
## References

1. Hoffman, P. F. & Schrag, D. P. The Snowball Earth hypothesis: Testing the limits of global change. *Terra Nova* **14**, 129–155 (2002).

2. Donnadieu, Y., Goddéris, Y. & Le Hir, G. *Treatise on Geochemistry* 2nd edn, Vol. 6, 217–229 (Elsevier, 2014).

3. Hoffman, P. F. Strange bedfellows: Glacial diamictite and cap carbonate from the Marinoan (635 Ma) glaciation in Namibia. *Sedimentology* **58**, 57–119 (2011).

4. Allen, P. A. & Etienne, J. L. Sedimentary challenge to Snowball Earth. *Nature Geosci.* **1**, 817–825 (2008).

5. Rieu, R., Allen, P. A., Plötze, M. & Pettke, T. Climatic cycles during a Neoproterozoic “snowball” glacial epoch. *Geology* **35**, 299–302 (2007).

6. Le Heron, D. P., Busfield, M. E. & Kamona, F. An interglacial on snowball Earth? Dynamic ice behaviour revealed in the Chuos Formation, Namibia. *Sedimentology* **60**, 411–427 (2013).

7. Fairchild, I. J. & Hambrey, M. J. Vendian basin evolution in East Greenland and NE Svalbard. *Precambrian Res.* **73**, 217–233 (1995).

8. Halverson, G. P. in *Neoproterozoic Geobiology and Paleobiology* (eds Xiao, S. & Kaufman, A. J.) 231–271 (Springer, 2006).

9. Li, X.-X., Evans, D. A. & Halverson, G. P. Neoproterozoic glaciations in a revised global palaeogeography from the breakup of Rodinia to the assembly of Gondwanaland. *Sedim. Geol.* **294**, 219–232 (2013).

10. Petronis, M. S. et al. in AGU Fall Meeting 2013 Abstract #GP41A-1107 (American Geophysical Union, 2013).

11. Harland, W. B. *The Geology of Svalbard* (Geological Society, 1997).

12. Creveling, J. R. & Mitrovica, J. X. The sea-level fingerprint of a Snowball Earth deglaciation. *Earth Planet. Sci. Lett.* **399**, 74–85 (2014).

13. Lyons, W. B. et al. The McMurdo Dry Valleys long-term ecological research program: New understanding of the biogeochemistry of the Dry Valley lakes: A review. *Polar Geogr.* **25**, 202–217 (2001).

14. Hoffman, P. et al. Are basal Ediacaran (635 Ma) basal ‘cap dolostones’ diachronous? *Earth Planet. Sci. Lett.* **258**, 114–131 (2007).

15. Fairchild, I. J., Hambrey, M. J., Spiro, B. & Jefferson, T. H. Late Proterozoic glacial carbonates in northeast Spitsbergen: New insights into the carbonate-tillite association. *Geol. Mag.* **126**, 469–490 (1989).

16. Bao, H., Fairchild, I. J., Wynn, P. M. & Spötl, C. Stretching the envelope of past surface environments: Neoproterozoic glacial lakes from Svalbard. *Science* **323**, 119–122 (2009).

17. Cao, X. & Bao, H. Dynamic model constraints on oxygen-17 depletion in atmospheric O₂ after a snowball Earth. *Proc. Natl Acad. Sci. USA* **110**, 14546–14550 (2013).

18. Le Hir, G. et al. The snowball Earth aftermath: Exploring the limits of continental weathering processes. *Earth Planet. Sci. Lett.* **277**, 453–463 (2009).

19. Pierrehumbert, R., Abbot, D. S., Voigt, A. & Koll, D. Climate of the Neoproterozoic. *Ann. Rev. Earth Planet. Sci.* **39**, 417–460 (2011).

20. Pollard, D. & Kasting, J. F. in *The Extreme Proterozoic: Geology, Geochemistry, and Climate* (eds Jenkins, G. S., McMenamin, M. A. S., McKay, C. P. & Sohl, L.) 91–105 (Geophysical Monograph Series 146, American Geophysical Union, 2004).

21. Le Hir, G., Ramstein, G., Donnadieu, Y. & Goddéris, Y. Scenario for the evolution of atmospheric $p_{CO_2}$, during a snowball Earth. *Geology* **36**, 47–50 (2008).

22. Hoffman, P. F. & Li, Z. X. A palaeogeographic context for Neoproterozoic glaciation. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **277**, 158–172 (2009).

23. Pierrehumbert, R. T. Climate dynamics of a hard Snowball Earth. *J. Geophys. Res.* **110**, D01111 (2005).

24. Spiegl, T. C., Paeth, H. & Frimmel, H. E. Evaluating key parameters for the initiation of a Neoproterozoic Snowball Earth with a single Earth System Model of intermediate complexity. *Earth Planet. Sci. Lett.* **415**, 100–110 (2015).

25. Donnadieu, Y., Ramstein, G., Fluteau, F., Besse, J. & Meert, J. Is high obliquity a plausible cause for Neoproterozoic glaciations? *Geophys. Res. Lett.* **29**, 2127 (2002).

26. Paillard, D. Quaternary glaciations: From observations to theories. *Quat. Sci. Rev.* **107**, 11–24 (2015).

27. Rooney, A. D. et al. A Cryogenican chronology: Two long-lasting synchronous Neoproterozoic glaciations. *Geology* **43**, 459–462 (2015).

28. Abbot, D. S. & Pierrehumbert, R. T. Mudball: Surface dust and Snowball Earth deglaciation. *J. Geophys. Res.* **115**, D03104 (2010).

29. Abbot, D. S., Voigt, S. & Koll, D. The Jormungand global climate state and implications for Neoproterozoic glaciations. *J. Geophys. Res.* **116**, D18103 (2011).

30. Rose, B. E. J. Stable “Waterbelt” climates controlled by tropical ocean heat transport: A nonlinear coupled climate mechanism of relevance to Snowball Earth. *J. Geophys. Res.* **120**, 1404–1423 (2015).

## Acknowledgements

This work was supported by the NERC-funded project GR3/NE/H004963/1 Glacial Activity in Neoproterozoic Svalbard (GAINS). Logistical support was provided by the University Centre in Svalbard. This work was granted access to the HPC resources of CCRT under allocation 2014-017013 made by GENCI (Grand Equipement National de Calcul Intensif). We also thank D. Paillard and P. Hoffman for stimulating discussions and valuable insights.

## Author contributions

Field data were collected and analysed by I.J.F., D.I.B., E.J.F., M.J.H., E.A.McM., M.S.P., P.M.W. and C.T.E.S. Geochemical analyses were conducted by H.B. and P.M.W. Model experiments were designed and conducted by G.L.H., Y.D., C.D. and G.R. The manuscript and figures were drafted by D.I.B., I.J.F. and G.L.H., with contributions from the other authors.

## Additional information

Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to D.I.B.

## Competing financial interests

The authors declare no competing financial interests.
**Methods**

**Sedimentology.** Lithofacies were classified based on grain size, internal sedimentary structures and deformation structures, and bounding surfaces. Detailed stratigraphic logs were made in the field, supplemented by drawings and photographs of key features. Samples were taken for polishing and thin sectioning, to allow detailed examination of microstructures in the laboratory. In addition, data were collected on clast lithology, shape, surface features and fabric. Diamictites of the Wilsonbreen Formation are commonly very friable, allowing included clasts to be removed intact from the surrounding matrix, allowing measurement of both clast morphology and orientation, using methods developed for unlithified sediments. Clast morphology (shape, roundness and surface texture) was measured for samples of 50 clasts to determine transport pathways. Clast fabric analysis was performed by measuring a-axis orientations of samples of 50 clasts with a compass-clinometer, and data were summarized using the eigenvalue or orientation tensor method. Oriented samples for measurement of anisotropy of magnetic susceptibility (AMS) were collected using a combination of field-drilling and block sampling. AMS was measured using an AGICO KLY-3 Kappabridge operating at 875 Hz with a 300 A m⁻¹ applied field at the University of Birmingham and an AGICO MFK-1A Kappabridge operating at 976 Hz with a 200 A m⁻¹ applied field at New Mexico Highlands University.

**Geochemistry.** Laboratory procedures for extracting, purifying and measuring the triple oxygen (δ¹⁸O and Δ¹⁷O) and sulphur (δ³⁴S) isotope composition of carbonate associated sulphate (CAS) in bulk carbonates are detailed in ref. 16. Briefly, fresh carbonate-bearing rock chips were crushed into fine grains and powders using mortar and pestle. Rinsing the fines with 18 MΩ water revealed little water-leachable sulphate in any of the Wilsonbreen carbonates. Subsequently, about 10 to 30 g carbonates were slowly digested in 1–3 M HCl solutions. The solution was then centrifuged, filtered through a 0.2 µm filter, and acidified before saturated BaCl₂ droplets were added. BaSO₄ precipitates were collected after >12 h and purified using the DDARP method (see Supplementary Information). The purified BaSO₄ was then analysed for three different isotope parameters: Δ¹⁷O, by converting to O₂ using a CO₂-laser fluorination method; δ¹⁸O, by converting to CO through a Thermal Conversion Elemental Analyzer (TCEA) at 1,450 °C; and δ³⁴S, by converting to SO₂ by combustion in tin capsules in the presence of V₂O₅ through an Elementar Pyrocube elemental analyser at 1,050 °C. The Δ¹⁷O was run in dual-inlet mode, whereas the δ¹⁸O and δ³⁴S were run in continuous-flow mode.

Both the Δ¹⁷O and δ¹⁸O were run on a MAT 253 at Louisiana State University, whereas the δ³⁴S was determined on an Isoprime 100 continuous-flow mass spectrometer at the University of Lancaster, UK. The Δ¹⁷O was calculated as Δ¹⁷O = δ¹⁷O - 0.52 × δ¹⁸O in which δ' ≡ 1,000ln (R/sample/R/standard) and R is the molar ratio of ¹⁸O/¹⁶O or ¹⁷O/¹⁶O. All δ values are in Vienna Standard Mean Ocean Water VSMOW and Vienna Canyon Diablo Troilite (VCDT) for sulphate oxygen and sulphur, respectively. The analytical standard deviation (1σ) for replicate analysis associated with the Δ¹⁷O, δ¹⁸O and δ³⁴S are ±0.05‰, ±0.5‰ and ±0.2‰, respectively. As the CAS is heterogeneous in hand specimens, the standard deviation is for laboratory procedures. δ³⁴S values were corrected against VCDT using within-run analyses of international standard NBS-127 (assuming δ³⁴S values of +21.1‰). Within-run standard replication (1 s.d.) was <0.3‰. All geochemical data are included in Supplementary Table 1.

**Numerical modelling.** Model runs were conducted with a coupled atmospheric general circulation model (LMDz) and ice sheet model (GRISLI: GRenoble Ice Shelf and Land Ice model). LMDz (spatial resolution 4° in latitude × 5° in longitude with 38 vertical levels) was run with prescribed continental ice to climatic equilibrium. GRISLI has a 40 km grid size and is driven with downscaled climatic fields of surface air temperature, precipitation and evaporation. To capture ice sheet–climate feedbacks, LMDz is rerun using the new ice sheet distribution and topography. This procedure was repeated each 10 kyr to investigate orbital forcing.

Surface mass balance (accumulation minus sublimation and melting) was computed from monthly mean temperature, precipitation and evaporation rate. Melt rate is calculated using the positive-degree-day method.

No sea ice dynamics treatment is specified, the sea ice cover is prescribed and a thickness of 10 m is imposed. Ice albedo is fixed at 0.6, whereas snow albedo varies from 0.9 from 0.55 as a function of the zenith and ageing process. Land ice/snow free surface has the characteristic of a bare soil (rocky regolith) with an albedo of 0.3.

**Code availability.** Code for the GCM LMDz can be accessed at: http://lmdz.lmd.jussieu.fr. Code for the ISM GRISLI (GRenoble Ice Shelf and Land Ice model) is not available.

Further details of the methods and modelling procedures are provided in the Supplementary Information in the online version of the paper.
<figure><img src="image_104.png" /></figure>

RESEARCH ARTICLE

10.1002/2016PA002928

**Key Points:**

* Earth system model providing the first detailed simulation of Middle to Late Ordovician land ice growth

* Model/data comparison suggests a Darriwilian age for glacial onset

* A single ice sheet of large extent covered Gondwana during the Hirnantian glacial maximum

**Supporting Information:**

* Supporting Information S1

* Supporting Information S2

* Figure S1

**Correspondence to:**

A. Pohl,
pohl.alexandre@gmail.com

**Citation:**

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

Received 20 JAN 2016

Accepted 21 MAY 2016

Accepted article online 28 MAY 2016

Glacial onset predated Late Ordovician climate cooling

Alexandre Pohl¹, Yannick Donnadieu¹, Guillaume Le Hir², Jean-Baptiste Ladant¹, Christophe Dumas¹, Jorge Alvarez-Solas³, and Thijs R. A. Vandenbroucke⁴,⁵

¹Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France, ²Institut de Physique du Globe de Paris, Université Paris7-Denis Diderot, 1 rue Jussieu, Paris, France, ³Departamento Astrofísica y Ciencias de la Atmósfera, Universidad Complutense de Madrid, Madrid, Spain, ⁴Department of Geology, Ghent University, Krijgslaan, Belgium, ⁵Evo-Eco-Paléo UMR 8198, Université de Lille, Villeneuve d'Ascq CEDEX, France

**Abstract** The Ordovician glaciation represents the acme of one of only three major icehouse periods in Earth's Phanerozoic history and is notorious for setting the scene for one of the "big five" mass extinction events. Nevertheless, the mechanisms that drove ice sheet growth remain poorly understood and the final extent of the ice sheet crudely constrained. Here using an Earth system model with an innovative coupling method between ocean, atmosphere, and land ice accounting for climate and ice sheet feedback processes, we report simulations portraying for the first time the detailed evolution of the Ordovician ice sheet. We show that the emergence of the ice sheet happened in two discrete phases. In a counterintuitive sequence of events, the continental ice sheet appeared suddenly in a warm climate. Only during the second act, and set against a background of decreasing atmospheric CO₂, followed steeply dropping temperatures and extending sea ice. The comparison with abundant sedimentological, geochemical, and micropaleontological data suggests that glacial onset may have occurred as early as the Middle Ordovician Darriwilian, in agreement with recent studies reporting third-order glacioeustatic cycles during the same period. The second step in ice sheet growth, typified by a sudden drop in tropical sea surface temperatures by ~8°C and the further extension of a single, continental-scale ice sheet over Gondwana, marked the onset of the Hirnantian glacial maximum. By suggesting the presence of an ice sheet over Gondwana throughout most of the Middle and Late Ordovician, our models embrace the emerging paradigm of an "early Paleozoic Ice Age."

**1. Introduction**

The Middle to Late Ordovician (470–444 Ma) has long been considered an enigmatic and unique period in Earth's history. Much of this mystery revolved around our inability to explain how a glacial pulse [Finnegan et al., 2011], high pCO₂ levels potentially exceeding 20 times the preindustrial atmospheric level [Berner, 1990; Yapp and Poths, 1996] (PAL, 1 PAL = 280 ppm), a global biodiversity crisis [Sheehan, 2001], and major perturbations of the carbon cycle [Brenchley et al., 1994] could have coincided. Over the last few years, however, evidence has accumulated that some of these apparent conflicts may have resulted from erroneous interpretations of the relatively poorly preserved and incomplete Ordovician sedimentary record [Ghienne et al., 2014]. Recent studies suggest that the atmospheric CO₂ content was probably much lower than had been previously proposed [Rothman, 2002; Herrmann et al., 2004; Berner, 2006; Vandenbroucke et al., 2010a; Pancost et al., 2013] and that ocean temperatures reached modern-like values as early as the Middle Ordovician, distinctly contrasting with earlier visions of a perennial Ordovician supergreenhouse [Trotter et al., 2008; Vandenbroucke et al., 2009]. Despite a more comprehensive analysis of the physical evidence, major uncertainties persist about the mechanisms that drove ice sheet growth and about its final extent. Rooted in the disperse nature of the Late Ordovician glacial sedimentary records in North Africa [Ghienne et al., 2007], South Africa, and South America [Díaz-Martínez and Grahn, 2007], one of the major questions remains whether those glacial outcrops are indicative of a single, large ice sheet that developed from North Africa (then at a polar latitude) to South Africa (at tropical latitudes in latest reconstructions) [Torsvik and Cocks, 2013] or if independent ice sheets of smaller extent coexisted over the South Pole and in the tropics [Ghienne et al., 2007; Le Heron and Dowdeswell, 2009].
It is also still unclear how glacial onset could have possibly occurred in the midst of a prolonged period of hot climate. Evidence for warm shallow seas in the tropics (up to 37°C) [Finnegan et al., 2011], as suggested by most recent Ordovician geochemical data [Trotter et al., 2008; Finnegan et al., 2011], is in apparent contradiction with Middle to Upper Ordovician deposits where glacioeustatic events are identified [Loi et al., 2010; Turner et al., 2012; Elrick et al., 2013; Dabard et al., 2015; Rasmussen et al., 2016] or suggested [Nielsen, 2004] (see discussion in Amberg et al. [2016]). These eustatic cycles support the presence of ice sheets since the Middle Ordovician Darriwilian (467 Ma) rather than a sudden emplacement of a massive ice sheet restricted to the Hirnantian stage (445–444 Ma).

To investigate how these apparently contradictory interpretations can be reconciled, we employ an Earth system model to propose the first realistic simulations of Middle to Late Ordovician ice sheet growth. The last modeling attempts of this kind date back to two papers published by Herrmann et al. [2003, 2004]. Building on the seminal studies of Crowley and Baum [1991, 1995], they provided major insights into the climatic conditions required to trigger Late Ordovician glacial onset. However, they are more than 10 years old, and the experimental setup that is available to us has significantly advanced. Herrmann et al.'s [2003, 2004] models notably did not resolve ice sheet feedbacks on global climate. The simulated ice sheet did not extend beyond 60°S and therefore poorly compared to the Ordovician glacial sedimentary record [e.g., Ghienne et al., 2007]. More recently, Lowry et al. [2014] asynchronously coupled a slab mixed-layer ocean-atmosphere general circulation model (GENESIS) with an ice sheet model to investigate the impact of the changes in continental configuration on the potential for glacial inception throughout the Paleozoic. However, they did not use specific Ordovician boundary conditions in their coupled simulations—in terms of topography, solar luminosity, and CO₂ level in particular—which resulted in significant discrepancies between Ordovician geological data and modeled ice sheet extent and flow. The flat topography led to an unexpected ice sheet geometry with the maximum ice sheet height simulated along the coasts of Gondwana at 440 Ma. The resulting ice sheet flowed from the coast to the continent interior (Lowry et al. [2014]; see their discussion section), at odds with the sedimentary record of the Hirnantian glaciation [e.g., Ghienne et al., 2007]. They ran additional experiments considering the evolution of the solar luminosity value throughout the Paleozoic, but they did not account for the feedbacks of the ice sheet on global climate in these simulations. To capture processes leading to glacial onset and ensure both reliable ice sheet growth mechanisms and robust steady state ice sheet extent, we here use suitable boundary conditions together with an innovative modeling procedure that for the first time in the Ordovician considers ocean-atmosphere-ice sheet interplays as well as a fully dynamic ocean component. Our model accounts for Ordovician orbital variations and for the wide range of suggested CO₂ levels, i.e., 3–24 PAL (see section 2.2.2).

We compare our modeling results with abundant geochemical, sedimentological, and paleontological data, which provides additional insights into both the extent and duration of the glaciation. However, current ice sheet models exhibit severe issues of hysteresis, making land ice too stable and requiring unreasonably high atmospheric CO₂ levels to cause continental-scale ice sheet retreats [Pollard and DeConto, 2005; Horton and Poulsen, 2009]. Therefore, we did not attempt the simulation of the Late Ordovician deglaciation but identify this as a target for the future. We are confident that ongoing advances in modeling of the cryosphere [Ganopolski et al., 2010; Pollard et al., 2015; Gasson et al., 2016] will help resolve this issue in the years to come.

## 2. Methods

### 2.1. Asynchronous Coupling Design

The deep-time Earth system model used throughout the study is described below. Our aim is to simulate ice sheet-climate equilibria under decreasing pCO₂, representative of the well-established Ordovician cooling trend [Trotter et al., 2008]. This requires coupling an ocean-atmosphere general circulation model with an ice sheet model. The most straightforward way to couple an ocean-atmosphere climate model to an ice sheet model is the purely synchronous coupling method, which would account for exchanges between the ocean, atmosphere, and land ice components at each model time step. The long computational time required by an ocean-atmosphere general circulation model (GCM), however, forbids the use of synchronous coupling to solve deep-time issues, which usually require long ice sheet runs lasting hundreds of thousands of model years [e.g., Ladant et al., 2014]. Several alternative methods have therefore been proposed [Pollard, 2010], including asynchronous coupling with a GCM lookup table, which remains the only solution applicable on very long
<figure><img src="image_105.png" /></figure>

**Figure 1.** Ice sheet results obtained at 12 PAL, with our Earth system model accounting for orbital variations and ice sheet feedbacks on climate, throughout the successive asynchronous coupling steps. Red curves represent the evolution of the land ice cover through ice sheet model integration time (200 kyr per iteration). The thin blue line stands for the 30 kyr obliquity forcing used to interpolate between the CSO and the HSO. The associated right y axis gives the Southern Hemisphere mean summer insolation received at 65°S. The ice sheet model (GRISLI) is first run for 200 kyr from its initial, ice-free state (iter 1), by linearly interpolating between the CSO and HSO temperature and precipitation forcing fields according to the 30 kyr Ordovician obliquity period. Then, the largest ice sheet simulated during this first iteration is returned as a boundary condition to the climate models. New climatic forcing fields are simulated, which therefore take the ice-albedo local radiative effect into account. The latter climatic snapshots are added to the "ice cover" dimension of the GCM lookup table, providing new climatic fields to interpolate with. GRISLI is once again run for 200 kyr from its ice-free state (iter 2), and the ice sheet model climatic forcing fields are interpolated as a function of time both between the CSO and the HSO and between the ice-free state and the glaciated state predicted from previous iteration. The largest ice sheet simulated is used to force the third iteration (new climatic snapshots are added to the GCM lookup table), and the loop is repeated until the ice sheet is in equilibrium with global climate (iter 4 to iter 9), i.e., an additional iteration does not induce a much larger ice sheet extent (compare iter 8 and iter 9).

sheet sizes are not prescribed, hence evolving freely, instead of imposing successive ice sheet extents for the interpolation. As a consequence we build the matrix step by step.

time scales [e.g., Ladant et al., 2014]. In this method, GCM snapshots are gathered into a multidimensional matrix. Each dimension of the matrix corresponds to a given external forcing (e.g., the orbital configuration or the $pCO_2$) and contains surface temperature and precipitation fields simulated for several values of this forcing, which provide the entire information required during the subsequent ice sheet model integration. Then, at any point in the long-term ice sheet model run, the surface mass balance of the ice sheet is calculated based on mean monthly temperature and precipitation fields interpolated between the appropriate cells of the matrix [Pollard, 2010]. In the present case study, the GCM snapshots are collected in a matrix with three dimensions, i.e., the Earth's orbital parameters, the ice sheet size, and the $pCO_2$, following the methodology successfully applied to the Antarctic ice sheet onset at the Eocene-Oligocene transition [Ladant et al., 2014]. However, unlike the Antarctic ice sheet, which has limited spatial extent, the maximum ice sheet extent over the vast Gondwana paleocontinent is virtually unlimited, and as such, there is no reason to preconceive any maximum ice sheet extent prior to the buildup of the 3-D matrix. Here the ice

1. Climatic snapshots representative of ice-free conditions are first obtained: for a given $pCO_2$ (first dimension of the matrix), two climate states are computed by the GCM with an ice-free bedrock topography (second dimension), for two opposed orbital configurations (third dimension) yielding either a cold or a warm summer over the southern high latitudes (cold and hot summer orbits—CSO and HSO, respectively—see section 2.2.2).

2. A first ice sheet is modeled: the ice sheet model is run, for the same $pCO_2$ level, using a sinusoidal forcing in obliquity weighted between the two opposite climates from step 1 (iter 1 in Figure 1). The period of the sinusoidal variation in obliquity is set to 30 kyr, based on spectral analysis of Ordovician halite deposits [Williams, 1991]. Orbital variations are known as a major driver of Pleistocene glacial cycles [e.g., Paillard, 1998], and they have been shown to significantly impact the growth of the ice sheet during the Ordovician as well [Herrmann et al., 2003; Ghienne et al., 2014; Dabard et al., 2015].

3. An ice sheet of larger extent is simulated, provided that climate is favorable to ice sheet growth: in order to account for the feedbacks of the ice sheet on climate, the largest ice sheet simulated previously is reintroduced as a boundary condition into the climate model (still keeping the $CO_2$ constant). Two additional glacial climatic states are simulated, one for each of the two opposite orbital configurations (CSO and HSO). The climatic snapshots are gathered into the matrix of climatic states, each snapshot corresponding to a specific land ice mask and orbital configuration (i.e., a CSO snapshot and a HSO snapshot for each ice sheet extent, including the ice-free state). The ice sheet model is then run, starting from an ice-free state, forced in temperature and precipitation using climatic fields interpolated, throughout model integration time,
between the climatic snapshots of the matrix as a function of both the orbital forcing and the extent of the ice sheet (iter 2 in Figure 1).

4. Step 3 is repeated until steady state is reached between land ice and global climate: a new member is thus added to the "ice sheet extent" dimension of the matrix as many times as necessary, and the resulting ice sheet is neither limited to a given geometry nor artificially pulled toward a maximal extent, which might have been the case with prescribed ice sheets (iter 3 to iter 9 in Figure 1).

5. The whole procedure is repeated for the other CO₂ values.

Climatic snapshots are obtained using successively the ocean-atmosphere general circulation model FOAM (as a sea surface temperature generator) and the state-of-the-art atmospheric model LMDZ, with the same boundary conditions. The model FOAM is first integrated until deep-ocean equilibrium is reached (≥2000 years). Given that the atmospheric fields of temperature and precipitation simulated with FOAM are not sufficiently well resolved to directly force the ice sheet model GRISLI, the sea surface temperatures simulated with FOAM are subsequently used to force the atmospheric model LMDZ, the outputs of which are in turn used to force GRISLI. The alternate use of FOAM and LMDZ has been shown to perform well, as testified by previously published paleoclimate studies [Ladant et al., 2014; Licht et al., 2014].

## 2.2. Model Setup

### 2.2.1. Model Details

FOAM (Fast Ocean Atmosphere Model) version 1.5 is a fully coupled mixed-resolution ocean-atmosphere general circulation model with no flux corrections [Jacob, 1997]. Its quick turnaround time allows for long millennium-scale integrations. FOAM is consequently well designed for paleoclimate studies, for which purpose it has been routinely used in the past, including for the Ordovician [Nardin et al., 2011; Pohl et al., 2014, 2015]. The atmospheric component of the model is the National Center for Atmospheric Research's (NCAR) Community Climate Model version 2 (CCM2) benefitting from the upgraded radiative and hydrologic physics from CCM3 version 3.2. It is run at a R15 spectral resolution (4.5° × 7.5°) with 18 vertical levels (14 in the troposphere). The ocean module is the Ocean Model version 3 (OM3). It is a finite difference, 24-level z coordinate ocean general circulation model providing a higher resolution than the atmospheric module (1.4° × 2.8°) on a rectangular grid. Sea ice is simulated using the thermodynamics of the NCAR's Climate System Model 1.4 sea ice module, which is based on the Semtner [1976] three-layer thermodynamic snow/ice model.

The LMDZ (Laboratoire de Météorologie Dynamique Zoom) model is an IPCC-class (Intergovernmental Panel on Climate Change) atmospheric general circulation model benefitting, as the atmospheric component of the state-of-the-art IPSL (Institut Pierre-Simon Laplace) Earth system model, from the latest physical and dynamical refinements [Hourdin et al., 2013]. A midresolution version of the model is used in this study, providing a resolution of 1.9° × 3.75° with 39 vertical levels (21 in the troposphere) on a rectangular grid.

GRISLI (GRenoble Ice-Shelf and Land Ice) is a three-dimensional ice sheet model accounting for thermodynamical coupling between ice velocities and temperatures [Ritz et al., 2001]. GRISLI runs at 40 km × 40 km resolution and simulates both grounded and floating ice, by considering three types of ice flow: ice sheets, ice streams, and ice shelves. At each time step, the model determines the areas where each type of ice flow applies. Ice shelves are large floating ice plates fed by inland ice. Their location is determined based on a flotation criterion, and ice velocities in these grid points are computed using the shallow-shelf approximation with a basal stress set to zero. The FOAM model does not have the resolution needed to simulate oceanic melting under the ice shelves. This requires specific experimental setups, involving high-resolution (regional) climate models accounting for ocean temperature and circulation patterns under the shelves [e.g., Heimbach and Losch, 2012]. In this study, we therefore used a parameterization considering a basic dependence of the melt rate on the ocean bathymetry. We imposed a melt rate of 0.2 m yr⁻¹ on the continental shelf (water depth ≤600 m) and 2 m yr⁻¹ in open ocean conditions (water depth >600 m). These values are well within the range defined in previous studies on Antarctica [Pollard and DeConto, 2012, their Figure 4b]. Calving occurs when the thickness of the ice shelf drops below 150 m. All areas that are not ice shelves are grounded ice, including ice sheets and ice stream zones. Ice sheets are regions of grounded ice where ice flow is mainly due to deformation, basal sliding occurring only where the base is at the pressure melting point. Inland ice velocity fields are computed using the shallow-ice approximation. Because the latter is inaccurate to simulate fast-flowing ice [Bueler and Brown, 2009], such regions—the ice stream zones—are specifically treated by the model. Ice streams are corridors of fast flow with widths of a few kilometers, usually smaller than the
model grid size. It is not possible to simulate them individually. Their effect is parameterized in the model by treating areas which have the large-scale characteristics of ice streams independently. Ice stream zones are defined based on both the local slope of the terrain and the effective water pressure below the ice sheet, the latter as a function of the subglacial water depth. Ice flow in these ice stream zones is computed using both the shallow-ice and the shallow-shelf approximations. Contrary to what is prescribed for ice shelves, basal stress is not set to zero. GRISLI simulates the evolution of ice sheets in response to climatic forcing. The model calculates the surface mass balance of the ice sheet as the sum of accumulation minus ablation, using mean monthly surface temperature and precipitation. Ablation is in turn computed with the positive degree day (PDD) method, using a surface melt refreezing rate of 60% and a standard deviation (σ) of 5°C to represent the air temperature variability over the month [Reeh, 1991]. The PDD model is based on an empirical relation between air temperature and melting rate, assuming that the near-surface temperature transmits the majority of the climatic forcing to the ice sheet. Despite its simplicity, this ablation scheme has been shown to perform well [e.g., Ladant et al., 2014], although some recent studies suggest that explicitly accounting for the melt contribution from insolation may significantly impact the calculation of surface melt [e.g., Robinson and Goelzer, 2014]. GRISLI also accounts for the isostatic adjustment of the bedrock in response to ice load. The common "elastic lithosphere/relaxing asthenosphere" (ELRA) model is used. In this parameterization, the bed underneath the ice sheet is described as a thin, elastic lithosphere governing the spatial shape of the deflection by allowing a regional response to the ice load. The underlying, viscous asthenosphere governs the time-dependent characteristics of the deformation through a standard relaxation time of 3 kyr. GRISLI was initially developed for and validated over Antarctica by Ritz et al. [2001]. It has since then been successfully applied to study the mechanisms that drove Heinrich events (last 80 ka) [Alvarez-Solas et al., 2011, 2013], the inception of the Fennoscandian ice sheet during the last glacial period (90 ka) [Peyaud et al., 2007], the partial melt of the Greenland ice sheet during the last interglacial period (130–115 ka) [Quiquet et al., 2013], and the continental-scale initiation of the Antarctic ice sheet at the Eocene-Oligocene boundary (34 Ma) [Ladant et al., 2014].

### 2.2.2. Boundary Conditions

We use the Late Ordovician (450 Ma) continental configuration from Blakey [2016]. Given that continental vegetation was restricted to nonvascular plants during the Ordovician [Steemans et al., 2009; Rubinstein et al., 2010]—the coverage of which is difficult to estimate—we follow previous studies about Ordovician climate [Nardin et al., 2011; Pohl et al., 2014, 2015] in imposing a bare soil (rocky desert) on every continental grid point. Five realistic topographic classes are used, ranging from the deep ocean (–4000 m) to the Caledonian orogen (3000 m) [Pohl et al., 2014]. Simulations are conducted for two end-member orbital configurations (CSO and HSO), respectively, extremely favorable and unfavorable to south polar ice sheet growth. The HSO (CSO) is defined with an obliquity of 24.5° (22.5°), an eccentricity of 0.05 (0.05), and a longitude of perihelion of 90° (270°) [DeConto and Pollard, 2003]. Following the models of solar physics [Gough, 1981] and previous studies about Ordovician climate [Pohl et al., 2014, 2015], the solar constant is decreased by 3.5% compared to its present value. The values for Ordovician $p\text{CO}_2$ remain a major uncertainty. They were necessarily higher than today, compensating for a weaker insolation, but climatic proxies [Yapp and Poths, 1992; Rothman, 2002; Tobin and Bergström, 2002; Tobin et al., 2005; Vandenbroucke et al., 2010a; Pancost et al., 2013], climate models [Herrmann et al., 2003, 2004], and long-term carbon cycle models [Berner, 1990; François et al., 2005; Berner, 2006; Young et al., 2009; Nardin et al., 2011; Goddéris et al., 2014] give scattered values ranging from less than 5 PAL to more than 20 PAL. In order to cover this wide spectrum, experiments are conducted for $p\text{CO}_2$ values between 3 PAL and 24 PAL.

## 2.3. Sensitivity Tests

### 2.3.1. One-Way Coupling Versus Asynchronous Coupling Method

Most experiments discussed in this study are conducted with the Earth system model described previously (i.e., the asynchronous coupling method; see section 2.1). Additional experiments are conducted with a simpler experimental design to allow comparison with previously published work and to test the sensitivity of our results to various parameters. Similar to previous studies [Herrmann et al., 2003, 2004], we unidirectionally force the ice sheet model (GRISLI) with atmospheric outputs representative of ice-free conditions. As in the asynchronous coupling method, GRISLI updates its surface temperature and precipitation fields as a function of ice sheet height, through an adiabatic lapse rate for the temperature (5 K km$^{-1}$) and an exponential law for the precipitation (Clausius Clapeyron dependency), thus capturing the effect of the ice sheet on local climate. The one-way forcing method, however, differs from our asynchronously coupled model by lacking
ice sheet feedbacks on regional and global climates. The ice sheet model is not able to account, by itself, i.e., without being provided with additional climatic forcing fields, for larger-scale effects such as the impact of the ice sheet on the atmospheric circulation. It also significantly underestimates the regional cooling induced by the ice sheet by neglecting the strong ice-albedo positive feedback. For these one-way simulations, climatic forcing fields are computed using the method described previously (FOAM and then LMDZ). We use the cold summer orbital configuration (CSO), which is particularly favorable to Southern Hemisphere ice sheet growth. The resulting steady state ice sheets therefore represent the maximum land ice extent and volume that can be reached without considering land ice feedbacks on global climate.

### 2.3.2. Response of Ice Dynamics to Continental Topography

Published continental reconstructions do not necessarily concur with each other regarding the location of Ordovician highlands over Gondwana [e.g., Blakey, 2016; Scotese, 2016], and none of them provide numerical values for the altitude of these mountainous areas. In our baseline runs, we define mountainous areas based on Blakey's [2016] Late Ordovician reconstruction (section 2.2.2), and we assign constant altitude of 2000 m to these highlands. To estimate the dependence of our modeling results on the large uncertainties in Ordovician topographic reconstructions, we run two additional sets of experiments involving a flat Gondwana super-continent over the South Pole (constant altitude of 300 m) and a lowered topography where mountainous areas are set to 600 m altitude. For each of these tests, the one-way forcing method is used. LMDZ is run for each pCO₂ in a cold summer orbital configuration, using the appropriate topography, but keeping the sea surface temperatures previously simulated in our one-way baseline runs (section 2.3.1). GRISLI is subsequently forced using the LMDZ climatic fields. Using fixed sea surface temperatures allows to specifically investigate the impact of mountainous areas on the atmospheric circulation and thus on the patterns of the precipitation and temperature forcing fields provided to GRISLI, by avoiding ocean dynamics complex responses to changing topography.

### 2.3.3. Response of Ice Dynamics to Glacier/Bedrock Coupling

In spite of their limited spatial extent, ice streams today represent the most dynamic components of ice sheets, notably accounting for more than 90% of all the ice and sediment discharged by the Antarctic ice sheet [Bennett, 2003]. Ice streams consequently constitute a central component of the dynamic behavior of present ice sheets. In the model GRISLI, they are simulated using the shallow-shelf approximation (see section 2.2.1). In this context, the driving term for the ice flow computation is the basal sliding, which is computed as a function of basal shear stress and effective pressure, the latter as a function of the subglacial water pressure. Because shear stress and effective pressure are highly dependent on the geology of the subglacial environment, a sediment map is usually provided as a boundary condition to the model. However, the nature of the Ordovician continental surface is unknown over virtually the entire Gondwana supercontinent, challenging the accurate reconstruction of sediment thickness. The strength of the glacier/bed coupling therefore remains a major uncertainty in our deep-time modeling study, potentially biasing the modeled land ice volume and extent. In order to quantify the amplitude of this possible bias, we conduct a sensitivity test of the glacier/bed coupling strength by varying the basal dragging in the model. The basal dragging is the force generated by the friction between the flowing ice sheet and its bedrock, which hampers ice motion. In the model, the direction of basal dragging is opposite to that of basal velocity [Ritz et al., 2001], and its value is a function of the latter [Alvarez-Solas et al., 2011]:

$$ \tau_b = -v^2 N U_b \quad (1) $$

where $τ_b$ is the basal dragging, $U_b$ is the basal velocity, N represents the effective pressure, and $v^2$ is an empirical parameter representing the basal friction coefficient in ice streams, the value of which will be varied. For these tests, the climatic forcing fields used in GRISLI are the same in all the experiments: we use the one-way forcing method, the cold summer orbit, and Blakey's topography [2016] with mountainous areas set to 2000 m.

## 3. Results

### 3.1. One-Way Experiments

To facilitate the comparison of our results with previous modeling studies [Herrmann et al., 2003, 2004; Lowry et al., 2014], we first conduct experiments that account neither for ice sheet feedback processes nor for orbital variations (see section 2.3.1). As expected, ice sheets of various extent and height are simulated, depending on the pCO₂ (Figure 2). The lower the pCO₂, the larger the ice sheet. The relationship between pCO₂ and land
<figure><img src="image_106.png" /></figure>

**Figure 2.** South polar projection maps of Late Ordovician steady state ice sheet heights (m) simulated at various CO₂ levels with the “one-way” forcing method, i.e., with no land ice feedbacks on global climate. For these simulations, the cold summer orbit is used. The thick white line corresponds to a mean annual sea ice fraction equal to 50%. Background shading is topography/bathymetry. Latitude is shown at 30° intervals.

ice volume is linear in a log-axis graph (Figure 3a, blue line). Herrmann et al. [2004] employed a similar experimental setup to run simulations on two Late Ordovician continental configurations from Scotese and McKerrow [1991]. They investigated the impact of several factors on glacial onset, i.e., pCO₂, sea level, and heat transport. Their set of experiments including a low sea level and a “normal” (versus “decreased”) oceanic heat transport was typified by boundary conditions very close to the ones used in this study. As summarized in their Table 2, their simulation conducted at 15 PAL did not lead to glacial onset, whereas a continental-scale ice sheet was simulated at 8 PAL. Figure 2 reveals that in our models, (i) an ice sheet of limited spatial extent develops on polar highlands at 16 PAL, apparently indicating a glacial onset at higher CO₂ levels than Herrmann et al. [2004]. Herrmann et al. [2004], however, prescribed a flat continent, with altitudes set to 250 m on coastal grid points and 500 m everywhere else, compared to 2000 m for ancient orogens and 300 m elsewhere in our case. Interestingly, additional simulations with mountainous areas lowered to 600 m (see section 3.3.1) are typified by a glacial onset occurring at 8 PAL in our models, eventually providing a CO₂ threshold for glacial onset very comparable to the one proposed by Herrmann et al. [2004]; (ii) The ice sheet simulated at 8 PAL extends to ~60°S (Figure 2), which is in reasonable agreement with the extent simulated by Herrmann et al. [2004] for the same pCO₂ on their latest Ordovician reconstruction (see their Figure 6).

The experimental setup used by Lowry et al. [2014] was very comparable to the one of Herrmann et al. [2004]. Similar to Herrmann et al. [2004], climatic forcing fields were simulated using the GENESIS model on a flat topography under a cold Southern Hemisphere summer orbital configuration. In their simulations using a solar constant set to -2.5% below present, Lowry et al. [2014] obtained a threshold CO₂ level for glacial inception below 3 PAL. Further accounting for early Paleozoic solar luminosity values raised the pCO₂ threshold to 4 PAL. Land ice growth did not occur at 8 PAL. This discrepancy between the results of Herrmann et al. [2004] and ours, on the one hand, and the ones of Lowry et al. [2014], on the other hand, may at least partly result from the choice of the vegetation cover. Here we use a barren soil, the albedo of which corresponds to a rocky desert (section 2.2.2). Herrmann et al. [2004] used no vegetation and intermediate soil color values. Lowry et al. [2014], on the contrary, coupled GENESIS to BIOME4, an equilibrium terrestrial vegetation model simulating the changing distribution of 28 modern biome types in response to climatic forcing. Sensitivity tests demonstrated that ice sheet initiation was sensitive to vegetation (supporting information of Lowry et al. [2014]).
<figure><img src="image_107.png" /></figure>

**Figure 3.** Ice sheet results at various CO₂ levels. (a) Ice sheet extent simulated under a cold summer orbit without taking ice sheet feedbacks on climate into account (one way, blue line) and with our asynchronous coupling procedure accounting for orbital variations and ice sheet feedbacks (red line). Global mean annual surface air temperature (black line, right y axis) is displayed for the asynchronous coupling method only, once the ice sheet obtained steady state with global climate. (b–g) South Pole projection maps of Late Ordovician equilibrium ice sheet heights (m) simulated at each CO₂ level, using the asynchronous coupling method accounting for orbital variations and ice sheet feedbacks. Background shading is topography/bathymetry. The thick white line corresponds to a mean annual sea ice fraction equal to 50%. Latitude is given at 30° intervals. See Figure 8 for a time series of land ice growth. The grey rectangle in Figure 3a indicates the domain where ice sheets sustain in an otherwise relatively warm climate; see section 4.2.
The use of a bare ground induced larger ice volumes and surface area. In some time slices of the Phanerozoic, a barren soil raised the CO₂ threshold for glaciation.

It is more difficult to directly compare our results with the work from Herrmann et al. [2003], who conducted their experiments for a narrower pCO₂ range (8–12 PAL) and focused on the sensitivity of the modeled ice sheet volume to the obliquity forcing. Orbital variations are not considered in our one-way experiments, though they are in our asynchronously coupled runs (see section 2.1).

It appears that using the one-way forcing method, our models compare relatively well with the results of Herrmann et al. [2004]. In the simulations of Lowry et al. [2014], glacial onset occurs for slightly lower CO₂ levels, possibly due to their modern vegetation cover.

## 3.2. Accounting for Ice Sheet Feedback Processes

Simulations are run with the same models and boundary conditions, using the asynchronous coupling method accounting for orbital variations and, for the first time, for the feedbacks of the ice sheet on global climate. Changes in ice surface elevation and extent affect surface air temperatures and atmospheric circulation and vice versa (see section 2.1). Figure 3a displays the equilibrium surface air temperature (black line) and ice sheet extent (red line) simulated at various pCO₂ levels. At 24 PAL, the climate is warm (global, mean annual surface air temperature of ~25.7°C, Figure 3a) and no land ice is observed over the South Pole (Figure 3b). When CO₂ is decreased to 16 PAL, the climate cools slightly (by 1.8°C) and small ice caps appear on the polar highlands (Figure 3c). When CO₂ is lowered to 12 PAL, an ice sheet suddenly expands to the midlatitudes (~45°S, Figure 3d) and global temperature decreases by 3.9°C (from 23.9 at 16 PAL to 20.0°C at 12 PAL, Figure 3a). After this first step of sudden ice sheet growth, the position of the ice front stabilizes and ice sheets of similar size and volume are simulated at 12, 10, and 8 PAL (Figures 3a, 3d–3f). When CO₂ is further lowered to 3 PAL, global temperature suddenly plummets by ~14°C (Figure 3a) inducing a second-step latitudinal expansion of the land ice by a further 15° (Figures 3a and 3g), whereby the continental ice sheet reaches the tropical realm (30°S). Sea ice—previously confined to higher paleolatitudes—expands to 45°S (Figure 3g).

During the first step, global climate suddenly shifts from a greenhouse (pCO₂> 12 PAL) to an icehouse mode with the land ice front resting at the midlatitudes (12–8 PAL, Figures 3d–3f). Despite the continental configuration favoring a strong continentality over the southern high latitudes, analysis of the mechanisms of incipient glaciation reveals that accumulation is never the limiting factor for land ice growth (Figure 4). The latter is, on the contrary, driven by lack of ablation (Figure 4), as reported for the Laurentide ice sheet during the Last Glacial Maximum [Bonelli et al., 2009]. The Ordovician ice sheet nucleates near the South Pole on highlands where freezing persists during the austral summer season, i.e., where ablation is virtually zero throughout the year. Once the ice cover is sufficiently large, the ice-albedo and height mass balance feedbacks amplify the extension of the zero-ablation zone (Figure 4). The ice sheet expands equatorward until it reaches the midlatitudes, where it stabilizes from 12 to 8 PAL (Figures 3d–3f). This apparent contradiction between pCO₂ decrease and unchanged land ice extent, between 12 PAL and 8 PAL (grey domain in Figure 3a), is caused by higher ablation rates in the warmer midlatitudes (supporting information), limiting the ice line close to 45°S. The nonlinear extension of the land ice simulated in this first step is comparable, to some extent, to the Pleistocene glacial cycles, during which interstadials (when only Greenland was covered by ice) were succeeded by stadials when the Laurentide ice sheet suddenly extended to 40°N [Ehlers and Gibbard, 2007]. The ice sheet extension of the second step, simulated at 3 PAL, relies on a specific mechanism. Pohl et al. [2014, 2015] showed that there is a threshold CO₂ value, below which Ordovician climate abruptly shifts from a warm climatic equilibrium with no sea ice in the Northern Hemisphere to a much colder state with sea ice extending down to the midlatitudes. They attributed this climatic instability to the particular ocean dynamics developing in the mainly oceanic Ordovician Northern Hemisphere. At 3 PAL, the sudden spread of sea ice to the midlatitudes (Figure 3g), associated with the crossing of this climatic tipping point, induces a global climate cooling by 14°C (Figure 3a). This abrupt cooling leads in turn to the advance of the continental ice front to the tropical latitudes (Figure 3g).

The quasi-absence of ice shelves in our model runs (Figure 3) results from the very strong ablation outside the grounded ice sheet on the one hand (Figure 4h) and from the poorly constrained Ordovician bathymetry on the other hand. The South Pole, in particular, is the single area typified by lower ablation rates (Figure 4h). There, however, open ocean conditions suddenly succeed to the Gondwana peneplain continent in our
<figure><img src="image_108.png" /></figure>

**Figure 4.** Respective contributions of accumulation and ablation to the growth of the Ordovician ice sheet, illustrated by three snapshots taken from the simulation conducted with ice sheet feedbacks and orbital variations at 12 PAL, at (a–c) 22 ka, (d–f) 44 ka, and (g–i) 173 ka. Latitude is shown at 15° intervals, and land ice cover is systematically displayed through ice sheet height isolines. (j) The three snapshots (22, 44, and 173 ka) are placed on the time series of ice sheet growth (similar to iter 9 in Figure 1), where the red curve represents the increasing extent of the ice sheet through the ice sheet model integration time (200 ka) and the blue curve the oscillations of the orbital forcing between the CSO and HSO end-members (see section 2.2.2). Comparing the three plots provided for each selected model year reveals that accumulation is never the limiting factor for ice growth, because the positive accumulation zone always exceeds the ice sheet extent (Figures 4a, 4d, and 4g). On the contrary, potential ablation is very strong outside of the ice sheet edges, thus limiting further ice sheet extension. The enlargement of the zero-ablation zone (Figures 4b, 4e, and 4h)—driven by local cooling due to the ice-albedo and height mass balance positive feedbacks—allows the ice sheet to develop through the elapsing simulation (Figures 4c, 4f, and 4i).

bathymetric reconstruction, with no fine-scale slope in between (see Figure 4c), which implies strong oceanic melt rates under the shelves in the model (see section 2.2.1). Moreover, there are very few embayments to favor the development of ice shelves. Given that floating ice does not contribute to sea level variations when it forms or breaks up, uncertainties regarding the extent of simulated ice shelves do not alter our estimates of the glacioeustatic fall discussed in section 4.4.

### 3.3. Sensitivity Analysis

We here quantify the impact of the basal dragging and of the bedrock topography on the geometry of the simulated ice sheet. Because running our asynchronously coupled model for each of these additional experiments would be time consuming and expensive, we here use the one-way forcing method and we conduct simple, first-order sensitivity tests under the cold summer orbital configuration. We previously demonstrated that using the one-way method, a very low CO₂ level (3 PAL) is required to simulate a continental-scale ice sheet similar to the ice sheet modeled with our more sophisticated setup (compare Figure 2f with Figures 3d–3f). We therefore conduct most of our sensitivity tests at 3 PAL, in order to ensure that the results of these tests are instructive regarding the geometry of such a continental-scale ice sheet.

#### 3.3.1. Sensitivity to Topography

The most straightforward effect of a mountainous area is to allow the continental surface to intercept the 0°C isotherm, where a flat topography would not, thus making the nucleation of ice centers easier. Conversely, considering lowered relief is expected to shift the threshold for glacial onset toward lower pCO₂ values. With the default topography used in our baseline runs, glacial onset takes place at 12 PAL (Figure 2). When mountains are lowered to 600 m, the threshold pCO₂ value for glaciation is effectively shifted to 8 PAL. With a flat topography, land ice nucleates only at 3 PAL over Gondwana.
<figure><img src="image_109.png" /><figcaption>**Figure 5.** Sensitivity of simulated ice sheet height and extent to topography. Each row corresponds to an experiment. (a, c, and e) The topography used in the model. (b, d, and f) The ice sheet height simulated at equilibrium together with simulated ice sheet volume (V, $1 \text{ Mkm}^3 = 10^6 \text{ km}^3$) and extent (E, $1 \text{ Mkm}^2 = 10^6 \text{ km}^2$). For these simulations, we use the one-way forcing method and the cold summer orbital configuration.</figcaption></figure>

Figure 5 compares equilibrium ice sheets simulated at 3 PAL, using each topography. Although topography significantly impacts the ice sheet geometry, equilibrium ice sheet extent is almost insensitive to topography. Lowered mountains (Figures 5c and 5d), or even a flat Gondwana supercontinent (Figures 5e and 5f), alter the modeled land ice area by no more than 4.5%. The maximum latitudinal extent of the ice sheet, in particular, is remarkably stable. The regional advance of the ice front toward lower latitudes in the center of Gondwana, when topography is flattened (compare Figure 5b with Figures 5d and 5f), results from the changing atmospheric circulation when highlands are removed. Zonal winds are not hampered by relief any more, and they
<figure><img src="image_110.png" /></figure>

**Figure 6.** Austral summer (December, January, and February average) surface air temperature simulated using the one-way forcing method under CSO, land ice-free conditions using (a) the default topography reconstructed from Blakey (2000 m) and (b) the flat topography. Continental mass outlines are shown with the black, stepped, thin black lines. In Figure 6a, main mountains typifying the topography used over Gondwana during the simulation are shown with red, stepped lines. In Figure 6b, the topography is set to 300 m everywhere over Gondwana during the simulation but the location of former relief is mapped in order to facilitate the comparison with Figure 6a. In both cases, the thick black line stands for the 0°C isotherm.

allow the austral summer 0°C isotherm to advance toward the tropics (Figure 6). There is, however, quasi-cancelation between this regional extension of the ice sheet and the retreat of the ice sheet from surrounding mountainous areas due to the warming effect of the lapse rate for lowering topography (Figures 5b, 5d, and 5f). As a result, the extent of the ice sheet only slightly increases (+1 Mkm², Figures 5b, 5d, and 5f) and the land ice volume follows (+24 Mkm³, Figures 5b, 5d, and 5f).

### 3.3.2. Sensitivity to Basal Dragging

Land ice thickness simulated in our baseline runs reaches 6000 m over extensive areas, and the ice sheet systematically displays very steep flanks (Figure 7a). These characteristics are comparable with the results from Herrmann et al. [2004] (see their Figure 6), but they are not realistic when compared to present-day Antarctic [e.g., Fretwell et al., 2013] or LGM Laurentide [e.g., Heinemann et al., 2014] ice sheets. Furthermore, the model simulates virtually no ice stream zones over the entire Gondwana supercontinent (Figure 7b), which is at odds with the sedimentary record of the Ordovician glaciation [e.g., Denis et al., 2007; Ghienne et al., 2007]. These discrepancies arise from weakly dynamic ice sheet, notably resulting from the absence of topographic undulations in our poorly constrained Ordovician topographic reconstruction. The latter includes topographic classes with constant elevation, with no valleys to promote fast-flowing ice streams (see section 2.2.1). In order to quantify this bias, we conduct additional experiments using the basal dragging parameterization developed by Alvarez-Solas et al. [2011] to account for the ice sheet instability. This basal dragging code was successfully applied to study the Heinrich event 1 [Alvarez-Solas et al., 2011]. It uses, in the ice stream zones, a basal dragging coefficient $v^2$ (see section 2.3.3) that is lower than the one used in our baseline runs. It is also less restrictive regarding the conditions required to generate ice stream zones. In contrast to the basal dragging parameterization used in our baseline runs, it uses a sediment map and allows ice streams to develop over areas with very limited basal water content if the bedrock is covered by sediments. For the Ordovician, no detailed lithologic maps of Gondwana are available, but some areas were undoubtedly covered by sediments [Denis et al., 2007]. In the absence of better constraints, we here impose a uniform sediment cover over the supercontinent, and we set the basal friction coefficient for ice streams to $v^2 = 2 \times 10^{-5}$ (dimensionless), following Alvarez-Solas et al. [2011]. Results obtained with this setup are displayed in Figures 7c and 7d.

Ice sheet geometry simulated at 3 PAL with the Alvarez-Solas et al.'s [2011] parameterization shows large differences with our baseline simulation (Figures 7a and 7c). The equilibrium ice sheet extent does not change, but the ice sheet flattens and simulated ice volume decreases from 171 Mkm² to 99 Mkm², eventually representing only 58% of the volume simulated in the baseline experiment (Figures 7a and 7c). The constant ice sheet extent, on the one hand, is due to the unchanged climatic forcing. Decreased ice volume, on the other hand, results from high land ice velocities over extensive ice stream zones (Figures 7b and 7d). Figure 7e (black points) shows the results from additional experiments with various values of the basal dragging coefficient $v^2$. It confirms that the strength of the glacier/bedrock coupling impacts the simulated ice volume but not the ice sheet extent. These results are in agreement with the studies of Clark and Pollard [1998] and Tabor and Poulsen [2016], who simulated similar Laurentide ice sheet extents, but different volumes, as sediment cover changed across the middle Pleistocene transition (~0.9 Ma). We also conduct an experiment featuring a
<figure><img src="image_111.png" /></figure>

**Figure 7.** Sensitivity of simulated land ice thickness and extent to basal dragging coefficient. These simulations are conducted using the one-way forcing method and a cold summer orbit. (a, b) Results from our baseline experiments. (c, d) Results obtained using the basal dragging parameterization developed by *Alvarez-Solas et al.* [2011] with the standard basal dragging coefficient: $v^2 = 2 \times 10^{-5}$. Isolines correspond to ice sheet height. They are displayed every 500 m. (e) The dependence of simulated ice sheet extent (circles) and volume (squares, right y axis) on the basal dragging coefficient, using the parameterization from *Alvarez-Solas et al.* [2011]. Red circles and squares show the values obtained with the standard basal dragging coefficient of $2 \times 10^{-5}$, which results are displayed in Figures 7c and 7d. Grey circles and squares stand for the simulation using the same coefficient and a sediment map in addition. Black symbols stand for simulations where the value of the basal dragging coefficient is varied (x axis), with no sediment map.
<figure><img src="image_112.png" /><figcaption>**Figure 8.** Ice sheet height simulated at 12 PAL, under varying orbital forcing, throughout the successive asynchronous coupling steps (south polar projection maps). (a–i) The largest ice sheet simulated during each one of the successive iterations, which is subsequently returned as a boundary condition in the climate models to simulate the next ice sheet extent (see section 2.1). Background shading is topography/bathymetry. Latitude is shown at 30° intervals.</figcaption></figure>

partial sediment cover. This is a first-order test with sediments present everywhere in topographic lows (i.e., a penplain continent) but absent where the elevation of the terrain is above 300 m. Highlands are considered as rocky bedrock, and they are characterized by the very strong glacier/bedrock coupling typifying inland ice in GRISLI. This spatialization of the basal dragging strength induces differential sliding over Gondwana. Results (Figure 7e, grey points) show that including a sediment map does not critically impact our results.

In summary, sensitivity tests show that uncertainties in topographic reconstructions impact the threshold CO₂ value for glacial onset. Accounting for ice streams—fast-moving ice regions, demonstrated to have been present from the Upper Ordovician record [Moreau et al., 2005; Denis et al., 2007; Le Heron et al., 2007; Ghienne et al., 2007; Denis et al., 2010; Ravier et al., 2015]—further demonstrates that the land ice volume simulated in our baseline runs has to be considered with caution. The two opposite and relatively extreme parameterizations tested here, leading to a weakly dynamic (in our baseline simulations) and especially dynamic (using the basal dragging dependence on sediments from Alvarez-Solas et al. [2011]) ice sheet, provide, respectively, maximum and minimum estimates for the volume of the Ordovician ice sheet. Irrespective of the topography and basal dragging parameterization, the modeled ice sheet extent, by contrast, turns out to be remarkably stable and thus reliable.

# 4. Discussion

## 4.1. Constraining Ice Sheet Extent

When accounting for the ice sheet feedbacks on Ordovician climate, the ice sheet expands to the low latitudes rather than being restricted in the polar area (≥ 60°S). This is an important result compared to the earlier models that omitted ocean-atmosphere-ice sheet interplays, such as those of Herrmann et al. [2003, 2004]. In our advanced models, a single ice sheet of large extent covers Gondwana from the South Pole to the middle-to-tropical paleolatitudes (Figure 3). The time series of ice sheet growth at 12 PAL, shown in Figure 8, reveals that highlands facilitate independent glacial centers to simultaneously nucleate over the Pole and at 60°S during the primary stages of glacial growth (Figures 8a–8c). However, such glacial centers are not initiated at middle-to-tropical paleolatitudes, where South Africa and South America were situated 445 Ma ago [Torsvik and Cocks, 2013], in any of our simulations. Furthermore, results from our Earth system model accounting for ice sheet feedback processes reveal that the local cooling associated with the growth of these glacial centers rapidly induces their coalescence into a single, continental-scale ice sheet (Figures 8c, 8d, and Figure 3). In our simulations, the persistence of independent ice sheets over Gondwana is not a sustainable land ice state. This result contradicts the hypothesized configuration that features multiple, individual ice sheets in North Africa, South Africa, and South America [e.g., Le Heron and Dowdeswell, 2009]. The sensitivity analysis we conducted (section 3.3) further reveals that the simulated ice sheet extent is a robust
<figure><img src="image_113.png" /><figcaption>**Figure 9.** Comparison of tropical sea surface temperatures (SSTs) reconstructed based on δ¹⁸O [Trotter et al., 2008] and Δ⁴⁷CO₂ [Finnegan et al., 2011] data with annual mean tropical (30°S–30°N) SSTs simulated for the Late Ordovician once ice sheets in steady state (red circles), using our asynchronous coupling method accounting for ice sheet feedback processes and for orbital variations. Modified after Trotter et al. [2008] and Finnegan et al. [2011]. *Dap.: Dapingian; Darriw.: Darriwilian; Sand.: Sandbian; H: Hirnantian; W: Wenlock.*</figcaption></figure>

feature that is critically dependent neither on the topographic reconstruction nor on a specific basal dragging parameterization. Our results are also supported by the study of *Lowry et al.* [2014], who simulate a single ice sheet over Gondwana with a different asynchronous coupling method and other climate and ice sheet models.

## 4.2. Glacial Onset in a Warm Climate

The buildup of the ice sheet, in response to a decrease in atmospheric forcing, is not linear. Between 16 PAL and 12 PAL, the ice sheet suddenly extends from the South Pole to the midlatitudes (Figure 3a), whereas tropical sea surface temperatures decrease by no more than 1.3°C (Figure 9). This nonlinear response suggests that the development of a large-scale ice sheet during the Ordovician could have occurred in conjunction with a moderate tropical sea surface temperatures signal, the amplitude of which would be within the error bars of routine geochemical analyses [e.g., Finnegan et al., 2011] and thus potentially below the detection limits. This further implies that warm tropical Ordovician seas [Trotter et al., 2008; Finnegan et al., 2011] and geological evidence for glacioeustatic events [Loi et al., 2010; Turner et al., 2012; Dabard et al., 2015; Rasmussen et al., 2016] are not mutually exclusive (30.3°C at 12 PAL, Figure 9). Interestingly, experiments using a lowered topography (section 3.3.1) are typified by a glacial onset occurring at 8 PAL (instead of 12 PAL). They show that the unexpected glacial inception in a warm climate does not result from the high relief characterization of our baseline runs.

## 4.3. Toward an “Early Paleozoic Ice Age”

*Trotter et al.* [2008] demonstrated that tropical SSTs steadily cooled through the Early Ordovician, reaching present-day levels that were subsequently sustained during the Middle and Late Ordovician, followed by a sudden cooling by ~7°C during the latest Ordovician (Figure 9). Several arguments support the validity of these data. *Nardin et al.* [2011] in particular, using a coupled climate-carbon model, showed that changes in Ordovician paleogeography, along with variations in the surface area of outcropping fresh volcanic rocks, induce a long-term decrease in atmospheric CO₂. This provides a potential mechanism for the cooling pattern demonstrated by *Trotter et al.* [2008]. *Trotter et al.'s* [2008] data are also in good agreement with other published geochemical data for the Late Ordovician and for the Silurian, respectively, from *Buggisch et al.* [2010] and *Wenzel et al.* [2000]. Moreover, recent clumped-isotope data from tropical latitudes
<figure><img src="image_114.png" /></figure>

**Figure 10.** Comparison of the Hirnantian glacial sedimentary record with the ice sheet extent simulated at 3 PAL using our coupling method accounting for ice sheet feedbacks on global climate and for orbital variations. Data from North Africa are from Ghienne et al. [2007], and data from South America and South Africa are from Díaz-Martínez and Grahn [2007]. Continent outlines are taken from Torsvik and Cocks [2009].

[Finnegan et al., 2011], although they are differently calibrated against the absolute temperature scale, display a similar temperature drop of 5 to 7°C during the latest Ordovician Hirnantian (Figure 9).

Figure 9 shows that the absolute values of tropical SSTs reconstructed by Finnegan et al. [2011] are systematically higher than those modeled in this study; this applies for our entire $pCO_2$ range, except for the Hirnantian where modeled and proxy-based SSTs align and correspond to a $pCO_2$ drop to 8 PAL in the model. Nonetheless, the generally unaligned SSTs (Figure 9) constitute a major discrepancy between Finnegan et al.'s [2011] data and our model runs. Finnegan et al. [2011] suggest the presence of large land ice volumes well before the Hirnantian, during the Katian, whereas the comparison between their SSTs and our simulations implies $CO_2$ levels that would have remained above 24 PAL and thus associated with ice-free conditions until the latest Ordovician, i.e., the Hirnantian. The modeled "Step 2" glaciation never occurred, based on these data. Simulated tropical SSTs are, however, in much better agreement with the range of absolute values presented by Trotter et al. [2008]. Comparison with Trotter's data suggests that ice sheet growth may have occurred as early as the Middle Ordovician (Figure 9). Taken at face value, our results, combined with these data, indicate a Darriwilian age for the glacial onset (between 16 and 12 PAL, step 1 in Figure 9). The second step in ice sheet expansion would then correspond to the Hirnantian glacial maximum

(3 PAL, step 2 in Figure 9). Our modeled results are also in agreement with the amplitude of drop in SSTs during the Hirnantian that is suggested in both of these proxy data studies [Trotter et al., 2008; Finnegan et al., 2011], when we decrease $CO_2$ from 8 PAL to 3 PAL.

Such a straight model-data comparison has to be interpreted cautiously. In addition to potential bias in the data (e.g., dating or preservation issues), simulated absolute values are highly model dependent and may also be biased due to large uncertainties in Ordovician boundary conditions (see section 3.3). Comparing our modeling results with the long-term Ordovician climatic trend proposed by Trotter et al. [2008] is nevertheless instructive, as it suggests a Darriwilian glacial onset, which agrees with recently published studies that report glacioeustatic signals in the Darriwilian rock record [Turner et al., 2012; Dabard et al., 2015; Rasmussen et al., 2016] and with paleontological studies that suggest pre-Sandbian cooling [Vandenbroucke et al., 2009], and it provides physical mechanisms for the latest Ordovician glacial maximum and associated sudden climate cooling [Trotter et al., 2008; Finnegan et al., 2011]. Besides this interesting match with geochemically inferred tropical SSTs, our Katian estimates (8–12 PAL, Figure 9) match relatively well the ~8 PAL $pCO_2$ proposed by Pancost et al. [2013] based on the analysis of isotope fractionation between organic and inorganic carbon during photosynthesis ($\epsilon_p$). Furthermore, the modeled spatial ice sheet extent at 3 PAL is overall consistent with the Hirnantian glacial sedimentary record (Figure 10). Discrepancies with shallow-shelf deposits from Turkey can be explained, as glacioeustatic emersion in such environments is not accounted for in our modeling procedure. In addition, using micropaleontological assemblage data to reconstruct climatic belts, Vandenbroucke et al. [2010a, 2010b] demonstrated that the Polar Front shifted in latitude from the Sandbian (460 Ma) to the Hirnantian (~440 Ma), reaching 40°S during the latest Ordovician. Comparing these Hirnantian data with the modeled climatic belts (Figure 11) reveals that again, the models best fit the data at 3 PAL. Taken
<figure><img src="image_115.png" /></figure>

**Figure 11.** Comparison of Ordovician proxy-derived climatic estimates with values simulated with our asynchronous coupling method once the ice sheet reached steady state with global climate. Southern Hemisphere climatic belts reconstructed based on Hirnantian chitinozoan marine zooplankton data [Vandenbroucke et al., 2010a, 2010b] are compared with climatic belts inferred from modeled SSTs. Abbreviations are plankton province names: *P: polar*; *sP: subpolar*; *tr: transitional*; *sT: subtropical*, and *T: tropical*. In our model predictions, the boundary between two plankton provinces is defined as the latitude at which the mean annual, latitudinally averaged SSTs simulated in the Southern Hemisphere reach the temperature threshold separating the two provinces in the present-day ocean. The reader is referred to Vandenbroucke et al. [2010a] for more details about the methods (see their Figure 1).

as a whole, our results suggest a protracted glaciation beginning as early as the Middle Ordovician Darriwilian and featuring a glacial maximum during the Hirnantian. By supporting the existence an ice sheet over the supercontinent Gondwana during most of the Middle and Late Ordovician, our results support the emerging paradigm of a protracted "Early Paleozoic Ice Age" [Page et al., 2007].

### 4.4. Ice Volume Analysis

We calculate the glacioeustatic sea level fall associated with growing ice sheets using the methods employed by Horton et al. [2007] and Horton and Poulsen [2009]. First, the equivalent volume of water is estimated by multiplying modeled volume of ice above the flotation level by the ice density value (0.917). Then, we convert the water equivalent to an isostatically adjusted sea level equivalent (IASLE), by dividing the water equivalent by the Late Ordovician ocean surface area (~413 Mkm² at 450 Ma) and subsequently correcting the result for the response of the oceanic lithosphere to seawater loading/unloading. To that end, we multiply the result by (1 - *k*), where *k* is the ratio of seawater density to oceanic lithosphere density (*k* = 0.284). Using results from our baseline runs, this simple calculation provides estimates of ~220 m between 12 and 8 PAL (step 1, Middle Ordovician) and ~250 m at 3 PAL (step 2, latest Ordovician Hirnantian). The sensitivity tests conducted in section 3.3.2 reveal that these values constitute uppermost estimates. The calculated minimum estimates are, respectively, 125 and 145 m.

Loi et al. [2010], using a backstripping procedure on the margin of Gondwana (Bou Ingarf section, Anti-Atlas, Southern Morocco), estimated an upper Hirnantian sea level fall of ~70 m. However, they explain that this sea level fall does not correctly represent the total ice volume stored over Gondwana during the latest Ordovician Hirnantian. It does not account for land ice stored over Gondwana before the Hirnantian climax. Loi et al. [2010] therefore estimate that the amplitude of the sea level fall would have reached nearly 150 m relative to the preglaciation (late Katian) datum [see also Haq and Schutter, 2008], and they note that this is a lower end estimate due to the incompleteness of the glacial sedimentary record. Moreover, recent studies [Raymo et al., 2011; Stocchi et al., 2013] demonstrated that the relationship between ice sheet growth and relative sea level fall is less linear than thought previously. Results obtained by Stocchi et al. [2013] on the onset of the Antarctic glaciation at the Eocene-Oligocene transition notably suggest that the near-field sea level change accompanying the growth of a continental-scale ice sheet is expected to strongly deviate from the global eustatic signal. Due to crustal deformation and gravitational perturbations, relative sea level around Antarctica rose in response to the buildup of the ice sheet. In the same way, the sea level fall documented by Loi et al. [2010] along the margin of Gondwana (150 m) may significantly underestimate Late Ordovician eustatic variations and should be considered as a lower end-member. Furthermore, recent clumped-isotope data from Finnegan et al. [2011] suggest that the Hirnantian ice volume may have reached twice that of the LGM, potentially
representing a water volume stored on land equivalent to 240 m sea level fall. It appears that our simulation run at 3 PAL eventually provides a total ice volume equivalent to 145 to 250 m sea level fall relative to ice-free conditions during the Hirnantian, thus in the same order of magnitude as most recent estimates inferred from sedimentological and geochemical data.

Ice volumes modeled between 12 and 8 PAL (125 to 220 m), on the contrary, are clearly much higher than the Middle Ordovician estimates reported by *Dabard et al.* [2015] using a backstripping procedure in the Armorican Massif in France: 50 to 80 m. The latter values seem relatively robust. They confirm previous estimates of 15–50 m (from Jordan, *Turner et al.* [2012]) and are in the same order of magnitude as the eustatic variations documented by *Haq and Schutter* [2008] during the same period of time, though somewhat smaller than the value ~150 m recently reported by *Rasmussen et al.* [2016] based on facies analysis of outcrop sections from Baltica. These discrepancies, i.e., our simulated land ice volume that is too high relative to the geological data, even when considering a slippery bed for the ice sheet (see section 3.3.2), are probably best explained by a significant overestimation of the modeled ice sheet extent. Studies on the last glaciation provide interesting insights into this issue. During the LGM, the Laurentide ice sheet suddenly extended to 40°N [*Ehlers and Gibbard, 2007*], similarly to the Gondwana ice sheet during the first step of Ordovician ice sheet growth. Simultaneously, an ice sheet developed in northern Europe. Northern Asia, on the other side, mysteriously remained ice free at the same latitudes. *Krinner et al.* [2006] suggested that high deposition rates of mineral dust in Asia, due to climate aridity under glacial conditions combined with the proximity to dust sources such as the Gobi desert, prohibited any perennial land ice cover over almost the entire Asian continent during the LGM. The deposition of dust may have similarly restricted the extension of the ice sheet over some regions of Gondwana. An additional cause for the overestimation may be found in our simplified modeling procedure, which imposes a steady $pCO_2$ forcing until ice sheet-climate equilibrium is reached. *Heinemann et al.* [2014] forced a coupled climate-ice sheet model with realistic orbital and greenhouse gas variations to simulate the Earth system evolution from 78 ka to 0 ka. They successfully simulated, in their transient runs, the Laurentide and Fennoscandian ice sheets evolution over this period of time and also the land ice extent that typifies the LGM. They subsequently ran a simulation by maintaining the LGM orbital configuration and $pCO_2$ for 90 kyr in order to reach climate-ice sheet equilibrium. Under these conditions, ice sheets extended down to the midlatitudes over both North America and the Eurasian continent, where it was previously restricted to Scandinavia. In the same way, it is reasonable to suppose that the Darriwilian ice sheet never reached full equilibrium due to orbital and/or $pCO_2$ variations. Such a smaller pre-Hirnantian ice sheet would be consistent with all the existing geological data. It is noteworthy that similar caveats were also suggested by *Horton et al.* [2007] in their simulation of the Late Paleozoic Ice Age. Last but not the least, strong uncertainties remain in paleogeographical and paleoenvironmental reconstructions in deep time [e.g., *Lees et al.*, 2002], and the base maps used in this study might, for instance, overestimate the continental surface situated between the South Pole and 45°S. This would lead to an overestimation of the area available for ice sheet growth and thus of the extent of the ice sheet developing over these latitudes.

### 4.5. Challenges

Modeling climate and land ice during the early Paleozoic involves numerous processes, which are sometimes difficult to constrain.

Recent studies demonstrated that orbital [*Ghienne et al.*, 2014] and greenhouse gas concentration [*Young et al.*, 2010] variations drove Late Ordovician ice sheet waxing and waning and also Middle Ordovician glacial inception [*Turner et al.*, 2012; *Dabard et al.*, 2015]. These phenomena have been shown to control Pleistocene glacial cycles as well [e.g., *Paillard*, 1998], but they cannot yet be precisely reconstructed for the Ordovician (e.g., *Laskar et al.* [2004] for orbital variations). Lack of constraints on these critical boundary conditions has been demonstrated (in the previous section) to significantly impact simulated ice sheet extent during the Middle Ordovician Darriwilian. This also prevents us to capture the high-frequency glacioeustatic variations recently reported in the Middle [*Turner et al.*, 2012; *Dabard et al.*, 2015] and Late [*Ghienne et al.*, 2014] Ordovician. In addition, the simulated ice sheet shows very little sensitivity to orbital and atmospheric forcing once it gets relatively large (Figure 1). Such a behavior was also demonstrated by *Herrmann et al.* [2003] during the Ordovician—although in a lesser extent due to the smaller land ice volume in play in their study—and by *Horton and Poulsen* [2009] in a study about the Late Paleozoic Ice Age. More generally, this strong hysteresis is a major issue in coupled ice sheet-climate studies [e.g., *Huybrechts*, 1993; *Pollard and DeConto*, 2005]. Numerous studies suggest that current ice sheet models may not be sensitive enough to correctly respond to climatic forcing. Melting Last Glacial Maximum ice sheets, for instance, sometimes requires
ad hoc processes, such as input of mineral dust [Peltier and Marshall, 1995; Bonelli et al., 2009] or tuning of the PDD coefficients [Charbit et al., 2013]. Ganopolski et al. [2010] and Robinson and Goelzer [2014] further demonstrated a significant contribution of insolation in surface melt processes, suggesting that the usual PDD melt scheme may constitute a major limitation to simulate the waning of paleo ice sheets. Moreover, recent studies suggest that critical mechanisms may be missing in the building of current ice sheet models, preventing them to simulate correctly the response of land ice to climate changes. The studies recently published by Pollard et al. [2015] and Gasson et al. [2016], in particular, show that supplementary ice sheet instability mechanisms, such as ice cliff failure and hydrofracturing, are required to correctly simulate the periodic melting of the East and West Antarctic ice sheets during the last ~25 million years. Austermann et al. [2015] even suggest that vertical deflections of the continental crust by mantle convective flow, i.e., dynamic topography, may have played a significant role in the Antarctic ice sheet fluctuations during the mid-Pliocene (~3 Ma). Horton et al. [2010], in a study about the Late Paleozoic Ice Age, demonstrate that the simulation of glacial-interglacial fluctuations is possible only when they include terrestrial vegetation feedbacks from high-latitude ecosystem changes. They show that changes in vegetation cover amplify the orbitally driven temperature variations. The latter mechanism, however, is not tenable in the Ordovician, when only lichens and mosses covered the land surface [Steemans et al., 2009; Rubinstein et al., 2010].

Beyond the difficulties to simulate the high-frequency temporal evolution of land ice, caveats in modeling the whole deep-time Earth system may induce important errors in continental ice cover numerical reconstruction. Our setup, in particular, does not account for the eustatic sea level fall associated with the growing Gondwana ice sheet, which requires up-to-date bathymetric reconstructions allowing to flood continental areas according to eustatic variations. Such reconstructions exist today [e.g., Vérard et al., 2015], but generally, they are not publicly available.

In summary, the up-to-date deep-time Earth system model employed in this study constitutes an important step forward in modeling Paleozoic climate and ice sheets. Our models notably provide relatively reliable first-order estimates of the temporal and spatial distribution of Ordovician land ice. Uncertainties in the boundary conditions, along with caveats in current ice sheet models, however, prohibit numerical investigation of the high-resolution patterns recently reported based on the Ordovician sedimentary record.

**5. Conclusions**

In this study, we employed a pioneering deep-time Earth system model to investigate the mechanisms driving Middle to Late Ordovician ice sheet growth, along with the equilibrium ice sheet extent and geometry. Using methods more advanced than earlier reconstructions, our model accounts for ice sheet feedback processes and for orbital variations at the Ordovician obliquity frequency (30 kyr). It was run across the wide spectrum of suggested Ordovician $pCO_2$ values, 24-3 PAL, resulting in the first simulation portraying the detailed evolution of Middle to Late Ordovician ice sheet growth. Forcing the well-established Ordovician cooling trend with decreasing $pCO_2$, the emergence of the ice sheet occurs in two steps. In a counterintuitive sequence of events, the vast continental ice sheet appears first and suddenly, covering Gondwana from the South Pole to the midpaleolatitudes. Only during the second act, temperatures drop steeply and sea ice extends. Our models therefore suggest that a single ice sheet extended from the South Pole (North Africa) to the tropics (South Africa) during the latest Ordovician. Scenarios with individual, separate ice centers are not stable and are rejected.

In spite of the particular Ordovician continental configuration potentially inducing a strong continentality over the South Pole, land ice growth is controlled by ablation rather than accumulation. Lack of constraints on Gondwana-wide surface lithology patterns and sediment cover leads to significant uncertainties in modeled land ice volumes. Sensitivity tests demonstrate that the simulated ice sheet spatial extent, however, is reliable and overly dependent neither on the topographic reconstruction nor on the basal dragging parameterization in use.

The sequence of events that we model indicates that glacial onset predated Late Ordovician climate cooling. The comparison with abundant sedimentological, geochemical, and micropaleontological data suggests that the first step in ice sheet growth possibly occurred as early as the Middle Ordovician Darriwilian, in agreement with recent studies reporting third-order eustatic cycles during the same period. Orbital and greenhouse gas variations would have prohibited the ice sheet from reaching long-term equilibrium at this time, explaining the moderate amplitude of reported glacioeustatic variations. The second step in ice sheet growth
is typified by an abrupt global climate cooling by ~14°C, leading to the further extension of the ice sheet into the tropical realm (30°S) during the Hirnantian glacial maximum. By demonstrating that the onset of a pre-Hirnantian continental-scale ice sheet over Gondwana requires no sudden climate cooling, our modeling results supply physical mechanisms for the previously enigmatic Late Ordovician stepped clumped-isotope data, which first display changes in ice volume and only later a drop in tropical sea surface temperatures. Our findings eventually support recent visions of a land ice cover over the supercontinent Gondwana throughout most of the Middle and Late Ordovician, embracing the emerging paradigm of an "Early Paleozoic Ice Age".

This study constitutes a major increment over previously published work that did not account for ice sheet feedback processes. Still, major shortcomings inherent to deep-time studies prevent us from capturing the high-frequency waxing and waning of the ice sheet. Better constrained Ordovician boundary conditions and extensive advances in modeling the cryosphere now constitute a prerequisite for further investigation.

## Acknowledgments

We thank David Pollard and an anonymous reviewer for careful and helpful reviews and Heiko Pälike for editorial handling. A.P. also thanks J.-F. Buoncristiani for fruitful discussions and valuable inputs. This work was granted access to the HPC resources of TGCC under the allocation 2014-012212 made by GENCI. T.R.A.V. acknowledges ANR support through grant ANR-12-B506-0014-SeqStrat-Ice. This research was funded by a CEA PhD grant CFR. This is a contribution to IGCP 591. The code of the FOAM model is available here http://www.mcs.anl.gov/research/projects/foam/index.html. LMDZ can be downloaded from the dedicated website as well http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele. The code of the GRISLI model is available on request to Catherine Ritz (LGGE, Grenoble). Requests for the model output can be sent to A.P. (alexandre.pohl@lsce.ipsl.fr) or Y.D. (yannick.donnadieu@lsce.ipsl.fr).

## References

* Alvarez-Solas, J., M. Montoya, C. Ritz, G. Ramstein, S. Charbit, C. Dumas, K. Nisancioglu, T. Dokken, and A. Ganopolski (2011), Heinrich event 1: An example of dynamical ice-sheet reaction to oceanic changes, *Clim. Past*, 7(4), 1297–1306.

* Alvarez-Solas, J., A. Robinson, M. Montoya, and C. Ritz (2013), Iceberg discharges of the last glacial period driven by oceanic circulation changes, *Proc. Natl. Acad. Sci. U.S.A.*, 110(41), 16350–16354.

* Amberg, C. E. A., T. Collart, W. Salenbien, L. M. Egger, A. Munnecke, A. T. Nielsen, C. Monnet, Ø. Hammer, and T. R. A. Vandenbroucke (2016), The nature of Ordovician limestone-marl alternations in the Oslo-Asker District (Norway): Witnesses of primary glacio-eustasy or diagenetic rhythms?, *Sci. Rep.*, 6, 18787.

* Austermann, J., D. Pollard, J. X. Mitrovica, R. Moucha, A. M. Forte, R. M. DeConto, D. B. Rowley, and M. E. Raymo (2015), The impact of dynamic topography change on Antarctic ice sheet stability during the mid-Pliocene warm period, *Geology*, 43(10), 927–930.

* Bennett, M. R. (2003), Ice streams as the arteries of an ice sheet: Their mechanics, stability and significance, *Earth-Sci. Rev.*, 61(3–4), 309–339.

* Berner, R. A. (1990), Atmospheric carbon dioxide levels over Phanerozoic time, *Science*, 249(4975), 1382–1386.

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

* Blakey, R. C. (2016), Colorado Plateau Geosystems. [Available at http://cpgeosystems.com.]

* Bonelli, S., S. Charbit, M. Kageyama, M. N. Woillez, G. Ramstein, C. Dumas, and A. Quiquet (2009), Investigating the evolution of major Northern Hemisphere ice sheets during the last glacial-interglacial cycle, *Clim. Past*, 5(3), 329–345.

* Brenchley, P. J., J. D. Marshall, G. Carden, D. Robertson, D. Long, T. Meidla, L. Hints, and T. F. Anderson (1994), Bathymetric and isotopic evidence for a short-lived Late Ordovician glaciation in a greenhouse period, *Geology*, 22(4), 295–298.

* Bueler, E., and J. Brown (2009), Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, *J. Geophys. Res.*, 114, F03008, doi:10.1029/2008JF001179.

* Buggisch, W., M. M. Joachimski, O. Lehnert, S. M. Bergström, J. E. Repetski, and G. F. Webers (2010), Did intense volcanism trigger the first Late Ordovician icehouse?, *Geology*, 38(4), 327–330.

* Charbit, S., C. Dumas, M. Kageyama, D. M. Roche, and C. Ritz (2013), Influence of ablation-related processes in the build-up of simulated Northern Hemisphere ice sheets during the last glacial cycle, *Cryosphere*, 7(2), 681–698.

* Clark, P. U., and D. Pollard (1998), Origin of the middle Pleistocene transition by ice sheet erosion of regolith, *Paleoceanography*, 13(1), 1–9.

* Crowley, T. J., and S. K. Baum (1991), Toward reconciliation of Late Ordovician (~440 Ma) glaciation with very high CO₂ levels, *J. Geophys. Res.*, 96(D12), 22,597–22,610.

* Crowley, T. J., and S. K. Baum (1995), Reconciling Late Ordovician (440 Ma) glaciation with very high (14X) CO₂ levels, *J. Geophys. Res.*, 100(D1), 1093–1101.

* Dabard, M. P., A. Loi, F. Paris, J.-F. Ghienne, M. Pistis, and M. Vidal (2015), Sea-level curve for the Middle to early Late Ordovician in the Armorican Massif (western France): Icehouse third-order glacio-eustatic cycles, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 436, 96–111.

* DeConto, R. M., and D. Pollard (2003), Rapid Cenozoic glaciation of Antarctica induced by declining atmospheric CO₂, *Nature*, 421(6920), 245–249.

* Denis, M., J.-F. Buoncristiani, M. Konaté, J.-F. Ghienne, and M. Guiraud (2007), Hirnantian glacial and deglacial record in SW Djado Basin (NE Niger), *Geodin. Acta*, 20(3), 177–195.

* Denis, M., M. Guiraud, M. Konaté, and J. F. Buoncristiani (2010), Subglacial deformation and water-pressure cycles as a key for understanding ice stream dynamics: Evidence from the Late Ordovician succession of the Djado Basin (Niger), *Int. J. Earth Sci.*, 99(6), 1399–1425.

* Diaz-Martinez, E., and Y. Grahn (2007), Early Silurian glaciation along the western margin of Gondwana (Peru, Bolivia and northern Argentina): Palaeogeographic and geodynamic setting, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 245(1), 62–81.

* Ehlers, J., and P. L. Gibbard (2007), The extent and chronology of Cenozoic global glaciation, *Quat. Int.*, 164–165, 6–20.

* Elrick, M., D. Reardon, W. Labor, J. Martin, A. Desrochers, and M. Pope (2013), Orbital-scale climate change and glacioeustasy during the early Late Ordovician (pre-Hirnantian) determined from δ¹⁸O values in marine apatite, *Geology*, 41(7), 775–778.

* Finnegan, S., K. Bergmann, J. M. Eiler, D. S. Jones, D. A. Fike, I. Eisenman, N. C. Hughes, A. K. Tripati, and W. W. Fischer (2011), The magnitude and duration of Late Ordovician-Early Silurian glaciation, *Science*, 331(6019), 903–906.

* François, L., A. Grard, and Y. Goddéris (2005), Modelling atmospheric CO₂ changes at geological time scales, in *Pre-Cambrian to Palaeozoic Palaeopolynology and Palaeobotany*, edited by P. Steemans and E. Javaux, Carnets de Géologie/Notebooks on Geology, Memoir 2005/02, Brest, France.

* Fretwell, P., et al. (2013), Bedmap2: Improved ice bed, surface and thickness datasets for Antarctica, *Cryosphere*, 7(1), 375–393.

* Ganopolski, A., R. Calov, and M. Claussen (2010), Simulation of the last glacial cycle with a coupled climate ice-sheet model of intermediate complexity, *Clim. Past*, 6, 229–244.

* Gasson, E., R. M. DeConto, D. Pollard, and R. H. Levy (2016), Dynamic Antarctic ice sheet during the early to mid-Miocene, *Proc. Natl. Acad. Sci. U.S.A.*, 113(13), 3459–3464, doi:10.1073/pnas.1516130113.

* Ghienne, J.-F., D. P. Le Heron, J. Moreau, M. Denis, and M. Deynoux (2007), The Late Ordovician glacial sedimentary system of the North Gondwana platform, in *Glacial Sedimentary Processes and Products, Spec. Publ. Int. Assoc. of Sedimentologists*, edited by M. Hambrey et al., pp. 295–319, Blackwells, Oxford, U. K.
Ghienne, J.-F., et al. (2014), A Cenozoic-style scenario for the end-Ordovician glaciation, Nat. Commun., 5, 4485.

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

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

Haq, B. U., and S. R. Schutter (2008), A chronology of paleozoic sea-level changes, Science, 322(5898), 64–68.

Heimbach, P., and M. Losch (2012), Adjoint sensitivities of sub-ice-shelf melt rates to ocean circulation under the Pine Island Ice Shelf, West Antarctica, Ann. Glaciol., 53(60), 59–69.

Heinemann, M., A. Timmermann, O. Elison Timm, F. Saito, and A. Abe-Ouchi (2014), Deglacial ice sheet meltdown: Orbital pacemaking and CO₂ effects, Clim. Past, 10(4), 1567–579.

Herrmann, A. D., M. E. Patzkowsky, and D. Pollard (2003), Obliquity forcing with 8–12 times preindustrial levels of atmospheric pCO₂ during the Late Ordovician glaciation, Geology, 31(6), 485–488.

Herrmann, A. D., M. E. Patzkowsky, and D. Pollard (2004), The impact of paleogeography, pCO₂, poleward ocean heat transport and sea level change on global cooling during the Late Ordovician, Palaeogeogr. Palaeoclimatol. Palaeoecol., 206(1), 59–74.

Horton, D. E., and C. J. Poulsen (2009), Paradox of late Paleozoic glacioeustasy, Geology, 37(8), 715–718.

Horton, D. E., C. J. Poulsen, and D. Pollard (2007), Orbital and CO₂ forcing of late Paleozoic continental ice sheets, Geophys. Res. Lett., 34, L19708, doi:10.1029/2007GL031188.

Horton, D. E., C. J. Poulsen, and D. Pollard (2010), Influence of high-latitude vegetation feedbacks on late Palaeozoic glacial cycles, Nat. Geosci., 3(8), 572–577.

Hourdin, F., et al. (2013), Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model, Clim. Dyn., 40, 2167–2192.

Huybrechts, P. (1993), Glaciological modelling of the Late Cenozoic East Antarctic ice sheet: Stability or dynamism?, Geogr. Ann., 75(A), 221–238.

Jacob, R. L. (1997), Low frequency variability in a simulated atmosphere ocean system, PhD thesis, Univ. of Wisconsin-Madison, Madison, Wisconsin.

Krinner, G., O. Boucher, and Y. Balkanski (2006), Ice-free glacial northern Asia due to dust deposition on snow, Clim. Dyn., 27(6), 613–625.

Ladant, J. B., Y. Donnadieu, V. Lefebvre, and C. Dumas (2014), The respective role of atmospheric carbon dioxide and orbital parameters on ice sheet evolution at the Eocene-Oligocene transition, Paleceanography, 29, 810–823, doi:10.1002/2013PA002593.

Laskar, J., P. Robutel, F. Joutel, M. Gastineau, A. Correia, and B. Levrard (2004), A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428(1), 261–285.

Le Heron, D. P., and J. A. Dowdeswell (2009), Calculating ice volumes and ice flux to constrain the dimensions of a 440 Ma North African ice sheet, J. Geol. Soc., 166(2), 277–281.

Le Heron, D. P., J.-F. Ghienne, M. El Houicha, Y. Khoukhi, and J.-L. Rubino (2007), Maximum extent of ice sheets in Morocco during the Late Ordovician glaciation, Palaeogeogr. Palaeoclimatol. Palaeoecol., 245(1–2), 200–226.

Lees, D. C., R. A. Fortey, and L. R. M. Cocks (2002), Quantifying paleogeography using biogeography: A test case for the Ordovician and Silurian of Avalonia based on brachiopods and trilobites, Paleobiology, 28(3), 343–363.

Licht, A., et al. (2014), Asian monsoons in a late Eocene greenhouse world, Nature, 513(7519), 501–506.

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

Loi, A., et al. (2010), The Late Ordovician glacio-eustatic record from a high-latitude storm-dominated shelf succession: The Bou Ingarf section (Anti-Atlas, Southern Morocco), Palaeogeogr. Palaeoclimatol. Palaeoecol., 296(3–4), 332–358.

Moreau, J., J.-F. Ghienne, D. P. L. Heron, J.-L. Rubino, and M. Deynoux (2005), 440 Ma ice stream in North Africa, Geology, 33(9), 753–756.

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

Nielsen, A. T. (2004), Ordovician sea-level changes: A Baltoscandian perspective, in The Great Ordovician Biodiversification Event, edited by B. D. Webby et al., pp. 84–93, Columbia Univ. Press, New York.

Page, A., J. A. Zalasiewicz, M. Williams, and L. E. Popov (2007), Were transgressive black shales a negative feedback modulating glacioeustasy in the Early Palaeozoic Icehouse?, in Deep-Time Perspectives on Climate Change: Marrying the Signal from Computer Models and Biological Proxies, edited by M. Williams et al., pp. 123–156, The Micropalaeontological Soc., Spec. Publ., The Geolog. Soc., London, U. K.

Paillard, D. (1998), The timing of Pleistocene glaciations from a simple multiple-state climate model, Nature, 391, 378–381.

Pancost, R. D., K. H. Freeman, A. D. Herrmann, M. E. Patzkowsky, L. Ainsaar, and T. Martma (2013), Reconstructing Late Ordovician carbon cycle variations, Geochim. Cosmochim. Acta, 105, 433–454.

Peltier, W. R., and S. Marshall (1995), Coupled energy-balance/ice-sheet model simulations of the glacial cycle: A possible connection between terminations and terrigenous dust, J. Geophys. Res., 100(D7), 14,269–14,289.

Peyaud, V., C. Ritz, and G. Krinner (2007), Modelling the Early Weichselian Eurasian Ice Sheets: Role of ice shelves and influence of ice-dammed lakes, Clim. Past, 3, 375–386.

Pohl, A., Y. Donnadieu, G. Le Hir, J. F. Buoncristiani, and E. Vennin (2014), Effect of the Ordovician paleogeography on the (in)stability of the climate, Clim. Past, 10(6), 2053–2066.

Pohl, A., E. Nardin, T. Vandenbroucke, and Y. Donnadieu (2015), High dependence of Ordovician ocean surface circulation on atmospheric CO₂ levels, Palaeogeogr. Palaeoclimatol. Palaeoecol., doi:10.1016/j.palaeo.2015.09.036, in press.

Pollard, D. (2010), A retrospective look at coupled ice sheet-climate modeling, Clim. Change, 100(1), 173–194.

Pollard, D., and R. M. DeConto (2005), Hysteresis in Cenozoic Antarctic ice-sheet variations, Global Planet. Change, 45(1–3), 9–21.

Pollard, D., and R. M. DeConto (2012), Description of a hybrid ice sheet-shelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1275.

Pollard, D., R. M. DeConto, and R. B. Alley (2015), Potential Antarctic Ice Sheet retreat driven by hydrofracturing and ice cliff failure, Earth Planet. Sci. Lett., 412, 112–121.

Quiquet, A., C. Ritz, H. J. Punge, and D. Salas y Mélia (2013), Greenland ice sheet contribution to sea level rise during the last interglacial period: A modelling study driven and constrained by ice core data, Clim. Past, 9(1), 353–366.

Rasmussen, C. M. Ø., et al. (2016), Onset of main Phanerozoic marine radiation sparked by emerging Mid Ordovician icehouse, Sci. Rep., 6, 18884.

Ravier, E., J. F. Buoncristiani, J. Menzies, M. Guiraud, S. Clerc, and E. Portier (2015), Does porewater or meltwater control tunnel valley genesis? Case studies from the Hirnantian of Morocco, Palaeogeogr. Palaeoclimatol. Palaeoecol., 418, 359–376.

Raymo, M. E., J. X. Mitrovica, M. J. O'Leary, R. M. DeConto, and P. J. Hearty (2011), Departures from eustasy in Pliocene sea-level records, Nat. Geosci., 4(5), 328–332.
Reeh, N. (1991), Parameterization of melt rate and surface temperature on the Greenland ice sheet, *Polarforschung*, 59(1), 113–128.

Ritz, C., V. Rommelaere, and C. Dumas (2001), Modeling the evolution of Antarctic ice sheet over the last 420,000 years—Implications for altitude changes in the Vostok region, *J. Geophys. Res.*, 106(D23), 31,943–31,964.

Robinson, A., and H. Goelzer (2014), The importance of insolation changes for paleo ice sheet modeling, *Cryosphere*, 8(4), 1419–1428.

Rothman, D. H. (2002), Atmospheric carbon dioxide levels for the last 500 million years, *Proc. Natl. Acad. Sci. U.S.A.*, 99(7), 4167–4171.

Rubinstein, C. V., P. Gerrienne, G. S. de la Puente, R. A. Astini, and P. Steemans (2010), Early Middle Ordovician evidence for land plants in Argentina (eastern Gondwana), *New Phytol.*, **188**(2), 365–369.

Scotese, C. R. (2016), PALEOMAP Project. [Available at http://www.scotese.com.]

Scotese, C. R., and W. S. McKerrow (1991), Ordovician plate tectonic reconstructions, in *Advances in Ordovician Geology*, edited by C. R. Scotese and W. S. McKerrow, pp. 271–282, Geol. Surv. of Canada, Ottawa, Ontario.

Semtner, A. J. (1976), A model for the thermodynamic growth of sea ice in numerical investigations of climate, *J. Phys. Oceanogr.*, 6(3), 379–389.

Sheehan, P. M. (2001), The late Ordovician mass extinction, *Annu. Rev. Earth Planet. Sci.*, 29(1), 331–364.

Steemans, P., A. Le Herissé, J. Melvin, M. A. Miller, F. Paris, J. Verniers, and C. H. Wellman (2009), Origin and radiation of the earliest vascular land plants, *Science*, 324(5925), 353–353.

Stocchi, P., C. Escutia, A. J. P. Houben, B. L. A. Vermeersen, P. K. Bijl, H. Brinkhuis, R. M. DeConto, S. Galeotti, S. Passchier, and D. Pollard (2013), Relative sea-level rise around East Antarctica during Oligocene glaciation, *Nat. Geosci.*, 6(5), 380–384.

Tabor, C. R., and C. J. Poulsen (2016), Simulating the mid-Pleistocene transition through regolith removal, *Earth Planet. Sci. Lett.*, 434(C), 231–240.

Tobin, K. J., and S. M. Bergström (2002), Implications of Ordovician (460 Myr) marine cement for constraining seawater temperature and atmospheric $pCO_2$, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 181(4), 399–417.

Tobin, K. J., S. M. Bergström, and P. De La Garza (2005), A mid-Caradocian (453 Ma) drawdown in atmospheric $pCO_2$ without ice sheet development?, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 226(3–4), 187–204.

Torsvik, T. H., and L. R. M. Cocks (2009), BugPlates: Linking biogeography and palaeogeography, software manual. [Available at http://www.geodynamics.no.]

Torsvik, T. H., and L. R. M. Cocks (2013), New global palaeogeographical reconstructions for the Early Palaeozoic and their generation, in *Early Palaeozoic Biogeography and Palaeogeography, Geol. Soc. of London Memoirs*, edited by D. A. T. Harper and T. Servais, pp. 5–24, The Geol. Soc., London, U. K.

Trotter, J. A., I. S. Williams, C. R. Barnes, C. Lécuyer, and R. S. Nicoll (2008), Did cooling oceans trigger Ordovician biodiversification? Evidence from conodont thermometry, *Science*, 321, 550–554.

Turner, B. R., H. A. Armstrong, C. R. Wilson, and I. M. Makhlouf (2012), High frequency eustatic sea-level changes during the Middle to early Late Ordovician of southern Jordan: Indirect evidence for a Darriwilian Ice Age in Gondwana, *Sediment. Geol.*, 251, 34–48.

Vandenbroucke, T. R. A., H. A. Armstrong, M. Williams, J. A. Zalasiewicz, and K. Sabbe (2009), Ground-truthing Late Ordovician climate models using the paleobiogeography of graptolites, *Paleoceanography*, 24, PA4202, doi:10.1029/2008PA001720.

Vandenbroucke, T. R. A., H. A. Armstrong, M. Williams, F. Paris, J. A. Zalasiewicz, K. Sabbe, J. Nølvak, T. J. Challands, J. Verniers, and T. Servais (2010a), Polar front shift and atmospheric $CO_2$ during the glacial maximum of the Early Paleozoic Icehouse, *Proc. Natl. Acad. Sci. U.S.A.*, 107(34), 14,983–14,986.

Vandenbroucke, T. R. A., H. A. Armstrong, M. Williams, F. Paris, K. Sabbe, J. A. Zalasiewicz, J. Nølvak, and J. Verniers (2010b), Epipelagic chitinozoan biotopes map a steep latitudinal temperature gradient for earliest Late Ordovician seas: Implications for a cooling Late Ordovician climate, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 294(3–4), 202–219.

Vérard, C., C. Hochard, P. O. Baumgartner, G. M. Stampfli, and M. Liu (2015), 3D palaeogeographic reconstructions of the Phanerozoic versus sea, *J. Palaeogeogr.*, 4(1), 64–84.

Wenzel, B., C. Lécuyer, and M. M. Joachimski (2000), Comparing oxygen isotope records of Silurian calcite and phosphate—$\delta^{18}O$ compositions of brachiopods and conodonts, *Geochim. Cosmochim. Acta*, 64(11), 1859–1872.

Williams, G. E. (1991), Milankovitch-band cyclicity in bedded halite deposits contemporaneous with Late Ordovician-Early Silurian glaciation, Canning Basin, Western Australia, *Earth Planet. Sci. Lett.*, **103**(1), 143–155.

Yapp, C. J., and H. Poths (1992), Ancient atmospheric $CO_2$ pressures inferred from natural goethites, *Nature*, 355, 342–344.

Yapp, C. J., and H. Poths (1996), Carbon isotopes in continental weathering environments and variations in ancient atmospheric $CO_2$ pressure, *Earth Planet. Sci. Lett.*, **137**(1), 71–82.

Young, S. A., M. R. Saltzman, K. A. Poland, J. S. Linder, and L. R. Kump (2009), A major drop in seawater $^{87}Sr/^{86}Sr$ during the Middle Ordovician (Darriwilian): Links to volcanism and climate?, *Geology*, 37(10), 951–954.

Young, S. A., M. R. Saltzman, W. I. Ausich, A. Desrochers, and D. Kaljo (2010), Did changes in atmospheric $CO_2$ coincide with latest Ordovician glacial-interglacial cycles?, *Palaeogeogr. Palaeoclimatol. Palaeoecol.*, 296(3–4), 376–388.
# A warm or a cold early Earth? New insights from a 3-D climate-carbon model

Benjamin Charnay<sup>a,b</sup>, Guillaume Le Hir<sup>c</sup>, Frédéric Fluteau<sup>c</sup>, François Forget<sup>d</sup>, David C. Catling<sup>b,e</sup>

<sup>a</sup>LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 Place Jules Janssen, 92195 Meudon, France

<sup>b</sup>Virtual Planetary Laboratory, University of Washington, Seattle, WA 98125, USA.

<sup>c</sup>Institut de Physique du Globe de Paris, Sorbonne Paris Cit, Universit Paris Diderot, 75238 Paris cedex 05, France

<sup>d</sup>Laboratoire de Météorologie Dynamique, IPSL/CNRS/UPMC, Paris 75005, France.

<sup>e</sup>Department of Earth and Space Science, University of Washington, Seattle, WA 98125, USA.

## Abstract

Oxygen isotopes in marine cherts have been used to infer hot oceans during the Archean with temperatures between 60°C (333 K) and 80°C (353 K). Such climates are challenging for the early Earth warmed by the faint young Sun. The interpretation of the data has therefore been controversial. 1D climate modeling inferred that such hot climates would require very high levels of CO₂ (2-6 bars). Previous carbon cycle modeling concluded that such stable hot climates were impossible and that the carbon cycle should lead to cold climates during the Hadean and the Archean. Here, we revisit the climate and carbon cycle of the early Earth at 3.8 Ga using a 3D climate-carbon model. We find that CO₂ partial pressures of around 1 bar could have produced hot climates given a low land fraction and cloud feedback effects. However, such high CO₂ partial pressures should not have been stable because of the weathering of terrestrial and oceanic basalts, producing an efficient stabilizing feedback. Moreover, the weathering of impact ejecta during the Late Heavy Bombardment (LHB) would have strongly reduced the CO₂ partial pressure leading to cold climates and potentially snowball Earth events after large impacts. Our results therefore favor cold or temperate climates with global mean temperatures between around 8°C (281 K) and 30°C (303 K) and with 0.1-0.36 bar of CO₂ for the late Hadean and early Archean. Finally, our model suggests that the carbon cycle was efficient for preserving clement conditions on the early Earth without necessarily requiring any other greenhouse gas or warming process.

**Keywords:** early Earth, climate, carbon cycle, Hadean, Archean, Late Heavy Bombardment

## 1. Introduction

Life likely emerged on Earth before 3.5 Ga, during the Hadean or the early Archean [Buick et al., 1981; Nisbet and Sleep, 2001; Bell et al., 2015]. Studying the atmosphere and the climates of the Earth during its first billion year is critical for understanding the environment in which life emerged and developed. Oxygen isotopic ratio of Archean marine cherts have been interpreted to suggest hot oceans with temperatures between 60°C (333 K) and 80°C (353 K) [Knauth and Lowe, 2003; Robert and Chaussidon, 2006], compatible with the thermophiles inferred from evolutionary models of the ancient life [Gaucher et al., 2008] and possibly with the indications for a low ocean water viscosity [Fralick and Carter, 2011]. However, such an interpretation has been strongly debated [Kasting and Howard, 2006; van den Boorn et al., 2007; Marin-Carbonne et al., 2014; Tartese et al., 2016] and other analyses suggest temperate oceans with temperatures lower than 40°C [Hren et al., 2009; Blake et al., 2010]. Moreover, at 3.5 Ga, the presence of glacial rocks at 20-40° latitude implies temperatures below 20°C [de Wit and Furnes, 2016], at least episodically.

In addition, previous modeling of the carbon cycle on the early Earth by Sleep and Zahnle [2001] and Zahnle and Sleep [2002] suggested that the Archean was cold unless another strong greenhouse gas was present. They also suggested that the Hadean was likely very cold because of the weathering of impact ejecta, particularly during the Late Heavy Bombardment (LHB). Therefore, the temperature of the early oceans remains an open question.

Here, we use a 3D climate-carbon model with pCO₂ ranging from 0.01 to 1 bar in order to determine physically plausible climates of the late Hadean and early Archean Earth. In section 2, we describe the climate simulations. In
section 3 and 4, we analyze the carbon cycle and the effect of the LHB with respect to carbon cycle responses. We finish with a conclusion in section 5.

## 2. Climate modeling

### 2.1. Description of the model

We simulated the atmosphere of the early Earth using the Generic LMDZ GCM (Global Climate Model). This model solves the primitive hydrostatic equations of meteorology using a finite difference dynamical core on an Arakawa C grid. It uses robust and general parameterizations in order to simulate planets very different from the present-day Earth. The Generic LMDZ GCM has already been used for studying the temperate and cold climates of the Archean Earth [Charnay et al., 2013], the climates of early Mars with high CO₂ pressures [Forget et al., 2013], and the runaway greenhouse effect on terrestrial planets with high amounts of water vapor [Leconte et al., 2013]. It is thus adapted for simulating the early Earth with high amounts of CO₂ and potentially hot and moist climates.

The radiative scheme is based on the correlated-k method. At a given pressure and temperature, correlated-k coefficients in the GCM are interpolated from a matrix of coefficients stored in a 7 × 9 temperature and log-pressure grid: T=100, 150, 200, 250, 300, 350, 400 K, p = 10⁻¹, 10⁰, 10¹, ..., 10⁷ Pa. We used 36 spectral bands in the thermal infrared and 38 at solar wavelengths. Sixteen points were used for the g-space integration, where g is the cumulative distribution function of the absorption data for each band.

Simulations were performed with a horizontal resolution of 64×48 (corresponding to resolutions of 3.75° latitude by 5.625° longitude) and with 25 vertical layers with the lowest midlayer level at 5 m and the top level at 0.5 hPa. For the cloud microphysics, we either fixed the radii of water cloud particles (e.g. 12 µm for liquid droplets and 35 µm for icy particles for present-day Earth) or we fixed the density of cloud condensation nuclei (CCN) by mass of air (e.g. 5×10⁶ particles/kg for liquid droplets and 2×10⁴ particles/kg for icy clouds for present-day Earth). As in Charnay et al. [2013], the oceanic transport and the sea ice formation were computed with the 2-layer oceanic model from Codron [2012].

To investigate the early Earth at 3.8 Ga, we used a solar constant of 1024 W/m² corresponding to 75% of the present value (1361 W/m²). We ran simulations with no land. This hypothesis is valid for small continental surface fractions (e.g. <30% of the present-day fraction or <10% of the total Earth's surface) as we assumed at this time [Flament et al., 2008] (see also section 3.1), although the land fraction is a matter of dispute [Viehmann et al., 2014]. Such small continental fraction should indeed have a negligible direct impact on the global climate. We also assumed that Earth's rotation period was 14h. We used an atmospheric composition with 1 bar of nitrogen, a partial pressure of CO₂ (pCO₂) from 0.01 to 1 bar, and either no CH₄ or 2 mbar of CH₄, spanning the range of plausible methane concentrations estimated by a model for early ecosystems [Kharecha et al., 2005]. Wordsworth et al. [2017] showed that CO₂-CH₄ collision-induced absorption could have produced a strong warming on early Mars if the CH₄ mixing ratio was higher than 1%. With our CH₄ mixing ratio of 0.1-0.2 %, CO₂-CH₄ CIA produce a warming lower than 0.1 K with the 1D version of the model. We thus neglected CO₂-CH₄ CIA in our 3D simulations.

### 2.2. Results

Figure 1 shows the global mean and zonal mean surface temperatures obtained with the GCM for present-day properties of clouds (Global mean values are gathered in table 1). Using 3.8 Ga conditions, the Earth falls into a full glaciation for the present-day pCO₂ (see Charnay et al. [2013]), but it is never fully ice-covered with at least 0.01 bar of CO₂. For this case with no methane, the mean surface temperature is -11.8°C, well below the freezing point, but there is still a cold ice-free equatorial belt between -37°N and +37°N. With 0.1 bar of CO₂, the climate is temperate (global mean surface temperature around 7.4°C without CH₄ and 21.3°C with CH₄). Ice-free conditions are obtained for pCO₂ higher than around 0.3 bar without methane or slightly higher than 0.1 bar with methane, corresponding to a mean surface temperature higher than around 23°C. Our model has a climate sensitivity with pCO₂ very similar to the GCM from Wolf and Toon [2013] going up to 0.3 bar. With 0.5 bar of CO₂, the climate becomes warm (mean surface temperature around 34.7°C without CH₄ and 49.9°C with CH₄). With 1 bar of CO₂, the mean surface temperature is around 59.5°C without CH₄ and around 67.5°C with CH₄.

Ozak et al. [2016] showed that CO₂ collisional line mixing reduces the atmospheric opacity between 400 and 550 cm⁻¹ producing a significant cooling (~10°C) for a 1 bar CO₂ atmosphere for early Mars. Since we do not include
<table><caption>Table 1: Global mean values from the GCM. Columns from left to rights are: pCO₂ (in bar), pCH₄ (in mbar), the mean surface temperature (in °C), the mean planetary albedo, the mean precipitation rate (in mm/day), the shortwave cloud forcing (in W/m²) and the longwave cloud forcing (in W/m²).</caption><thead><tr><th>pCO<sub>2</sub><br>(bar)</th><th>pCH<sub>4</sub><br>(mbar)</th><th>Tsurf (°C)</th><th>Planetary<br>albedo</th><th>Precipitation<br>(mm/day)</th><th>SW forcing<br>(W/m<sup>2</sup>)</th><th>LW forcing<br>(W/m<sup>2</sup>)</th></tr></thead><tbody><tr><td>0.01</td><td>0</td><td>-11.8</td><td>0.4</td><td>1.8</td><td>-38.3</td><td>+24.5</td></tr><tr><td>0.01</td><td>2</td><td>2.7</td><td>0.33</td><td>2.2</td><td>-42.0</td><td>+20.2</td></tr><tr><td>0.1</td><td>0</td><td>7.4</td><td>0.34</td><td>2.8</td><td>-47.1</td><td>+24.9</td></tr><tr><td>0.1</td><td>2</td><td>21.3</td><td>0.28</td><td>3.5</td><td>-39.8</td><td>+18.1</td></tr><tr><td>0.3</td><td>0</td><td>22.1</td><td>0.31</td><td>3.6</td><td>-42.4</td><td>+20.8</td></tr><tr><td>0.3</td><td>2</td><td>34.7</td><td>0.26</td><td>4.4</td><td>-32.6</td><td>+13.8</td></tr><tr><td>0.5</td><td>0</td><td>34.7</td><td>0.27</td><td>4.6</td><td>-30.4</td><td>+16.0</td></tr><tr><td>0.5</td><td>2</td><td>49.9</td><td>0.23</td><td>5.1</td><td>-22.8</td><td>+12.3</td></tr><tr><td>1.0</td><td>0</td><td>59.5</td><td>0.24</td><td>5.7</td><td>-19.4</td><td>+14.2</td></tr><tr><td>1.0</td><td>2</td><td>67.5</td><td>0.23</td><td>6.1</td><td>-20.3</td><td>+14.0</td></tr></tbody></table>

collisional line mixing, our model might overestimate temperatures for high pCO₂. However, for the early Earth with a 1 bar CO₂ atmosphere, the peak of thermal emission from the surface would be around 1110 cm⁻¹ (compared to around 730 cm⁻¹ for early Mars), so quite far from the 400-550 cm⁻¹ spectral region. We therefore expect a weaker cooling effect by CO₂ collisional line mixing for the early Earth than for early Mars.

Fundamental changes to cloud properties on early Earth have been proposed as alternative solutions to the faint young Sun paradox. In particular, it has been proposed that cloud particles were larger on the early Earth because of less cloud condensation nuclei from biological sources [Rosing et al., 2010]. In figure 2, we tested the impact of larger liquid cloud particles (17 µm versus 12 µm). Larger particles warm the planet by decreasing the planetary albedo [Charnay et al., 2013]. This effect is particularly strong for temperate climate with a warming of around 8°C for pCO₂=0.1 bar. However, it becomes weak for warmer climates. We also simulated the climates with a fixed amount of CCN and varying cloud radii (black line in left panel in figure 2). The mean surface temperature is almost identical to the case with 12 µm particles.

The right panel in figure 2 shows the shortwave (albedo cooling effect) and longwave (greenhouse warming effect) radiative forcing by clouds depending on pCO₂ for the different particle size parametrizations. The trends are similar for all particles size parametrizations. The net cloud radiative forcing is higher for cold climates (e.g. pCO₂=0.01 bar) than for temperate climates (e.g. pCO₂=0.1 bar). This is due to the lower albedo cooling of clouds over sea ice and to a reduction of the cloud cover over the cold equatorial ocean, particularly close to the freezing line (see Charnay et al. [2013]). The net cloud radiative forcing increases with pCO₂ when pCO₂ exceeds 0.1 bar. The albedo cooling effect decreases because of a reduction of lower clouds (see right panel in figure 3). At the same time, the greenhouse warming effect decreases. This is mostly due to the gaseous infrared emission level (the altitude where most of the thermal emission is emitted) going up with more CO₂, which reduces the greenhouse effect of clouds. The cloud albedo cooling decreases faster than the cloud greenhouse warming, which leads to a net forcing increasing with pCO₂. Clouds therefore produce a destabilizing feedback on the climate. For hot climates, they tend to have a small net radiative effect (their cooling effect being reduced by around 30 W/m²), which depends little on the radii of cloud particles.

In response to clouds properties, our GCM predicts that hot climates can be obtained with 1 bar or less of CO₂ at 3.8 Ga. This required amount of CO₂ is significantly lower than previous 1D estimations by Kasting and Howard [2006]. They found that if there were no methane, around 3 bars of CO₂ would be required for getting a mean surface temperature of 60°C at 3.3 Ga. This discrepancy is due to cloud feedbacks, not simulated in the 1D model, and to a few other differences in the climate modeling. In particular, they used the present-day land cover while we assumed no continent, and we used more recent opacity data, likely increasing the visible and infrared gaseous absorption. From a purely climate modeling point of view (i.e., ignoring long-term geochemical feedbacks), a warm early Earth with temperature of 60-70°C is possible with sufficient CO₂ if there is no land or almost no land.
<figure><img src="image_116.png" /><figcaption>Figure 1: Top panel: global mean surface temperature without methane (blue line) and with 2 mbars of methane (red line). Bottom panel: zonal mean surface temperature for an atmosphere with no methane and with 0.01, 0.1, 0.3, 0.5 and 1 bar of CO₂ for lines from light gray to black respectively. The blue dashed line is the freezing temperature of water.</figcaption></figure>

# 3. Carbon cycle modeling

## 3.1. *Model set-up and boundary conditions*

In order to quantify the carbon cycle equilibrium response to simulated climate, we used the model from Le Hir *et al.* [2008a, b]. We generally assumed a classical plate tectonic regime as in *Sleep and Zahnle* [2001] (see Fig. 1 in Appendix A), although there is debate about when a modern style of plate tectonics began [Stern, 2007]. CO₂ reservoirs are 1) the atmosphere, 2) the surface/deep ocean, 3) the mantle, 4) carbonates on continental platforms and 5) calcite precipitation in the oceanic crust resulting from seafloor weathering. The inventory of carbon between
<figure><img src="image_117.png" /></figure>

Figure 2: Left panel: global mean surface temperature with no methane with 12 µm liquid cloud droplets (blue), 17 µm liquid cloud droplets (red) and varying cloud droplet radii with a fixed CCN concentration (black line). Right: shortwave cloud forcing (dotted lines), longwave cloud forcing (dashed lines) and net cloud forcing (solid lines). The colors correspond to the same cases as the left panel.

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

Figure 3: Left panel: vertical profile of global mean temperature. Middle panel: vertical profile of global mean water vapor mixing ratio (in kg/kg). Right panel: vertical profile of global mean condensed water mixing ratio (in $10^{-5}$ kg/kg). For the three panels, pCO₂ is 0.01 (dotted line), 0.1 (dashed-dotted line), 0.5 (dashed line) and 1 (solid line) bar and there is no methane.

these six reservoirs evolves through CO₂ outgassing by mid oceanic ridges and arc volcanoes, continental and seafloor weathering, carbonate formation/dissolution, subduction and recycling. The part of subducted carbon released by arc volcanoes depends on the recycling factor. Based on a steady state assumption, we considered an instantaneous carbon recycling in subduction zones. The model parametrizations for these processes are described in Appendix A.

We simulated the carbon cycle using the surface temperatures and precipitation-evaporation rates (P-E) from the GCM. In the absence of constraints about geography, we assumed that emerged continents are homogeneously spread over all the Earth as small continents (each grid cell has its own fraction of continents). That tends to maximize the continental weathering. We interpolated GCM temperatures and Precipitation-Evaporation (P-E) between pCO₂ values at each grid point. We parametrized the seafloor weathering with a dependence on the deep ocean temperature (see Coogan and Gillis [2013] for geological evidence and Krissansen-Totton and Catling [2017] for application for the last 100 Ma), assumed to be equal to sea surface temperature at 60° latitude. Most of the weathering occuring in the upper part of the oceanic crust (<1km), the seafloor weathering is explicitly described between 0 and 500 m by 6
layers with an ambient temperature increasing with depth (116°C/km). We assumed that the potential temperature of the mantle was around 1550°C at 3.8 Ga compared to 1330°C today. Formulas from Flament et al. [2008] (with an activation energy for the mantle viscosity of ~350-400 kJ/mol) lead to a oceanic crust mass production rate ~11 times higher at 3.8 Ga than today. This value is consistent with the estimation by Ohta et al. [1996] of a 10 times higher oceanic crust production rate during the Archean.

We assumed that the CO₂ flux outgassed by oceanic ridges is proportional to the oceanic crust mass production rate while the seafloor weathering, limited to the first km of the oceanic crust, is proportional to the spreading rate. For a fixed oceanic crust mass production rate, a thicker crust has no impact on the outgassed CO₂ flux but reduces the seafloor weathering in our model, leading to an enhancement of pCO₂. It is generally accepted that the oceanic crust was thicker during the Archean (~18 km versus 6 km in average) [Bickle, 1986; Ohta et al., 1996; Flament et al., 2008]. If we assume that the early oceanic crust was 3 times thicker at the mid-oceanic ridge, the spreading rate (SR) must have been 3.8 times higher than today, according to our estimation of the oceanic crust mass production rate. Since the oceanic crust thickness and the associated spreading rate during the Archean are highly uncertain [Bickle, 1986], we therefore tested cases with the present-day oceanic crust thickness (6 km) and SR=11.3.

For a warmer mantle and a higher spreading rate on the early Earth, the younger oceanic lithosphere was warmer and melted more easily when subducted. If arc volcanoes, associated with subduction zones, provide a quick return of materials to the surface, the direct melting of oceanic crust may increase the amount of fluid released, suggesting a recycling of suducted carbon higher than today. According to this assumption, the decarbonation efficiency may exceed Cenozoic values used as reference [Hoareau et al., 2015]. We thus explored the effect of the recycling (R) using the Cenozoic value as reference (R=0.1, corresponding to a recycling efficiency of 10%) and values 3 times higher (R=0.3, a recycling efficiency of 30%).

The fraction of emerged land during the Hadean/Archean is still debated [Roberts and Spencer, 2015]. We made the assumption that this fraction was low (<30% of the present-day fraction) during the early Archean as suggested by some geochemical data (e.g. [Dhuime et al., 2012]), some evidence for a more voluminous ocean during the Archean [Pope et al., 2012] and some models (e.g. [Flament et al., 2008]). This hypothesis allows us to use the outputs from the 3D land-free simulations for the carbon cycle model. We explored the effect of the emerged continental fraction considering four cases:

* 1) an aquaplanet (no land)

* 2) a very low emerged land fraction (1% of the present-day land fraction)

* 3) a moderate continental fraction (12% of the present-day land fraction)

* 4) a relatively high continental fraction (26% of the present-day fraction)

Although there is clear evidence of emerged land at 3.7-3.8 Ga [Nutman et al., 2015], an aquaplanet model is an idealized case, which is very useful for understanding the effect of seafloor weathering.

## 3.2. Results and discussion

We ran the climate-carbon model for the four land fractions (0, 0.01, 0.12 and 0.26 times present-day fraction), with a recycling factor R equal to 1 (present-day) or 3, with or without 2 mbars of atmospheric methane and with the present-day oceanic crust thickness (H=6 km, SR=11.3) or a 3 times thicker oceanic crust (H=18 km, SR=3.8). The equilibrium points span a temperature range from 8.5°C to 31.5°C and a pCO₂ range from 0.12 bar to 0.36 bar (see figure 4 and table 1). Ice-free or partly ice-free oceans are obtained in any cases.

### 3.2.1. Cases with no land

If there is no emerged land, the carbon cycle is controlled by the seafloor weathering and volcanic/mid-oceanic ridge outgassing. At equilibrium, the seafloor weathering is simply given by $F_{sfw} = \frac{F_{mor}+F_{deep}}{1-R}$ ($F_{mor}$ and $F_{deep}$ are the CO₂ fluxes from midocean ridges and mantle, respectively, see Appendix A) and depends on R. Because the carbon species partitioning between ocean and atmosphere is controlled by temperature, sea surface temperature (SST) implicitly determines the seawater pH and pCO₂ for a fixed value of seafloor weathering rate. Light blue and purple curves in figure 4 reveal the CO₂ partial pressure for a given SST with our model. They are computed
assuming no CO₂ climate feedback and thus represent the intrinsic response of the carbon model for Earth's climate in the absence of continents. The light blue curves (6 km thick crust) reveal a maximum pCO₂ of 0.14 bar for R=0.1 (left panel) and 0.26 bar for R=0.3 (right panel). The CO₂ solubility in the ocean decreases with temperature. At temperatures lower than 15°C, this feedback dominates and pCO₂ increases with temperature. When the temperature becomes higher than 15°C, the temperature dependence of the seafloor weathering dominates and pCO₂ decreases with temperature. The combination of both effects (solubility and seafloor weathering) leads to an upper limit for pCO₂ for the thin oceanic crust case. The reduced seafloor accretion rate associated to a thicker oceanic crust scenario combined with cold water strongly limits carbon consumption by seafloor weathering, hence the purple curves (18 km thick crust) do not show a maximum pCO₂ below 1 bar. For these cases without land, the equilibrium point of the climate-carbon model corresponds to the crossing between the intrinsic response of the carbon model (i.e. the light blue and purple curves) and the GCM curves showing the global mean surface temperature as a function of pCO₂ (i.e. the blue and red curves). For instance, the equilibrium point for a methane-free atmosphere with R=0.1 and a thick crust corresponds to a mean SST of 22.3°C and pCO₂=304 mbar.

### 3.2.2. Effects of land cover

According to our model, land cover has only a limited impact on the carbon cycle and the global climate for the early Earth (temperature change of ~1°C between the four land distributions). Indeed, while the carbon cycle was mostly controlled by continental weathering during the Proterozoic and the Phanerozoic (e.g. *Krissansen-Totton and Catling* [2017]), it was mostly controlled by seafloor weathering on the early Earth. That difference is mainly due to the higher spreading rate (high accretion rate provides a sufficient quantity of weatherable materials which enhances the seafloor weathering) and to the assumed lower land fraction on the early Earth. Under such circumstances, our assumption about the land distribution (i.e. small continents distributed homogeneously) is therefore of minor importance for the results. We also notice that the feedback produced by land does not have the same sign with a 6 or 18 km thick oceanic crust. With a 6 km oceanic crust, the seafloor weathering is slightly stronger than the continental weathering per unit of surface. In that case, the global mean surface temperature and pCO₂ are higher with a higher land fraction. In contrast, with a 18 km thick oceanic crust, the continental weathering is stronger than the seafloor weathering per unit of surface (as today) and a higher land fraction leads to lower pCO₂. The coldest case with our assumptions (mean SST of 8.5°C) is thus obtained with no land, a methane-free atmosphere, R=0.1 and a thin crust. The warmest case (mean SST of 31.5°C) is obtained with no land, a methane-rich atmosphere, R=0.3 and a thick oceanic crust, as shown by the blue point in Figure 4b.

### 3.2.3. Effects of methane, oceanic crust thickness and recycling rate

A thicker oceanic crust or a higher recycling efficiency always leads to an increase of pCO₂ and so to an increase of temperature. In contrast, adding methane in the atmosphere reduces pCO₂ but still increases the global mean surface temperature. If we use as a reference, the case with a moderate land cover, a 6 km thick crust, no methane and R=0.1 (mean surface temperature=8.5°C), adding 2 mbar of methane produces a 11.7°C warming, a 18 km crust produces a 12.9°C warming, and a 3 times higher recycling rate produces a 8.8°C warming. These three parameters therefore produce significant temperature and pCO₂ changes. In particular, the thicker oceanic crust is generally expected for the early Earth, with an effect on pCO₂. A thicker oceanic crust would have had a large impact, allowing warmer climate for the early Earth than today. The combination of a thicker crust, a high recycling rate, and the presence of methane leads to global mean surface temperatures around 30°C for the three land distributions. That temperature change (~20°C) is much lower than the sum of individual warming (~33°C), mostly because of the non-linear dependence of the seafloor weathering with the deep ocean temperature.

### 3.2.4. Absence of plate tectonic: the heat-pipe scenario

The onset time of plate tectonic is unknown but most geological and isotopic data suggest that subduction started during the Hadean or early Archean [*Sleep*, 2007; *Tang et al.*, 2016]. *Moore and Webb* [2013] suggested that before the onset of plate tectonic, Earth was in a heat-pipe regime in which volcanism dominated surface heat transport. We tested the consequences of this possible early regime for the carbon cycle with our model. We considered that all the heat flux is released by volcanism, that the spreading rate is null and that submarine basalt flows are weathered as seafloor basalts. In our simulation of the heat-pipe Earth, we assumed a very low land fraction (1% of the present-day land cover). The volcanic outgassing becomes around 4 times as high as with the plate tectonic regime, or 10 times as
<table><caption>Table 2: Equilibrium states from the carbon model coupled to the GCM. Columns from left to right are: the land fraction, the oceanic crust thickness (in km), the spreading rate, the recycling rate, the presence (2 mbar) or absence of methane, the outgassing rate, pCO₂ (in mbar), the mean surface temperature (in °C), the mean oceanic pH, the mean dissolved inorganic carbon (in mol/m³), the CO₂ flux outgassed from mid-oceanic ridges and from arc volcanoes, and the CO₂ sink by continental weathering and seafloor weathering (all fluxes in 1e12 mol/yrs). The first case, with land fraction equal to 1, is the present-day Earth. The last case, with SR=0, is the heat-pipe scenario with no plate tectonics.</caption><thead><tr><th colspan="6">Boundary conditions</th><th colspan="7">Results</th></tr><tr><th>Land fraction</th><th>Oceanic crust (km)</th><th>SR</th><th>Recycling</th><th>CH<sub>4</sub></th><th>rdeg</th><th>pCO<sub>2</sub> (mbar)</th><th>Tsurf (°C)</th><th>pH</th><th>DIC</th><th>F<sub>mor</sub></th><th>F<sub>arc</sub></th><th>F<sub>con</sub></th><th>F<sub>sfw</sub></th></tr></thead><tbody><tr><td>1</td><td>6</td><td>1</td><td>0.1</td><td>no</td><td>1</td><td>0.285</td><td>14.9</td><td>8.4</td><td>2.4</td><td>1.6</td><td>7.7</td><td>-10.5</td><td>-1.1</td></tr><tr><td>0.0</td><td>6</td><td>11.3</td><td>0.1</td><td>no</td><td>2.5</td><td>118</td><td>8.5</td><td>6.70</td><td>32.2</td><td>18</td><td>11.2</td><td>0.0</td><td>-29.2</td></tr><tr><td>0.0</td><td>6</td><td>11.3</td><td>0.3</td><td>no</td><td>2.5</td><td>225</td><td>16.5</td><td>6.60</td><td>39.8</td><td>18</td><td>19.5</td><td>0.0</td><td>-37.5</td></tr><tr><td>0.0</td><td>18</td><td>3.8</td><td>0.1</td><td>no</td><td>2.5</td><td>304</td><td>22.3</td><td>6.56</td><td>48.9</td><td>18</td><td>11.2</td><td>0.0</td><td>-29.2</td></tr><tr><td>0.0</td><td>18</td><td>3.8</td><td>0.1</td><td>yes</td><td>2.5</td><td>190</td><td>27.0</td><td>6.66</td><td>35.4</td><td>18</td><td>11.2</td><td>0.0</td><td>-29.2</td></tr><tr><td>0.0</td><td>18</td><td>3.8</td><td>0.3</td><td>no</td><td>2.5</td><td>359</td><td>25.8</td><td>6.53</td><td>52.4</td><td>18</td><td>19.5</td><td>0.0</td><td>-37.5</td></tr><tr><td>0.0</td><td>18</td><td>3.8</td><td>0.3</td><td>yes</td><td>2.5</td><td>243</td><td>30.6</td><td>6.61</td><td>39.6</td><td>18</td><td>19.5</td><td>0.0</td><td>-37.5</td></tr><tr><td>0.01</td><td>6</td><td>11.3</td><td>0.1</td><td>no</td><td>2.5</td><td>118</td><td>8.5</td><td>6.70</td><td>32.21</td><td>18</td><td>11.2</td><td>-0.05</td><td>-29.2</td></tr><tr><td>0.01</td><td>6</td><td>11.3</td><td>0.3</td><td>no</td><td>2.5</td><td>231</td><td>16.9</td><td>6.65</td><td>45.10</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.01</td><td>18</td><td>3.8</td><td>0.1</td><td>no</td><td>2.5</td><td>303</td><td>22.3</td><td>6.56</td><td>49.51</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.01</td><td>18</td><td>3.8</td><td>0.1</td><td>yes</td><td>2.5</td><td>186</td><td>26.8</td><td>6.67</td><td>35.39</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.01</td><td>18</td><td>3.8</td><td>0.3</td><td>no</td><td>2.5</td><td>358</td><td>25.8</td><td>6.53</td><td>52.73</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.01</td><td>18</td><td>3.8</td><td>0.3</td><td>yes</td><td>2.5</td><td>240</td><td>30.5</td><td>6.62</td><td>39.64</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>6</td><td>11.3</td><td>0.1</td><td>no</td><td>2.5</td><td>117</td><td>8.5</td><td>6.70</td><td>32.06</td><td>18</td><td>11.6</td><td>-0.4</td><td>-29.2</td></tr><tr><td>0.12</td><td>6</td><td>11.3</td><td>0.1</td><td>yes</td><td>2.5</td><td>97</td><td>20.2</td><td>6.80</td><td>25.52</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>6</td><td>11.3</td><td>0.3</td><td>no</td><td>2.5</td><td>235</td><td>17.3</td><td>6.66</td><td>46.69</td><td>18</td><td>23.7</td><td>-3.1</td><td>-38.7</td></tr><tr><td>0.12</td><td>6</td><td>11.3</td><td>0.3</td><td>yes</td><td>2.5</td><td>135</td><td>23.3</td><td>6.75</td><td>31.47</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>18</td><td>3.8</td><td>0.1</td><td>no</td><td>2.5</td><td>292</td><td>21.4</td><td>6.59</td><td>50.13</td><td>18</td><td>12.6</td><td>-4.1</td><td>-26.5</td></tr><tr><td>0.12</td><td>18</td><td>3.8</td><td>0.1</td><td>yes</td><td>2.5</td><td>163</td><td>25.2</td><td>6.71</td><td>34.38</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>18</td><td>3.8</td><td>0.3</td><td>no</td><td>2.5</td><td>356</td><td>25.6</td><td>6.57</td><td>56.31</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>18</td><td>3.8</td><td>0.3</td><td>yes</td><td>2.5</td><td>219</td><td>29.0</td><td>6.69</td><td>41.89</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.26</td><td>6</td><td>11.3</td><td>0.1</td><td>no</td><td>2.5</td><td>121</td><td>8.8</td><td>6.70</td><td>32.50</td><td>18</td><td>12.1</td><td>-0.64</td><td>-29.5</td></tr><tr><td>0.26</td><td>6</td><td>11.3</td><td>0.3</td><td>no</td><td>2.5</td><td>245</td><td>18.0</td><td>6.66</td><td>48.47</td><td>18</td><td></td><td></td><td></td></tr><tr><td>0.12</td><td>6</td><td>0</td><td>0</td><td>no</td><td>10</td><td>257</td><td>18.8</td><td>6.60</td><td>45.75</td><td>31.9</td><td>1.1</td><td>-0.9</td><td>-32.1</td></tr></tbody></table>

high as the present volcanic outgassing. At equilibrium, pCO₂ is 0.26 bar and the mean surface temperature is around 18.8°C. A heat-pipe regime would thus have led to a temperate early Earth.

### 3.2.5. Discussion

Our climate-carbon model never reaches warm climates for the early Earth because of the strong temperature feedback by the seafloor weathering. It predicts global mean surface temperatures around 8.5-30.5°C, with pCO₂ around 0.1-0.36 bar. This temperature range is compatible with the most recent measurements by Hren et al. [2009] and Blake et al. [2010], and also with the evolutionary model from Bousseau et al. [2008] for a mesophilic last universal common ancestor (LUCA). Thermophily near the roots of the tree of life may reflect local warm environments (like hydrothermal vents) or survival after hot climates produced just after large impacts [Bousseau et al., 2008; Abramov and Mojzsis, 2009]. The pCO₂ range of 0.1-0.35 bar is generally higher than the upper limits initially derived from the stability of different minerals from Archean samples [Rye et al., 1995; Sheldon, 2006; Rosing et al., 2010; Driese et al., 2011]. These constraints globally suggest a pCO₂ lower than 30 mbars. However, they are not all compatible together and were essentially obtained for the late Archean. Recent work by Kanzaki and Murakami [2015] gives an upper limit for pCO₂ of 0.14-0.7 bar at 2.77-2.46 Ga, compatible with our results.

## 4. Climate and carbon cycle during the Late Heavy Bombardment

### 4.1. Impact modeling

The Earth was likely struck by numerous of impactors during the Late Heavy Bombardment (LHB), which would have occurred between around 4.1 and 3.8 Ga [Gomes et al., 2005; Bottke and Norman, 2017]. These impactors may have been a sink of CO₂ by producing large amounts of ejecta, which can be easily weathered [Sleep and Zahnle, 2001]. They may also have been a source of CO₂ by vaporizing carbonates from the oceanic crust. For instance,
<figure><img src="image_119.png" /></figure>

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

Figure 4: Global mean surface temperature as figure 1 but including the equilibrium points of the carbon model assuming a land fraction of 0.01 (light blue and blue circles) and 0.12 (yellow and red circles) compared to the present-day land fraction. Panel a) shows equilibrium points for the present-day recycling (R=0.1). Light blue and yellow circles are for the present-day oceanic crust thickness while blue and red circles for a 3 times thicker oceanic crust. The black diamond corresponds to the heat pipe scenario with no plate tectonic on an aquaplanet. Panel b) is similar to panel a) but with a 3 times higher recycling (R=0.3). The light blue and the purple curves are the intrinsic response of the carbon cycle model. They are obtained by changing the surface temperature from 10°C to 40°C with no feedback of pCO₂ on the temperature, with no land, for the present-day oceanic crust (light blue) and for a 3 times thicker oceanic crust (purple).

Chixulub impact may have outgassed 350-3500 Gt of CO₂ [Pierazzo et al., 1998], or the equivalent of 0.06-0.6 mbar of CO₂. We investigated the long term effects of impacts on the carbon cycle to understand their consequences for the mean climate during the LHB, in particular if they had a global net warming or cooling effect.

We assumed that the frequency of impact follows a scaling law given in Anbar et al. [2001] and Sleep and Zahnle [2001]:
$$f(\ge m) = f_0 \times \left(\frac{m}{m_{max}}\right)^{-b} \qquad (1)$$

where f is the frequency of impactors having a mass higher than m. For the LHB, we chose b = 0.7 [Anbar et al., 2001]. We considered that the maximal mass was $m_{max}$ = 10¹⁹ kg (i.e. diameter around 190 km). This upper limit is constrained by the largest crater observed on the Moon. We neglected impacts with a mass lower than 10¹³kg (i.e. diameter lower than around 1.8 km), which produce little ejecta because of the strong deceleration by the ocean. We fixed $f_0 = 2.9 \times 10^{-8}$ impacts/year in order to have a cumulative impactor mass equal to $2 \times 10^{20}$ kg for the 300 millions of years of the LHB. That corresponds to an impact flux of $6.7 \times 10^{11}$ kg/year.

As Sleep and Zahnle [2001] and Zahnle and Sleep [2002], we computed the weathering of ejecta, which provides a source of cations, precipitating as carbonates in the ocean. We assumed that the mass distribution of ejecta is given by $N(\ge m) \propto m^{-\gamma}$, where $\gamma=0.87$ (derived from scaling laws of Collins et al. [2005]). Zahnle and Sleep [2002] used a similar law with $\gamma=0.9-0.95$. We considered that the diameter of each particle decreases of 1 mm in 167 kyrs [Crovisier et al., 1987] and that cations are present at 0.005 mol/g. With the oceanic crust lifetime for the early Earth, ejecta with a diameter larger than 100 m contribute little to the weathering compared to smaller ejecta.

We included CO₂ delivered by the impactors as done by de Niem et al. [2012]. We assumed that the mass fraction of CO₂ is 1% in asteroids and 20% in comets. For comparison, Rosetta measurements of Comet 67P / Churyumov-Gerasimenko indicate a mass ratio of CO₂ to H₂O around 40% [Marty et al., 2016] for a mass fraction of H₂O to dust around 20%, implying a mass fraction of CO₂ around 8%. If all impactors were comets, the LHB would have provided the equivalent of 3.8 bar of CO₂ (taking into account the escaping mass of impactors). The fraction of comets during the LHB is not known and Bottke et al. [2012] suggest that comets were a minor player. We did simulations using fractions of comets of 0 and 100%. We also included the CO₂ outgassed by oceanic carbonates and lithosphere melting after impact. We considered that the oceanic crust is covered by a 40 m layer of carbonate, which outgassed CO₂ when shocked at a high pressure (above 60 GPa). We assumed that the lithosphere below the carbonate layer (considered as pure calcite) contains 100 ppmv of CO₂ and outgasses all CO₂ when melted. The complete assumptions for the ejecta weathering and CO₂ outgassing are described in Appendix B.

Figure 5 shows the amount of CO₂ outgassed from impactor/crust and consumed by ejecta as a function of impactor mass for asteroids (density=3000 kg/m² and initial velocity=21 km/s) and comets (density=1000 kg/m² and initial velocity=30 km/s). The solid red lines show the maximal amount of CO₂ that can be consumed by ejecta. The solid red line show the amount of CO₂ consumed by ejecta in 15 Ma (the lifetime of the oceanic crust in our model). With our size distribution for ejecta, most of ejecta are weathered during the lifetime of the oceanic crust. In all cases, impacts produce a net sink of CO₂ by ejecta weathering.

## 4.2. Results and discussion for the LHB

To simulate the effect of impacts during the LHB, we computed sequences of 10 millions of years with stochastic impacts following the statistical distribution described above. For simplicity, we assumed that impacts never overlap and that they all fall into the ocean. This assumption is valid because we consider small continental fractions and because the total surface of transient craters (diameter $D_{tc}$ in Appendix B) accumulated in 10 million years only represents 0.6% of the Earth's surface. These sequences provide a time-dependent CO₂ outgassing rate and alkaline source, which are used in the carbon cycle model to computed the climate and pCO₂ evolution. The top panel in figure 6 shows the cumulated CO₂ outgassed by impacts and consumed by ejecta during a typical 10 million year sequence with stochastic impacts. In this sequence, one 143 km impactor strikes the Earth at 7 Ma and dominates CO₂ fluxes. On average, the weathering of ejecta produces a cations flux of $16 \times 10^{12}$ mol/yr (compared to $300 \times 10^{12}$ mol/yr in Sleep and Zahnle [2001]), around half the CO₂ flux from mid oceanic ridges and arc volcanoes (around $30 \times 10^{12}$mol/year). The vaporization of carbonates produces a CO₂ flux of $0.05 \times 10^{12}$mol/year, negligible on the long term compared to the other fluxes. The bottom panel in figure 6 shows the evolution of pCO₂ during this sequence for the 0.01 land fraction, no methane, R=0.1, a thin or thick oceanic crust and with asteroids or comets. The giant impact at 7 Ma produces a strong decrease of pCO₂ by ejecta weathering triggering a full glaciation for thin oceanic crust cases and a very cold climate with a global mean surface temperature of -10°C for the thick oceanic crust case. Before this impact, the carbon cycle is stabilized with pCO₂ around 50 mbar for the thin crust and around 250 mbar
for the thick crust. As in Sleep and Zahnle [2001] and Zahnle and Sleep [2002], we thus find a strong cooling by ejecta weathering leading to cold climates, but not always completely ice covered. In particular, cases with a thick crust, methane and R=0.3 should remain temperate or at least partially ice-free most of the time during the LHB.

To better understand the role of each parameter, we derived a simple analytical expression of pCO₂. If we assume a case with no land, a recycling R₀ and an oceanic crust thickness H₀, the equilibrium of the carbon cycle can be expressed for a reference state with no impact as:

$$F_{mor}^{0} + F_{arc}^{0} + F_{sfw}^{0} = 0 \qquad (2)$$

where F₀, Farc and Fsfw are the CO₂ fluxes from mid-oceanic ridges, arc volcanoes and seafloor weathering from the (crust+mantle) to the (ocean+atmosphere). For a new state with other values of R and spreading rate and with impacts, the equilibrium is given by:

$$F_{mor} + F_{arc} + F_{sfw} + F_{ej} = 0 \qquad (3)$$

with F₀ = F₀ and Farc = Fdeep - R × (Fsfw + Fej) (see appendix A). We make the assumption that the seafloor weathering does not depend on oceanic temperature, what is valid for cold climates. By using the simple expression:

$$F_{sfw} = F_{sfw}^{0} \left( \frac{pCO_2}{pCO_2^{0}} \right)^{\alpha} \left( \frac{SR}{SR_0} \right) \text{ with } \alpha=0.5 \text{ [Haqq-Misra et al., 2016], we obtain an expression for pCO}_2\text{:}$$

$$pCO_2 = pCO_2^0 \times \left(\frac{SR_0}{SR}\right)^{1/\alpha} \left(\frac{1-R_0}{1-R}\right)^{1/\alpha} \left(1 - \frac{(1-R) \times F_{ej}}{F_{mor} + F_{deep}}\right)^{1/\alpha} \qquad (4)$$

The term $(1-R) \times F_{ej} / F_{mor} + F_{deep}$ simply corresponds to the ratio of the amount of CO₂ lost in the mantle (due to ejecta weathering) by the amount of CO₂ outgassed directly from the mantle (e.g. without taking into account recycling). When this ratio is higher than 1, there is a net sink of CO₂ toward the mantle and the Earth necessary falls into a full snowball Earth. That occurs for $F_{ej} \le -29 \times 12 \text{ mol/yr}$ for R=0.1 and $F_{ej} \le -37 \times 12 \text{ mol/yr}$ for R= 0.3 with our model values. With our impact distribution, the ejecta cations flux ($16 \times 10^{12} \text{ mol/yr}$) is below these limits, and the Earth should have avoided a permanent full glaciation during the LHB. Because of the large uncertainties on the carbon cycle parameters and the duration and intensity of the LHB, it is not possible to be definitive on that point. We conclude that the LHB likely had a strong impact on the carbon cycle and climate with potentially snowball Earth events after large impacts.

Large impactors were also frequent during all the Archean [Bottke et al., 2012; Johnson and Melosh, 2012; Bottke and Norman, 2017] but much less than during the LHB. According to the dynamical model from Bottke et al. [2012], the impact frequency was an order of magnitude lower at 3.7 Ga than at the peak of the LHB. Therefore, impacts probably had only a limited effect on the long term carbon cycle and climate during most of the Archean. However, our LHB impact sequence reveals that individual giant impacts could trigger episodic glaciations. Johnson and Melosh [2012] evaluated the diameters of ancient large impactors from spherule layers during the Archean and the Proterozoic. A few of them with diameters higher than around 50 km could have been large enough to trigger glaciations.

The existence of the Late Heavy Bombardment is still debated by some authors. For instance, Boehnke and Harrison [2016] suggest that the LHB could be an artifact of monotonically decreasing impact flux combined with episodic crust formation. If it did not occur, the impact flux at the end of the Hadean would have been too small to affect the climate on the long term, but strong cooling after episodic large impacts may still have been possible.

## 5. Conclusions

Our 3D climate simulations (taken in isolation) reveal that warm climates, with mean surface temperatures around 60°C, could have been obtained during the late Hadean/early Archean with around 1 bar of CO₂. However, our carbon model suggests that such states would not have been stable. Indeed, we found that warm climates can be maintained
<figure><img src="image_121.png" /><figcaption>Figure 5: Effect of an individual impact on atmospheric CO₂. Amount of CO₂ outgassed from impactor (black solid lines) and from the crust (black dashed lines) or consumed by ejecta weathering as a function of impactor mass for asteroids (top panel) and comets (bottom panel). The red solid line is the maximal amount of CO₂ that can be consumed by ejecta weathering, while the red dashed line is the amount of CO₂ consumed during 15 Ma (approximately the lifetime of the Archean oceanic crust).</figcaption></figure>

only if the seafloor weathering is assumed independent of temperature. Such an assumption is inconsistent with recent measurements indicating a higher weathering during the late Mesozoic, when the deep ocean was around 10°C warmer [Coogan and Gillis, 2013]. Stable very warm climates with a N₂-CO₂-CH₄-H₂O early atmosphere are therefore excluded.

Our results favor instead temperate climates with global mean surface temperatures around 8.5-30.5°C and with pCO₂ around 0.1-0.36 bar. With the increasing solar luminosity over the Archean, the carbon cycle would have led to a reduction of pCO₂. In parallel, with the emergence of land and the slowing seafloor spreading rate during the Archean, the carbon cycle would have progressively become controlled by continental weathering while initially controlled by seafloor weathering. However, a better understanding of the data from O isotopes in ancient cherts together with better constraints on pCO₂ are required to validate our conclusions about the long-debated climate of early Earth.

During the Late Heavy Bombardment, the weathering of ejecta would have dramatically decreased the partial pressure of CO₂ leading to cold climates but not necessary to a snowball Earth during all that period. We predict that large impactors (i.e. with diameters higher than around 50 km) could have trigger extreme glaciations during the Hadean and the Archean.
<figure><img src="image_122.png" /><figcaption>Figure 6: Effect of impacts on pCO₂. Top panel: impact scenario expressed in cumulative mol of CO₂ outgassed by impacts (black line) and consumed by ejecta (red line) with asteroids (solid lines) or comets (dotted lines). Bottom panel: evolution of pCO₂ for the impact scenario, assuming a low land fraction, R=0.1 and no methane. Yellow line is a case with a thin oceanic crust, asteroids and ejecta weathering. Black line is as yellow line but without ejecta weathering. Blue line is as yellow line but with comets. Red line is as yellow line but with a thick oceanic crust.</figcaption></figure>

Finally, our modeling work suggests that the CO₂-dependent carbon cycle could have maintained temperate climates on the early Earth, without requiring any additional greenhouse gas or warming process. This work also highlights the importance of the seafloor weathering, which can dominate the continental weathering. Under such circumstances, the climate becomes relatively insensitive to the presence of land. The strong temperature dependence expected for the seafloor weathering also has implications for the habitability of exoplanets, allowing in particular a stabilizing climate-weathering feedback on waterworlds. Waterworlds might therefore not be much more subject to water loss by moist greenhouse than partially ocean-covered Earth-like planets, as suggested by Abbot et al. [2012].

*Acknowledgments*: B.C. acknowledges support from Paris Sciences et Lettres and from an appointment to the NASA Postdoctoral Program, administered by USRA. This work was performed as part of the NASA Astrobiology Institute's Virtual Planetary Laboratory, supported by NASA under Cooperative Agreement No. NNA13AA93A. DC also acknowledges support from NASA Exobiology Program grant NNX15AL23G. We thank J. Kasting and an anonymous reviewer for helpful comments.
Appendix A. Carbon cycle modeling

<figure><img src="image_123.png" /><figcaption>Figure A.7: Representation of the carbon model used for the Hadean Earth.</figcaption></figure>

The carbon cycle model has been previously used for the Proterozoic and the Phanerozoic [Goddéris et al., 2006; Le Hir et al., 2008a, b]. Figure A.7 shows the different processes included in the model. The parametrizations for the different CO₂ flux are briefly described below.

CO₂ outgassing by oceanic ridges and arc volcanoes

We considered that the CO₂ flux from mid oceanic ridges F<sub>mor</sub> is proportional to the production of oceanic crust:

$$F_{\text{mor}} = F_{\text{mor}_0} \times \frac{SR}{SR_0} \times \frac{H}{H_0} \qquad (\text{A.1})$$

where F<sub>mor<sub>0</sub></sub>= 1.8×10¹² mol/yr, SR is the spreading rate and H is the oceanic crust thickness (H<sub>0</sub> = 6 km).

We assumed that the flux from arc volcanoes (F<sub>arc</sub>) is the total degassing flux due to carbonate recycling at subduction zones and mantle outgassing (with the exception of mid-oceanic ridge). Degassing over geologic times being assumed proportional to the heat flux from the interior of the Earth, the mantle outgassing rate is adjusted, as appropriate, using a forcing factor, rdeg=2.5 at 3.8 Ga [Goddéris and Veizer, 2000]):

$$F_{\text{arc}} = F_{\text{sub}} \times R + F_{\text{deep}} \qquad (\text{A.2})$$

where R is the recycling rate and F<sub>deep</sub> = C<sup>0</sup><sub>deep</sub> × rdeg with C<sup>0</sup><sub>deep</sub> = 3.3 × 10¹² mol/yr. For the heat-pipe case, we