<table><caption>**Table 1.** Simulations performed in the Sect. 3 ($\Delta F_{\text{trop}}$ is held constant and equals to +19 W m$^{-2}$).</caption><thead><tr><th rowspan="2">set of simulations</th><th colspan="2">Boundary conditions</th><th rowspan="2">LOD (hours)</th><th rowspan="2">Droplet size(µm)</th><th rowspan="2">seawater freezing point</th></tr><tr><th><i>Paleogeography</i></th><th><i>Solar constant</i></th></tr></thead><tbody><tr><td>modern clouds</td><td>3.5-1Ga (Fig. 2b)</td><td>3.5-1Ga (Fig. 2a)</td><td>24</td><td>5-10</td><td>-1.8 °C</td></tr><tr><td>no clouds</td><td>“</td><td>“</td><td>24</td><td>/</td><td>“</td></tr><tr><td>Large droplet clouds</td><td>“</td><td>“</td><td>24</td><td>20</td><td>“</td></tr></tbody></table>

<figure><img src="image_60.png" /><figcaption>**Fig. 2.** Boundary conditions for solar constant (upper part) and continental growth (lower part). $\eta_{\text{land}}$ is the emerged surface/Earth surface ratio (the present-day surface is $146 \times 10^6$ km$^2$ or 29% of the Earth's surface). Each red square represents the boundary conditions used in our sets of simulations.</figcaption></figure>

(Godderis and Veizer, 2000) and $\delta^{34}$S throughout the Archaean and early Proterozoic. Accordingly, the surface of the emerged continents increased from 4% at 3.5 Ga to 27.5% at 1 Ga. Palaeogeographies were reconstructed based on paleomagnetic data from 1 to 2.5 Ga (Pesonen et al., 2003). 3.5 to 2.75 Ga plate reconstructions being poorly constrained, continents were assumed to be located at low and mid latitudes. The elevation of land points was fixed at 200 m in all experiments.

### 2.2.2 Radiative forcing by greenhouse gases and ozone amounts

To set up boundary conditions correctly, we estimated the greenhouse radiative forcing provided by the combination of 0.9 mbar of CO$_2$ and CH$_4$ (i.e. CO$_2$ and CH$_4$ concentrations used by Rosing et al. (2010) and also tested by Charney et al. (2013)), rather than directly use $pCO_2$ and $pCH_4$ as boundary conditions. This was motivated because the assumption of such a low CO$_2$/high CH$_4$ is likely incorrect (Dauphas and Kasting, 2011; Reinhard and Planavsky, 2011). Moreover, Rosing et al. (2010) ignored the organic haze formation ($[CH_4]/[CO_2] > 0.2$ (Haqq-Misra et al., 2008)). Since it makes little difference whether radiative forcing is controlled by CO$_2$ or CH$_4$ (Hansen et al., 1997), this approach avoids doubtful assumptions about the Archaean/Proterozoic atmosphere composition. Hence, using the present-day atmosphere as a reference, the change in net radiative flux at the tropopause ($\Delta F_{\text{trop}}$, Fig. 1) equals $19 \text{ W m}^{-2}$ ($\Delta F_{\text{trop}} = 6$ and $13 \text{ W m}^{-2}$ for 0.9 mbar of CO$_2$ and CH$_4$, respectively).

In addition, due to the low oxygen content of the Archean atmosphere imposed by sulfur mass independent fractionation (S-MIF) of sedimentary sulfides and sulfates ($10^{-5}$ Present Atmospheric Level, PAL, Farquhar et al., 2001; Pavlov et al., 2000), the stratospheric ozone layer was removed in all experiments prior to 2.25 Ga and constrained to the present-day value in other experiments. The absence of an ozone layer precluded the formation of a stratospheric temperature inversion and resulted in a 40 km-thick atmosphere without temperature stratification (corresponding to the topmost atmospheric grid points in the AGCM). The global mean lapse rate from the surface to the 200 hPa level was not significantly affected by the absence of the ozone layer, and only a slight mean global cooling (0.7 °C) was simulated at the Earth's surface (Appendix Fig. A1). Induced cooling of surface temperature (1 to 2 °C) associated with the removal of the ozone layer is in good agreement with previous studies (Kasting, 1984; Kiehl and Dickinson, 1987; Wolf and Toon, 2013; Kunze et al., 2014).
### 2.2.3 Experimental design and assumptions for clouds

We performed sets of numerical experiments ranging from 3.5 to 1 Ga ago with a 0.25 Ga time increment. To avoid any effect of initial boundaries, simulations were initiated with warm conditions (present-day as starting point) and driven toward colder conditions (step-like decrease in the Sun's brightness). Each run lasted several tens of years until a steady state was reached for a decade or more. Climatic variables are averaged over the last decade.

To test the radiative impact of clouds, we considered three scenarios involving modern clouds, no cloud at all (cloud fraction fixed to 0), and clouds with large droplets (Table 1). A change in cloud droplet size affects scattering and therefore the cloud albedo (Pierrehumbert, 2010). For the parameterisation of Slingo (1989) implemented in FOAM, the cloud optical properties ($OP_{clds}$) for short-wave radiations were defined as

$$OP_{clds} = \text{CWP}[a + b/r_e], \quad (4)$$

where a and b are the cloud radiative coefficients (different for ice or liquid water clouds), CWP the cloud water path (g kg⁻<fnref n="1" />) and $r_e$ the effective droplet radius (droplet radius affected by its distribution).

In FOAM, liquid cloud droplet radii of modern clouds vary from 5 to 10 µm over landmasses and are fixed at 10 µm over oceans. This distinction is made to account for enhanced cloud drop nucleation that can occur over CCN-rich land areas. On the modern Earth, non-anthropogenic production of CCN is essentially related to soil dust formation, volcanic sulfate aerosols and dimethylsulfide (DMS) produced by planktonic eukaryotic algae (Charlson et al., 1987; Andree and Rosenfield, 2008). In the absence of planktonic algae during the Archaean, the scarcity of emerged terrestrial surfaces and subaerial volcanism were considered as the possible, yet speculative, causes for the low CCN availability in the early Earth atmosphere (Rosing et al., 2010). To quantify the influence of large liquid droplets on climate, we ran a set of experiments in which we used the mean radius of liquid droplets of 20 µm<fnref n="1" /> used by Rosing et al. (2010) over both the oceans and continents. Cloud droplets exceeding 15 µm being rarely observed (Bréon et al., 2002), this set of simulations should be regarded as an end member. Although large liquid water droplets of 20 µm diameter are expected to reduce cloud albedo (scattering effect), they also tend to rain out more rapidly (Boucher et al., 1995). According to the relationship established by Kump and Pollard (2008), the enhancement of cloud-to-precipitation rate varies as

$$\text{cloud-to-precipitation rate} = (r_e/r_{eo})^n, \quad (5)$$

where $r_e$ represents the new effective radius (20 µm), $r_{eo}$ the modern effective radius (assumed to be 10 µm, as oceans

cover a large fraction of the Earth's surface), and n the parameter of enhanced rain-out (here n = 3). The literature features a wide range of values for the enhanced rain-out parameter; here we used an intermediate value suggesting a mid-strength feedback between 1 and 5.37 (Goldblatt and Zahnle, 2011b; Penner et al., 2006; Charney et al., 2013). As the cloud lifetime cannot be changed in the absence of a cloud microphysics module, we increased the precipitable water by 8 to represent the enhanced rain-out. Albeit speculative (see Sect. 2.1), we showed that CWP fluctuations are only governed by the air moisture (Appendix Fig. A2).

## 3 Results: estimating pan-glaciation thresholds

Figure 3 shows the radiative deficit required to induce a pan-glaciation, the global mean surface temperature and the mean planetary albedo. It is worth emphasising that the results presented here are not representative of climate evolution over the Precambrian era because only part of the boundary conditions used are realistic (solar luminosity, palaeogeography, fraction of emerged land). Other boundary conditions including the greenhouse gas radiative forcing ($\Delta F_{\text{trop}} = +19 \text{ W m}^{-2}$) are unrealistic. However, these were implemented in our model to allow comparison with previous simulations (Charney et al., 2013; Rosing et al., 2010) and to determine the radiative forcing threshold values required for the onset of pan-glaciations. The experiments (Exp) are labelled from 3.5 to 1 Ga according to the boundary conditions used.

### 3.1 Influence of the solar constant

Figure 3b (blue curve) represents the global mean surface temperature (Ts) simulated using the present-day cloud parameterisation. Ts is as low as -75 °C in Exp(3.5 Ga) and increases linearly up to -65 °C in Exp(2.5 Ga). The planetary albedo always exceeds 0.65 (Fig. 3c, blue curve) when the Earth is ice covered. An abrupt change, expressed by a shift of both temperature and albedo, is observed between Exp(2.5 Ga) and Exp(2.25 Ga). In Exp(2.25 Ga), Ts reaches -20 °C and the albedo drops to 0.45. This corresponds to a partially frozen Earth where sea ice spreads up to 30° latitude in both hemispheres. Between Exp(2.25 Ga) and Exp(1.75 Ga), a slight global warming of about 5 °C is observed, followed by a marked change of some 24 °C between Exp(1.75 Ga) and Exp(1.25 Ga), and a final slight warming of 4 °C between Exp(1.25 Ga) and Exp(1 Ga).

From these experiments, three climate states can be defined: a snowball Earth between Exp(3.5 Ga) and Exp(2.5 Ga), a partially frozen Earth between Exp(1.75 Ga) and Exp(2.25 Ga), and a weakly frozen Earth between Exp(1 Ga) and Exp(1.25 Ga). Although the solar brightness increased sub-linearly between 3.5 and 1 Ga, the global mean surface temperature changed dramatically during this time

<footnote n="1">¹The effective radius of 20 µm used by Rosing et al. (2010) was finally revised to 17 µm (Goldblatt and Zahnle, 2011a).</footnote>
<figure><img src="image_61.png" /><figcaption>**Fig. 3.** 3-D climate simulations from 3.5 to 1 Ga with a greenhouse radiative forcing by CO₂ and CH₄ fixed to 19 W m⁻² (ΔF_trop; see Fig. 1). (a) Radiative deficit required for the onset of a pan-glaciation with RCM (no cloud runs, surface albedo fixed to 0.264), FOAM runs (* modern clouds and cloud-free runs and ** large-droplet clouds (20 µm)). Grey lines represent uncertainties related to temporal resolution (0.25 Ga). (b) Evolution of surface temperature versus time assuming increasing irradiance and linear continental growth (see Fig. 2). Grey line is from Rosing et al. (2010). Blue and red curves were obtained using liquid droplet sizes of 5 to 10 µm (modern clouds) and 20 µm (clouds with large droplets), respectively. Green squares represent the cloud-free scenario. (c) Evolution of planetary albedo versus time (α_ocean 0.06, α_rocky desert continents 0.32, α_seaice 0.5, α_snow 0.8). (d) Red stars represent sea ice extent (white) predicted with runs performed with large (20 µm) droplets (Fig. 3b, c shows the corresponding temperature/albedo).</figcaption></figure>

period (Fig. 3b). Between 1 and 1.25 Ga, the solar brightness is high enough to restrict the sea ice capping to the high latitudes (~60°) and to bring the equator-pole surface temperature gradient to the present-day level. Between Exp(1.25 Ga) and Exp(1.75 Ga), sea ice spreads towards the equator up to 30° of latitude due to the reduced incoming solar radiation. Between Exp(1.75 Ga) and Exp(2.25 Ga), the sea ice maximum extent remains located in the vicinity of the tropics. This is due to a latitudinal gradient of temperature two times higher than in the Exp(1 Ga), which arises from a strong meridional energy transport by the atmosphere and oceans. For a solar constant lower than 1148 W m⁻² (2.25 Ga), the solar brightness is low enough to drive climate towards snowball Earth conditions. Hence, for a fixed greenhouse-gas concentration of ΔF_trop = +19 W m⁻², the solar constant evolution and its interplay with the ice-albedo feedback appear as the predominant factor governing the Earth's climate. This contrast with the RCM results (Fig. 3a) where the albedo is held constant at 0.264 (Goldblatt and Zahnle, 2011b) due to the absence of the ice-albedo feedback.

## 3.2 Influence of continental surfaces

A lower surface albedo due to a reduced continental area ($α_{bare soil} = 0.32, α_{ocean} = 0.06$) has been suggested as a solution to the FYSP (Schatten and Endel, 1982; Rosing et al., 2010; but see Walker et al., 1982 and Goldbatt and Zahnle, 2011, for a controversy). No-cloud experiments show that the planetary albedo, which corresponds to the surface albedo, fluctuates with the growth of sea ice (Fig. 3b, green curve). Using the amount of energy reaching the Earth (expressed by $F = S/4(1-α)$, with $Δα = 0.10$ and $S = 1148$ W m⁻²), we estimated the radiative forcing induced by the ice-albedo feedback to -29 W m⁻² (cooling effect) between 1 and 2.25 Ga. This calculation shows that a slight change in the sea ice surface ($α_{sea-ice} = 0.6$) can easily overcome the influence of a smaller continental area on the surface albedo ($α_{bare soil} = 0.32, α_{ocean} = 0.06$). These results demonstrate that the surface albedo is mainly governed by the sea ice extension, the growth of continental surfaces being a second-order factor. This conclusion is to be challenged when the climate-carbon cycle feedbacks will be considered (Walker et al., 1981; Godderis and Veizer, 2000) (see companion paper).

## 3.3 Influence of clouds

Cloud experiments performed at 1 Ga show that the cloud radiative forcing is equal to -10 and -20 W m⁻² at high and mid-latitudes, respectively (Fig. 4b). Accordingly, the absence of clouds seems to limit the sea ice spreading at high latitudes and to prevent the increase in the planetary albedo (Fig. 3c). This mechanism can explain why cloud-free experiments result in a warmer climate than those using modern clouds (Fig. 3b, green curve). Surprisingly, the shift
<figure><img src="image_62.png" /><figcaption>Fig. 4. (a) Net cloud radiative forcing (CRF, unit W m⁻²) versus time. Blue squares represent runs with modern clouds. Large-droplet clouds (red squares) result in a positive forcing because the effect of infrared trapping overcomes that of scattering. (b) CRF as a function of latitude for simulations performed at 1 Ga. The black line represents large-droplet clouds and the dashed line the modern cloud parameterisation. (c) Simulations performed using 1 Ga boundary conditions. Dashed lines represent long-wave (LW) budget of clouds (greenhouse effect), continuous lines their short-wave budget (albedo effect).</figcaption></figure>

from a snowball state toward a partially ice-free Earth occurs between 2.5 and 2.25 Ga in both types of experiments. However, this apparent synchronicity is due to the relatively long duration increment of 250 Ma used between experiments. Since the incoming solar energy varies by as much as 5 W m⁻² between 2.5 and 2.25 Ga (Fig. 2), it is not possible to evaluate with precision when the pan-glaciation occurs. In addition, it cannot be excluded that in a cold climate environment, the cloud scheme of FOAM underestimates the cloud radiative forcing (Abbot et al., 2012).

## 3.4 Effects of the cloud droplet size

The use of large liquid droplets (20 µm) increases the global mean surface temperature by about 20 to 35 °C compared to present-day droplet sizes (Fig. 3, red curves). In the large-droplet experiments, the shift towards snowball state occurs between 3 and 3.25 Ga with a mean global surface temper-

ature of 13 °C. In contrast, modern clouds experiments result in a non-snowball state at 2.25 Ga with a mean global surface temperature of −19 °C. To explain why large-droplet experiments are in contrast to modern cloud simulations, we examined how clouds affect the Earth’s budget. The cloud radiative forcing (CRF) obtained with large droplets was positive in all experiments (Fig. 4a). As an example, in the 1 Ga run, a net CRF of 14 W m⁻² was obtained using large droplets instead of −4 W m⁻² for modern cloud simulation, a result in agreement with Goldblatt and Zahnle (2011b). In detail, clouds with large droplets generate a positive CRF at high latitudes (+40 W m⁻²), preventing sea ice formation (Fig. 4b). Inhibiting sea ice formation limits the ice-albedo feedback efficiency and maintains the planetary albedo at low values (Fig. 3c, red curve). This change in CRF results from the combination of a slightly lower negative CRF at short wavelength (scattering effect) and a much higher positive CRF at large wavelength (IR trapping; Fig. 4c). As a consequence, even low-elevation clouds (the most reflective clouds) yield a positive CRF, as already shown by Goldblatt and Zahnle (2011b). The increase in evaporation due to warmer oceans enhances the cloud water content. This mechanism increases in turn the greenhouse effect because it amplifies the cloud emissivity ($ε_{cld}$) expressed as

$$ε_{cld} = 1 - e^{(-D * Kabs * CWP)}, \qquad (6)$$

where D is a diffusivity factor (1.66), Kabs is the longwave absorption coefficient (m² g⁻¹) and CWP is the cloud liquid water path.

## 4 Resolution of the faint young Sun problem: sensitivity tests

Fixing the change in the net radiative flux at the tropopause to 19 W m⁻² (Fig. 1), none of the experiments can prevent the onset of pan-glaciation. The Earth’s surface temperature threshold before shifting into a snowball Earth varies from −19 °C, using modern clouds, to +13 °C, using large-droplet clouds. The estimated radiative deficit required for the onset of a pan-glaciation between −57.5 ± 2.5 and −72 ± 2.5 W m⁻² results in a threshold value between 21.5 ± 2.5 and 6.5 ± 2.5 W m⁻², respectively, to solve the FYSP (Fig. 3a). In the following, several warming parameters (greenhouse gas concentrations, ice albedo, a faster rotation rate and/or saltier oceans) able to generate an early Archaean climate at least as warm as today are tested at about 3.5 Ga.

### 4.1 Effects of greenhouse gases and ice albedo

Assuming a solar constant of 1053 W m⁻² (77% of the present-day value), an idealised geography, no ozone layer and two different cloud parameterisations (modern clouds and clouds with large droplets), the radiative forcing was enhanced either by increasing the concentrations of greenhouse
<table><caption>**Table 2.** Resolutions of the Faint Young Sun Problem: sensitivity tests (Sect. 4.1).</caption><thead><tr><th>set of simulations</th><th>Boundary conditions</th><th>pCO<sub>2</sub> (10<sup>-3</sup>bar)</th><th>pCH<sub>4</sub> (10<sup>-6</sup>bar)</th><th>Ice albedo Vis/nIR</th><th>ΔF<sub>trop</sub> (W m<sup>-2</sup>)</th><th>Droplet size (μm)</th><th>Global Ts (°C)</th></tr></thead><tbody><tr><td>Run1</td><td>3.5 Ga</td><td>2.4</td><td>300</td><td>0.7/0.5</td><td>+24</td><td>20</td><td>14.2</td></tr><tr><td>Run2</td><td>“</td><td>25.2</td><td>1.7</td><td>“</td><td>+34</td><td>“</td><td>18.5</td></tr><tr><td>Run3</td><td>“</td><td>16.8</td><td>“</td><td>“</td><td>+30</td><td>“</td><td>15.7</td></tr><tr><td>Run4</td><td>“</td><td>25.2</td><td>“</td><td>“</td><td>+34</td><td>5–10</td><td>-73</td></tr><tr><td>Run5</td><td>“</td><td>56.0</td><td>“</td><td>“</td><td>+41</td><td>“</td><td>15.5</td></tr><tr><td>Run6</td><td>“</td><td>112.0</td><td>“</td><td>“</td><td>+49</td><td>“</td><td>21.1</td></tr><tr><td>Run7</td><td>“</td><td>25.2</td><td>“</td><td>0.5/0.5*</td><td>+34</td><td>“</td><td>-10.5</td></tr></tbody></table>

gases ($pCO_2$ and/or $pCH_4$) or by reducing the ice albedo from 0.6 to 0.5 (e.g. Warren et al., 2002). A summary of the simulations and global mean surface temperatures is given in Table 2.

When clouds with large droplets are considered, a radiative forcing at the tropopause ($\Delta F_{trop}$) of +24 W m$^{-2}$ yields an Earth surface temperature not significantly different from its present-day value (14.2 °C). Such a radiative forcing, which corresponds to 2.4 mbar of $CO_2$ and 0.3 mbar of $CH_4$, is able to solve the FYSP (Fig. 1). These results, which are in agreement with Charney et al. (2013), and support the hypothesis defended by Rosing et al. (2010), are strongly dependent on the boundary conditions used. As an example, considering a modern $pCH_4$, the $pCO_2$ required to maintain a warm climate is about 16.8 mbar (run3, Table 2). When both $pCH_4$ and droplet sizes are fixed to their modern values, a $pCO_2$ of 56 mbar is required to simulate a mean temperature of 15 °C (run5, Table 2). In addition, ice albedo represents another factor that can affect the FYSP. When reducing the ice albedo by 0.1, the model simulates a colder climate, although the equator region remains ice free (run 7, Table 2). Climate response to varying ice albedo suggests the need for a good treatment of sea ice (notably for the sea ice dynamics), and also explains why the calculated greenhouse gas concentrations required for solving the FYSP are model dependent, most of them using different values for ice albedo (0.65 for LMDz; Charney et al., 2013, 0.68 (visible) to 0.3 (nearIR) for CAM3; Wolf and Toon, 2013, and 0.5–0.75 for EMAC, Kunze et al., 2014).

## 4.2 Effects of faster Earth rotation rate and saltier oceans

Palaeozoic rhythmites, mollusc shells and coral fossil records indicate that the Earth rotated faster in the deep past leading to a shorter Length Of the Day (LOD) (Walker and Zahnle, 1986). Between 3.8 and 2.5 Ga, the LOD increased from 14 to 18 h. Jenkins (1993a) and Jenkins et al. (1993) tested the effects of a shorter LOD on the early Earth climate using an atmospheric GCM coupled to a swamp ocean. He demonstrated that a faster rotation (corresponding to a LOD of

~ 14 h) reduces the mean global cloudiness by 20% and increases the mean global surface temperature by 2 °C. In fact this global warming hides local temperature changes, namely a cooling at the poles and a warming at the tropics due to the weakening of transient eddies carrying the heat from low to high latitudes. Hence, faster Earth rotation may have caused the climate to be more sensitive to the ice-albedo feedback (Kienert et al., 2012). To check this assumption, we ran a set of experiments assuming a LOD of 16 h (~ 3.5 Ga). Another peculiar feature of the Archaean Earth is the occurrence of saltier oceans. Knauth (2005) suggested a long-term decline in ocean salinity throughout the Precambrian in parallel with the development of sedimentary basins and climatic conditions able to form and preserve halite (NaCl) deposits. As a consequence, the salinity of Archaean oceans may have been close to 60‰ (in comparison with 35‰ at the present-day) (Knauth, 2005; Foriel et al., 2004). Because a higher salinity lowers the seawater freezing temperature (Millero and Poisson, 1981), sea ice formation could have been delayed, potentially reducing the influence of the ice-albedo feedback. Since the atmospheric GCM FOAM is coupled with a slab ocean, a higher salinity cannot be taken into account. To mimic saltier oceans, we fixed the seawater freezing point at -3.3 °C following the equation of state of seawater proposed by Millero and Poisson (1981). Note that heat diffusion in our mixed-layer ocean was not adjusted in the Earth's rotation rate and salinity experiments (Table 3).

Accounting for saltier oceans, FOAM simulates a warming at high latitudes (Fig. 5), restricting sea ice formation and reducing the Earth's sensitivity to the ice-albedo feedback. More intriguing is the response induced by a faster Earth rotation. In present-day conditions, a change in day length from 24 to 16 h causes a global mean surface warming of about 2 °C (from 14.5 to 16.3 °C, Appendix Fig. A1), in agreement with previous results (Jenkins, 1993b). However, the warming is not homogeneous. Due to the reduced latitudinal heat transport and change in the cloud pattern, a faster Earth rotation warms up low latitudes but cools down high latitudes (Appendix Fig. A1), the warming of low latitudes overcoming the cooling of high latitudes. However,

<footnote>* Visible snow albedo fixed to 0.7 (instead of 0.9).</footnote>
<table><caption>**Table 3.** Simulations carried out with a $\Delta F_{\text{trop}} = +24 \text{ W m}^{-2}$ (Sect. 4.2). Note that heat diffusion in our mixed-layer ocean was not adjusted in the rotation rate and salinity experiments.</caption><thead><tr><th rowspan="2">set of simulations</th><th colspan="2">Boundary conditions</th><th rowspan="2">LOD (hours)</th><th rowspan="2">Droplet size(µm)</th><th rowspan="2">Seawater freezing</th><th rowspan="2">Mean surface temperature</th></tr><tr><th><i>Paleogeography</i></th><th><i>Solar</i></th><th><i>constant</i></th></tr></thead><tbody><tr><td>Standard conditions</td><td>3.5 Ga</td><td>1053 W m<sup>-2</sup></td><td>24</td><td>20</td><td>-1.8 °C</td><td>14.2 °C</td></tr><tr><td>Saltier ocean</td><td>“</td><td>“</td><td>24</td><td>“</td><td>-3.3 °C</td><td>15.0 °C</td></tr><tr><td>Faster Rotation</td><td>“</td><td>“</td><td>16</td><td>“</td><td>-1.8 °C</td><td>14.2 °C</td></tr><tr><td>Early Earth</td><td>“</td><td>“</td><td>16</td><td>“</td><td>-3.3 °C</td><td>14.4 °C</td></tr></tbody></table>

<figure><img src="image_63.png" /><figcaption>**Fig. 5.** Sea surface temperature anomalies ($\Delta$SST) as a result of a faster Earth rotation rate (LOD = 16 h) and/or saltier ocean. Simulations performed with 3.5 Ga boundary conditions with large-droplet clouds and a $\Delta F_{\text{trop}} = +24 \text{ W m}^{-2}$ (standard conditions). The saltier ocean experiment uses standard conditions associated with an ocean two times saltier; the faster rotation experiment has a length of day of 16 h; the early Earth experiment combines both.</figcaption></figure>

in colder conditions, a shorter LOD leads to the opposite trend: the sea ice formation at high and mid latitudes is facilitated, increasing the albedo and amplifying the cooling (Fig. 5). Consequently, when both a shorter LOD and saltier oceans are considered, their effects are balanced (a shorter LOD increases the Earth's sensitivity to the ice-albedo feedback while saltier oceans reduce it). This explains why these processes do not affect the solution of the FYSP (Fig. 5). It is noteworthy that our faster rotation simulation agrees with LMDz, but not CLIMBER (Kienert et al., 2012). One of the main reasons for such a difference is discussed below.

## 5 Discussion

FOAM, like other atmospheric GCMs, relies on a physical basis, which means that mechanisms and feedbacks are explicitly accounted for. Goldblatt and Zahnle (2011b) challenged the use of GCMs given that sub-grid scale parameterisation of clouds in these models precludes their use in palaeoclimate modelling. Here, we discuss differences between AGCMs and other categories of climate models (RCMs and EMIC) and argue that AGCMs are relevant and accurate tools for investigating deep-time climates.

### 5.1 FYSP resolutions: FOAM versus CLIMBER (EMIC)

FOAM results largely disagree with CLIMBER (Kienert et al., 2012, 2013). For instance, using CLIMBER, Kienert et al. (2012, 2013) suggested that a $p\text{CO}_2$ of 0.6 bar is required to generate a mean global temperature of 15 °C (for a 3.8 Ga solar constant of $1023 \text{ W m}^{-2}$) and that a faster rotation rate strongly impacts the ice-albedo feedback. This is at odd with the results obtained with FOAM in the present study. One of the reasons why CLIMBER yields opposite results to atmospheric GCMs essentially because of the different lapse rate parameterisation. Kienert et al. (2013) used Outgoing Longwave Radiation in clear-sky diagnostic (OLR) to check the validity of longwave radiative transfer, fixing the thermal structure of the atmosphere according to the lapse rate given in Halevy et al. (2009). However, they omitted to quantify OLR fluctuations caused by CLIMBER lapse rate parameterisation. As shown by Pierrehumbert (2004), the lapse rate
<figure><img src="image_64.png" /><figcaption>Fig. 6. Lapse rates (annual mean) calculated by CLIMBER lapse rate formula and FOAM using 1 Ga conditions (no clouds). Black curve represents the standard formulation in CLIMBER (Petoukhov et al., 2000), the red line is the calculation adapted for a faster rotation in CLIMBER (Kienert et al., 2012, 2013). The green line represents the FOAM lapse rate.</figcaption></figure>

governs the sensitivity of climate to CO₂ so that the greenhouse effect (G) is expressed as

$$G(W m^{-2}) = \sigma T s^4 - OLR, \quad (7)$$

where Ts is the surface temperature (K), σ the Stefan-Boltzmann constant. OLR corresponds to the energy lost by the Earth as infrared radiation at the top of the atmosphere.

The CLIMBER lapse rate is based on a theoretical approach (Petoukhov et al., 2000) and was developed to fit empirical data such as

$$\text{lapse rate} = \Gamma_0 + \Gamma_1 (Ts - 273.16 \text{ K}) * (1 - a_q * Q^2) - \Gamma_2 * n_{\text{cld}}, \quad (8)$$

where Γ₀, Γ₁ and Γ₂ are lapse rate parameters (Γ₀ = 5.2 × 10⁻³ K m⁻¹, Γ₁ = 5.5 × 10⁻⁵ K m⁻¹, Γ₂ = 10⁻³ K m⁻¹), Q is the surface-specific humidity (kg kg⁻¹), ncld the cumulus cloud amount and aq = 10³ (kg kg⁻¹). Kienert et al. (2012, 2013) updated this formula to account explicitly for a faster Earth rotation rate such as Γ₀ = 5.1 × 10⁻³ K m⁻¹, Γ₁ = 6.1 × 10⁻⁵ K m⁻¹, Γ₂ = 0.9⁻³ K m⁻¹ and aq = 1.171 × 10³ (kg kg⁻¹). To illustrate CLIMBER lapse rate parameterisation better, we used FOAM outputs (i.e. Ts and Q). To suppress the cloud effects, we extracted outputs from a run performed at 1 Ga with no clouds; therefore ncld is fixed to 0 (lapse rates computed with clouds are available in Appendix Fig. A3). In CLIMBER the radiative transfer depends on Ht (troposphere height), so we calculated lapse rates for the entire troposphere (from the Earth's surface to an atmospheric pressure level of 254 mbar).

Figure 6 shows lapse rates for different case scenarios: (1) with the standard formula implemented in CLIMBER

(Petoukhov et al., 2000), (2) using the adapted formula in case of a faster rotation (Kienert et al., 2012, 2013) and (3) using FOAM. For the same values of Ts and Q, we show that the lapse rate is systematically underestimated with CLIMBER, hence inducing warmer temperature and a higher thermal radiation flux loss at the tropopause (the unusual shape of the lapse rate is explained by the northern hemisphere which is almost entirely oceanic, see Fig. 2b). Therefore, for the same pCO₂, CLIMBER results in a less efficient greenhouse effect. This artefact in the CLIMBER lapse rate explains why this EMIC needs such a large pCO₂ when a present-day Earth rotation is considered (Kienert at al., 2012). The required pCO₂ is even higher when a faster Earth rotation is imposed. Although the authors discussed the uncertainties linked to CLIMBER radiative module (Kienert et al., 2013), they omitted lapse rate/OLR interplays. In addition to possible effects of ocean and sea ice dynamics and an ice albedo above classical values, this omission could explain why Kienert et al. (2012, 2013) proposed such a high pCO₂ to solve the FYSP.

## 5.2 Pan-glaciation thresholds and ice-albedo feedback: AGCMs versus RCMs

It is noteworthy that atmospheric GCMs are able to maintain a large fraction of open oceans (i.e. above the seawater freezing point) with very moderate amounts of greenhouse gases, 0.9 mbar of CO₂ and CH₄ for LMDz (3.8 Ga), 3.7 mbar of CO₂ for EMAC (3.5 Ga) and 5 mbar of CO₂ for CAM3 (2.7 Ga), a result that contrasts with RCM estimates (0.04 bar of CO₂ at 3.5 Ga, Kasting, 1987). RCMs assume that a drop of the mean temperature below 0 °C leads to a runaway glaciation, and therefore induces high pCO₂ in order to maintain ice-free oceans during the Archaean. However, stable climates have been simulated using AGCMs while the global mean surface temperature falls well below 0 °C (Charney et al., 2013; Wolf and Toon, 2013; Kunze et al., 2014, this study). Wolf and Toon (2013) suggested that the difference in the temperature leading to a runaway glaciation was the reason for the high estimates of pCO₂ by RCMs.

To better understand why RCMs are not able to calculate the threshold temperature allowing the onset of runaway glaciations, we implemented an albedo dependency on temperature in the RCM developed by MacKay and Khalil (1991). It should be noted that Wang and Stone (1980) incorporated this process, but for unknown reasons, all RCM experiments kept a surface albedo held constant. This lack had biased the climate response and the Earth was maintained unfrozen for unrealistic reasons. To represent the ice-albedo feedback, we used the formula given by Pierrehumbert et al. (2011) where the Earth albedo is assumed to be

- Ice-free for $T > T_0 = 20$ °C, with an albedo ($\alpha_0$) = 0.2

- Totally ice covered for $T < T_{\text{freeze}}$, with an albedo ($\alpha_{\text{freeze}}$) = 0.6
<figure><img src="image_65.png" /><figcaption>Fig. 7. RCM simulations showing the evolution of the surface temperature versus time between 3.5 and 1 Ga, with and without the ice-albedo feedback, and assuming increasing irradiance in the absence of clouds and a greenhouse radiative forcing of 19 W m⁻². Blue, green and red curves were obtained using the ice-albedo feedback formula from Pierrehumbert et al. (2011) where T<sub>freeze</sub> values (−20 °C, −5 °C, 0 °C) represent the freezing temperature for the Earth (snowball Earth threshold). The black curve is obtained using a surface albedo fixed at 0.264 (absence of ice-albedo feedback). Rosing et al. (2010) and FOAM runs are described in Fig. 3.</figcaption></figure>

- Partially ice covered for T<sub>freeze</sub> < T < T<sub>0</sub>, the albedo is interpolated according to the formula:

$$ \alpha(T) = \alpha_o + (\alpha_{\text{freeze}} - \alpha_o)[(T - T_0)^2/(T_{\text{freeze}} - T_0)^2] \quad (9) $$

where T is the Earth's temperature and T<sub>freeze</sub> is the threshold temperature for the onset on global glaciations. The difference between T<sub>freeze</sub> and T<sub>0</sub> is a measure of the meridional temperature gradient and therefore represents the efficiency of latitudinal heat transport. Because Wolf and Toon (2013) proposed that the freezing temperature is the main factor explaining differences between RCMs and AGCMs, we tested several T<sub>freeze</sub> values (−20 °C, −5 °C, 0 °C). To avoid the difficulty of figuring out how to represent cloud changes between FOAM and RCM we removed clouds, so that the albedo corresponds to the planetary albedo incorporating the effects of Rayleigh scattering, atmospheric solar absorption and the reflective properties of the surface.

The RCM and FOAM are forced with the same boundary conditions, that is a net radiative forcing of +19 W m⁻² and no clouds. In the absence of the ice-albedo feedback (α = 0.264), the RCM simulates a linear increase in the mean temperature through time, the threshold value of 0 °C (supposed to represent the threshold temperature for runaway glaciations) being reached at about 2.6 Ga (Fig. 7). In contrast, RCM experiments that include the ice-albedo feedback show a drastic shift in the evolution of the mean temperature

through time, similar to the one observed in FOAM experiments. More interestingly, when the ice-albedo feedback is parameterised such as the threshold temperature for the onset of snowball Earth (T<sub>freeze</sub>) equals −20 °C, the timing of the shift from snowball Earth to ice-free conditions is similar in RCM and FOAM experiments. This similarity occurs because RCM mimics the latitudinal air surface temperature gradient of FOAM (ΔT = 40 °C from equator to 45° of latitudes – sea ice latitudes). This calculation illustrates an essential weakness of RCM: their inability to calculate T<sub>freeze</sub> in the absence of latitudinal heat transport. In AGCMs, T<sub>freeze</sub> varies with the prescribed boundary conditions (for example in FOAM, the maximum latitudinal extent of the ice cap before the shift to runaway glaciation varies from 30° to 45° and 70° for experiments using modern clouds, cloud-free, and large-droplet clouds, respectively). This extent is controlled by the change in ocean heat flux, a process which cannot be captured by RCMs. Faster rotation and saltier oceans also slightly affect T<sub>freeze</sub>. Consequently, AGCMs, in contrast to RCMs, more faithfully represent the ice-albedo feedback and hence the climatic transition towards snowball Earth.

## 5.3 Pan-glaciation thresholds and ice-albedo feedback: FOAM versus other AGCMs

The recent use of atmospheric GCMs for investigating the FYSP offered a great opportunity to compare FOAM with CAM3 (Wolf and Toon, 2013) LMDz (Charney et al., 2013) and EMAC (Kunze et al., 2014). All these studies suggested that Earth is less prone to full glaciations. Using a 2.7 Ga-old solar constant of 1093 W m⁻² and 1 mbar of CO₂ and CH₄ (other parameters being held at present-day values), CAM3 simulates a mean surface temperature close to 2 °C (Wolf and Toon, 2013). For the same conditions, FOAM yields a full glaciation. A small difference in ice albedo of 0.68(vis)/0.3 (nearIR) with CAM3 and 0.7(vis)/0.5(nearIR) with FOAM seems to explain the different behaviour. Although minor, this difference has important consequences on the impact of the ice-albedo feedback when Earth becomes ice-covered (ice line spreads down to 30–40°; run 7 of Table 2). Using a 3 Ga-old solar constant (1082 W m⁻²) and 0.9 mbar of CO₂, LMDz does not simulate a full glaciation (Charney et al., 2013). This behaviour arises from both the prescribed boundary conditions and heat transport. Charney et al. (2013) made the hypothesis that Earth was an aquaplanet at 3 Ga. In such conditions, LMDz simulates larger amount of water vapour, which enhances heat transport to high latitudes and decreases the equator-pole temperature gradient, so that the ice-albedo feedback becomes less efficient and the runaway glaciation is delayed. This result is in good agreement with the result of Kunze et al. (2014) who also assumed a water-covered world for the early Archaean (3.5 Ga, solar constant 1046 W m⁻²). Using EMAC, they showed that it was possible to get a stable ice-free ocean at low latitudes (surface temperature equals to −9 °C) with moderate CO₂ concentrations (3670 ppmv).
For comparison, LMDz prevents the early Earth (3.8 Ga, 1025 W m⁻²) from completely freezing with a greenhouse forcing corresponding to 0.9 mbar of a CO₂ combined to 0.9 mbar of CH₄ (Charney et al., 2013).

LMDz simulations also investigated the effect of CCN. Performing simulations with larger droplets (droplet radius set to 17 instead of 20 µm in our study), Charney et al. (2013) suggested a strong impact on the terrestrial radiative budget, in agreement with this study. More interestingly, both AGCMs simulate an increase in the cloud water path in low-altitude clouds located at high latitudes. Since the use of these two different cloud schemes shows similar behaviour, we suspect that the warming driven by high-latitude clouds inhibiting sea ice formation is not model dependent (see Sect. 3.4).

Other processes such as pressure broadening of CO₂ and H₂O by high atmospheric pressure (Goldblatt et al., 2009) or high concentrations of H₂ (Wordsworth and Pierrehumbert, 2013), which are not treated here, could also warm up the early Earth and affect the FYSP solutions. In addition, several processes are missing in this study, which are potential sources of uncertainties. The use of a slab ocean, with heat transport parameterised through diffusion, precludes the faithful representation of the meridional ocean heat transport. Indeed, during the Archaean the oceanic heat transport should supposedly be enhanced due to the smaller land fraction (Endal and Schatten, 1982). We also acknowledge that saltier oceans possibly lead to a weaker overturning (Williams et al., 2010). These conflicting assessments concerning ocean circulation, and their respective influence on the FYSP, could be explored using a fully coupled dynamic ocean model, an investigation which has not been conducted here. The last limitation relative to the use of a slab ocean is the absence of sea ice dynamics, which tends to slow the sea ice flow (Voigt and Abbot, 2012). This lack may change the radiative deficit triggering the onset of pan-glaciation and underestimate the greenhouse gas concentrations required to keep the early Earth unfrozen.

Obviously atmospheric GCMs still have limitations (Hack, 1998; Abbot et al., 2012; Stevens and Bony, 2013) and results can slightly differ according to their respective parameterisations or boundary conditions. However, as fully dynamic models, AGCMs represent a powerful tool for improving our understanding of deep past Earth climates, capturing cloud mechanisms involved in global warming, as well as time-dependent factors (e.g. palaeogeography, rotation rate). However, constraints on paleo-pCO₂ gained from climate modelling studies cannot be directly compared to geological data because all these studies omit to consider the carbon cycle-climate equilibrium. To estimate more accurately pCO₂ and the Earth's surface temperature evolution, we associated the atmospheric GCM simulations discussed in this paper with a geochemical model to investigate the carbon cycle climate during the Archaean (see companion paper).

## 6 Conclusions

FOAM, like other AGCMs, suggests a faint young Sun problem less problematic than initially assumed and solves the FYSP with a composition of greenhouse gases one order a magnitude lower than RCM or EMIC. Accounting for a 23 % weaker Sun at 3.5 Ga, we show that a moderate radiative forcing of 41 W m⁻² is sufficient to maintain the mean global surface temperature close to its present-day value. Assuming a high pCH₄ in the Archaean atmosphere (Pavlov et al., 2000; Catling et al., 2007), we show that such a moderate radiative forcing can be provided by the combination of 1 (0.3) mbar of methane and 11.5 (16) mbar of carbon dioxide. The presence of large-droplet clouds can theoretically provide significant warming so that lower CO₂ and CH₄ partial pressures are required for the maintenance of an ice-free surface. We demonstrated, however, that such changes in cloud properties cannot alone provide sufficient warming to solve the FYSP with low (0.9 mbar) pCO₂ as defended by Rosing et al. (2010).

This study highlights the central role of AGCMs over RCMs in investigating glacial thresholds. RCMs overestimate the early Archaean pCO₂ (i.e. 0.2–0.3 bar) required to keep equatorial waterbelts ice free due to their inability to calculate the correct threshold mean temperature for the onset of full glaciation (i.e. T freeze). Applying Archaean boundary conditions, our AGCM simulations show a wide range of possible values for T freeze, from –20 to 13 °C, according to the latitudinal temperature gradient. In contrast, RCMs use a fixed T freeze value (0 °C, the freezing point of seawater) and therefore inaccurately estimate the poleward heat transport. The weak sensitivity of CLIMBER to CO₂ could be explained by its lapse rate calculation, which overestimates the energy leaving the Earth as infrared radiation. Our AGCM supports the main conclusions arising from previous AGCMs, but predicts an early Earth more inclined to full glaciation. Higher ice albedo is a likely cause for the stronger ice-albedo feedback in FOAM compared to CAM3. LMDz and EMAC simulate warmer climate due to a weaker ice-albedo feedback inherited from the absence of continents, a boundary condition that enhances poleward moist heat fluxes. However, all these studies cannot simulate ocean and sea ice dynamics, a limitation which may underestimate the greenhouse gas concentrations required to keep the early Earth unfrozen. Following this step forward in the understanding of physical processes able to solve the FYSP, we now consider carbon cycle–climate interplays to provide a complete picture of the faint young Sun problem (companion paper).

*Acknowledgements.* We thank the two anonymous referees for their very insightful reviews which considerably improved the quality of the manuscript, and the editorial board for the quality of the review procedure. This work was granted access to the HPC resources of CCRT under allocation 2014-017013 made
by GENCI (Grand Equipement National de Calcul Intensif). This work was supported by ANR eLIFE2 and LabexUnivEarth (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). This is IPGP contribution 3494.

Edited by: M. Claussen

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

The publication of this article is financed by CNRS-INSU.

## References

Abbot, D. S., Voigt, A., Branson, M., Pierrehumbert, R. T., Pollard, D., Le Hir, G., and Koll, D. D. B.: Clouds and snowball Earth deglaciation, Geophys. Res. Lett., 39, L20711, doi:10.1029/2012GL052861, 2012.

Andreae, M. O. and Rosenfeld, D.: Aerosol–cloud–precipitation interactions. Part 1. The nature and sources of cloud-active aerosols, Earth-Sci. Rev., 89, 13–41, doi:10.1016/j.earscirev.2008.03.001, 2008.

Boucher, O., Le Treut, H., and Baker, M. B.: Precipitation and radiation modeling in a general circulation model: Introduction of cloud microphysical processes, J. Geophys. Res., 100, 16395–16414, doi:10.1029/95JD01382, 1995.

Bréon, F. M., Tanré, D., and Generoso, S.: Aerosol Effect on Cloud Droplet Size Monitored from Satellite, Science, 295, 834–838, doi:10.1126/science.1066434, 2002.

Catling, D. C., Claire, M. W., and Zahnle, K. J.: Anaerobic methanotrophy and the rise of atmospheric oxygen, Phil. Trans. R. Soc., 365, 1867–1888, doi:10.1098/rsta.2007.2047, 2007.

Charlson, R. J., Lovelock, J. E., Andreae, M. O., and Warren, S. G.: Oceanic Phytoplankton, atmospheric sulfur, cloud albedo and climate, Nature, 326, 655–661, 1987.

Charney, B., Forget, F., Wordworth, R., Leconte, J., Millour, E., Codron, F., and Spiga, A.: Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3-D GCM, J. Geophys. Res., 118, 10414–10431, doi:10.1002/jgrd.50808, 2013.

Dauphas, N. and Kasting, J. F.: Low pCO₂ in the pore water, not in the Archean air, Nature, 474, E2–E3, doi:10.1038/nature09960, 2011.

Driese, S. G., Jirsa, M. A., Ren, M., Brantley, S. L., Sheldon, N. D., Parker, D., and Schmitz, M.: Neoarchean paleoweathering of tonalite and metabasalt: implications for reconstructions of 2.69 Ga early terrestrial ecosystems and paleoatmospheric chemistry, Precambrian Res., 187, 1–17, doi:10.1016/j.precamres.2011.04.003, 2011.

Farquhar, J., Airieau, S., and Thiemens, M. H.: Observation of wavelength-sensitive mass-independent sulfur isotope effects during SO₂ photolysis: Implications for the early atmosphere, J. Geophy. Res.-Planets, 106, 32829–32839, 2001.

Feulner, G.: The faint young sun problem, Rev. Geophys., 50, RG2006, doi:10.1029/2011RG000375, 2012.

Foriel, J., Philippot, P., Rey, P., Somogyi, A., Banks, D., and Menez, B.: Biological control of Cl/Br and low sulfate concentration in a 3.5-Gyr-old seawater from North Pole, Western Australia, Earth Planet. Sci. Lett., 228, 451–463, 2004.

Godderis, Y. and Veizer, J.: Tectonic control of chemical and isotopic composition of ancient oceans: the impact of continental growth, Am. J. Sci., 300, 434–461, 2000.

Goldblatt, C. and Zahnle, K. J.: Faint Young Sun Paradox remains, Nature, 474, E1, doi:10.1038/nature09961, 2011a.

Goldblatt, C. and Zahnle, K. J.: Clouds and the Faint Young Sun Paradox, Clim. Past, 7, 203–220, doi:10.5194/cp-7-203-2011, 2011b.

Goldblatt, C., Claire, M. W., Lenton, T. M., Matthews, A. J., Watson, A. J., and Zahnle, K. J.: Nitrogen-enhanced greenhouse warming on early Earth, Nat. Geosci., 2, 891–896, doi:10.1038/NGEO692, 2009.

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

Hack, J. J.: Sensitivity of the simulated climate to a diagnostic formulation for cloud liquid water, J. Climate, 11, 1497–1515, 1998.

Halevy, I., Pierrehumbert, R. T., and Schrag, D. P.: Radiative transfer in CO₂-rich paleoatmospheres, J. Geophys. Res., 114, D18112, doi:10.1029/2009JD011915, 2009.

Hansen, J., Sato, M., and Ruedy, R.: Radiative Forcing and climate response, J. Geophys. Res., 102, 6831–6864, 1997.

Haqq-Misra, J. D., Domagal-Goldman, S. D., Kasting, P. J., and Kasting, J. F.: A revised, hazy methane greenhouse for the archean Earth, Astrobiology, 8, 1127–1137, 2008.

Jacob, R.: Low frequency variability in a simulated atmosphere ocean system, PhD Thesis, University of Wisconsin-Madison, 1997–310, 1997.

Jenkins, G. S.: The effects of reduced land fraction and solar forcing on the general circulation results from the NCAR CCM, Glob. Planet. Change, 7, 321–333, 1993a.

Jenkins, G. S.: A General Circulation Model Study of the Effects of Faster Rotation Rate, Enhanced CO₂ Concentration, and Reduced Solar Forcing Implications for the Faint Young Sun Paradox, J. Geophys. Res., 98, 20803–20811, 1993b.

Jenkins, G. S., Marshall, H. G., and Kuhn, W. R.: Precambrian climate – the effects of land area and Earth’s rotation rate, J. Geophys. Res.-Atmos., 98, 8785–8791, doi:10.1029/93JD00033, 1993.

Kasting, J. F.: The evolution of the prebiotic atmosphere, Origins Life Evol. B., 14, 75–82, 1984.

Kasting, J. F.: Theoretical constraints on oxygen and carbon-dioxide concentrations in the Pre-cambrian atmosphere, Precambrian Res., 34, 205–229, 1987.

Kasting, J. F.: Early Earth faint young sun redux, Nature, 464, 687–689, 2010.

Kiehl, J. T. and Dickinson, R. E.: A study of the radiative effects of enhanced atmospheric CO₂ and CH₄ on early Earth surface temperatures, J. Geophys. Res.-Atmos., 92, 2991–2998, 1987.

Kienert, H., Feulner, G., and Petoukhov, V.: Faint young sun problem more severe due to ice-albedo feedback and higher rotation rate of the early Earth, Geophys. Res. Lett., 39, L23710, doi:10.1029/2012GL054381, 2012.
Kienert, H., Feulner, G., and Petoukhov, V.: Albedo and heat transport in 3-D model simulations of the early Archean climate, Clim. Past, 9, 1841–1862, doi:10.5194/cp-9-1841-2013, 2013.

Knauth, L. P.: Temperature and salinity history of the Precambrian ocean: implications for the course of microbial evolution, Palaeogeogr. Palaeocl., 219, 53–69, 2005.

Kump, L. R. and Pollard, D.: Amplification of cretaceous warmth by biological cloud feedbacks, Science, 320, 195–195, 2008.

Kunze, M., Godolt, M., Langematz, U., Grenfell, J. L., Hamann-Reinus, A., and Rauer, H.: Investigating the early Earth faint young Sun problem with a general circulation model, Planet. Space Sci., doi:10.1016/j.pss.2013.09.011, in press, 2014.

Le Hir, G., Donnadieu, Y., Krinner, G., and Ramstein, G.: Toward the snowball earth deglaciation, Clim. Dynam., 35, 285–297, 2010.

Lowe, D. R. and Tice, M. M.: Tectonic controls on atmospheric, climatic, and biological evolution 3.5–2.5 Ga, Precambrian Res., 158, 177–197, 2007.

MacKay, R. M. and Khalil, M. A. K.: Theory and Development of a One Dimensional Time Dependant Radiative Convective Climate Model, Chemosphere, 22, 383–417, 1991.

Millero, F. J. and Poisson, A.: International one-atmosphere equation of state of seawater, Deep-Sea Res. Pt. I, 28, 625–629, 1981.

Pavlov, A. A., Kasting, J. F., Brown, L. L., Rages, K. A., and Freedman, R.: Greenhouse warming by CH₄ in the atmosphere of early Earth, J. Geophys. Res.-Planets, 105, 11981–11990, 2000.

Penner, J. E., Quaas, J., Storelvmo, T., Takemura, T., Boucher, O., Guo, H., Kirkevåg, A., Kristjánsson, J. E., and Seland, Ø.: Model intercomparison of indirect aerosol effects, Atmos. Chem. Phys., 6, 3391–3405, doi:10.5194/acp-6-3391-2006, 2006.

Petoukhov, V., Ganopolski, A., Brovkin, V., Claussen, M., Eliseev, A., Kubatzki, C., and Rahmstorf, S.: CLIMBER-2: a climate system model of intermediate complexity. Part I: model description and performance for present climate, Clim. Dynam. 16, 1–17, doi:10.1007/PL00007919, 2000.

Pesonen, L. J., Elming, S. A, Mertanen, S., Pisarevsky, S., D'Agrella-Filho, M. S., Meert, J. G., Schmidt, P. W., Abrahamsen, N., and Bylund, G.: Palaeomagnetic configuration of contaminants during the Proterozoic, Tectonophysics, 375, 289–324, 2003.

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

Pierrehumbert, R. T.: Principles of Planetary Climate, Cambridge University Press, 652 pp., 2010.

Pierrehumbert, R. T., Abbot, D. S., Voigt, A., and Koll, D.: Climate of the Neoproterozoic Annual, Rev. Earth Planet. Sci., 39, 417–460, doi:10.1146/annurev-earth-040809-152447, 2011.

Reinhard, C. T. and Planavsky, N. J.: Mineralogical constraints on Precambrian pCO₂, Nature, 474, E1–E2, doi:10.1038/nature09959, 2011.

Rosing, M. T., Bird, D. K., Sleep, N. H., and Bjerrum, C. J.: No climate paradox under the faint early Sun, Nature, 464, 744–747, doi:10.1038/nature08955, 2010.

Rye, R., Kuo, P. H., and Holland, H. D.: Atmospheric carbon dioxide concentrations before 2.2-billion years ago, Nature, 378, 603–605, 1995.

Sagan, C. and Mullen, G.: Earth and Mars – evolution of atmospheres and surface temperatures, Science, 177, 52–56, doi:10.1126/science.177.4043.52, 1972.

Schatten, K. H. and Endal, A. S.: The faint young Sun climate paradox – volcanic influence, Geophys. Res. Lett., 9, 1309–1311, doi:10.1029/GL009i012p01309, 1982.

Sheldon, N. D.: Precambrian paleosols and atmospheric CO₂ levels, Precambrian Res., 147, 148–155, 2006.

Slingo, A.: A GCM parameterization for the shortwave radiative properties of water clouds, J. Atmos. Sci., 46, 1419–1427, 1989.

Stevens, B. and Bony, S.: What Are Climate Models Missing?, Science, 340, 1053–1054, doi:10.1126/science.1237554, 2013.

Voigt, A. and Abbot, D. S.: Sea-ice dynamics strongly promote Snowball Earth initiation and destabilize tropical sea-ice margins, Clim. Past, 8, 2079–2092, doi:10.5194/cp-8-2079-2012, 2012.

von Paris, P., Rauer, H., Grenfell, J. L., Patzer, B., Hedelt, P., Stracke, B., Trautmann, T., and Schreier, F.: Warming the early earth-CO₂ reconsidered, Planet. Space Sci., 56, 1244–1259 doi:10.1016/j.pss.2008.04.008, 2008.

Walker, J. C. G.: Climatic factors on the Archean Earth, Palaeogr. Palaeocl., 40, 1–11, doi:10.1016/0031-0182(82)90082-7, 1982.

Walker, J. C. G. and Zahnle, K. J. Lunar nodal tide and distance to the Moon during the Pre-cambrian, Nature, 320, 600–602, 1986.

Walker, J. C. G., Hays, P. B., and Kasting, J. F.: A Negative Feedback Mechanism For The Long-Term Stabilization of Earth's Surface-Temperature, J. Geophys. Res.-Ocean. Atmos., 86, 9776–9782, 1981.

Wang, W. and Stone, P. H.: Effect of Ice-Albedo Feedback on Global Sensitivity in a One-Dimensional Radiative-Convective Climate Model, J. Atmosph. Sci., 37, 545–552, 1980.

Warren, S. G., Brandt, R. E., Grenfell, T. C., and McKay, C. P.: Snowball Earth: Ice thickness on the tropical ocean, J. Geophys. Res.-Ocean., 107, 31-1-31-18, doi:10.1029/2001JC001123, 2002.

Williams, P. D., Guilyardi, E., Madec, G., Gualdi, S., and Scoccimarro, E.: The role of mean ocean salinity in climate, Dynam. Atmos. Oceans, 49, 108–123, doi:10.1016/j.dynatmoce.2009.02.001, 2010.

Wolf, E. T. and Toon, O. B.: Hospitable Archean Climates Simulated by a General Circulation Model, Astrobiology, 13, 656–673, doi:10.1089/ast.2012.0936, 2013.

Wordsworth, R. and Pierrehumbert, R.: Hydrogen-nitrogen greenhouse warming in Earth's early atmosphere, Science, 339, 64–67, doi:10.1126/science.1225759, 2013.
Appendix A

<figure><img src="image_67.png" /><figcaption>**Fig. A1.** (a) Temperature gradient as a function of Earth rotation rates. (b) Atmospheric lapse rate as a function of ozone layer presence/removal. Simulations performed under present-day boundary conditions.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>**Fig. A2.** Specific humidity response to varying air temperature. The Earth's global temperature equals 12.9 °C (1 Ga, modern clouds), 30.4 °C (1 Ga, large-droplet clouds) and 12.3 °C (1 Ga, large-droplet clouds).</figcaption></figure>
<figure><img src="image_69.png" /><figcaption>**Fig. A3.** Lapse rates (annual mean) using 1 Ga boundary conditions (modern clouds parameterization), calculated by the CLIMBER lapse rate formula and FOAM using the same conditions. Black curve represents the standard formulation in CLIMBER (Petoukhov et al., 2000), the red line is the calculation adapted for a faster rotation in CLIMBER (Kienert et al., 2012, 2013). The green line represents the FOAM lapse rate.</figcaption></figure>
<figure><img src="image_70.png" /></figure>

# Effect of the Ordovician paleogeography on the (in)stability of the climate

A. Pohl¹, Y. Donnadieu¹, G. Le Hir², J.-F. Buoncristiani³, and E. Vennin³

¹LSCE – Laboratoire des Sciences du Climat et de l'Environnement, UMR8212 – CNRS-CEA-UVSQ, CEA Saclay, Orme des Merisiers, 91191 Gif-sur-Yvette Cedex, France

²IPGP – Institut de Physique du Globe de Paris, Université Paris7-Denis Diderot, 1 rue Jussieu, 75005 Paris, France

³Laboratoire Biogéosciences, UMR/CNRS 6282, Université de Bourgogne, 6 Bd Gabriel, 21000 Dijon, France

*Correspondence to:* A. Pohl (alexandre.pohl@lsce.ipsl.fr)

Received: 23 May 2014 – Published in Clim. Past Discuss.: 2 July 2014

Revised: 8 October 2014 – Accepted: 14 October 2014 – Published: 26 November 2014

**Abstract.** The Ordovician Period (485–443 Ma) is characterized by abundant evidence for continental-sized ice sheets. Modeling studies published so far require a sharp CO₂ drawdown to initiate this glaciation. They mostly used non-dynamic slab mixed-layer ocean models. Here, we use a general circulation model with coupled components for ocean, atmosphere, and sea ice to examine the response of Ordovician climate to changes in CO₂ and paleogeography. We conduct experiments for a wide range of CO₂ (from 16 to 2 times the preindustrial atmospheric CO₂ level (PAL)) and for two continental configurations (at 470 and at 450 Ma) mimicking the Middle and the Late Ordovician conditions. We find that the temperature-CO₂ relationship is highly non-linear when ocean dynamics are taken into account. Two climatic modes are simulated as radiative forcing decreases. For high CO₂ concentrations (≥ 12 PAL at 470 Ma and ≥ 8 PAL at 450 Ma), a relative hot climate with no sea ice characterizes the warm mode. When CO₂ is decreased to 8 PAL and 6 PAL at 470 and 450 Ma, a tipping point is crossed and climate abruptly enters a runaway icehouse leading to a cold mode marked by the extension of the sea ice cover down to the midlatitudes. At 450 Ma, the transition from the warm to the cold mode is reached for a decrease in atmospheric CO₂ from 8 to 6 PAL and induces a ~ 9°C global cooling. We show that the tipping point is due to the existence of a 95% oceanic Northern Hemisphere, which in turn induces a minimum in oceanic heat transport located around 40°N. The latter allows sea ice to stabilize at these latitudes, explaining the potential existence of the warm and of the cold climatic modes. This major climatic instability potentially brings a new ex-

planation to the sudden Late Ordovician Hirnantian glacial pulse that does not require any large CO₂ drawdown.

## 1 Introduction

The Ordovician Period (485–443 Ma) is characterized by fundamental changes in the living organisms of our planet. Following the Cambrian explosion, the Early Ordovician is characterized by a major radiation in marine life and critical changes in paleoecology during the Great Ordovician Bio-diversification Event (GOBE) through the rise of the Paleozoic evolutionary fauna (Servais et al., 2010). This time of marine diversification is yet interrupted by the second of the five biggest extinctions of the Phanerozoic Eon in terms of the percentage of genera and families lost (Sheehan, 2001). Trotter et al. (2008) proposed a critical role of the oceanic thermal state during the GOBE, showing that the long-term cooling from a very warm ocean to present-day-like temperatures during the Early to Middle Ordovician may have led to optimal conditions for widespread taxonomic radiations and the appearance of complex communities. Cooler conditions in the Late Ordovician are generally associated with the onset of large ice sheets over the southern Gondwana (e.g., Ghienne et al., 2007). The growth and the reduction of these ice sheets are concomitant with the two pulses of the global extinction (Sheehan, 2001). Thus, a very close binding relationship may have existed between global climate and organismal evolution during the Ordovician.
Atmospheric CO₂ concentration largely constrained climate evolution over the history of the Earth (Berner, 2006). The Ordovician glaciation stands out from the crowd since it occurred under high CO₂ values. Climate proxies (Yapp and Poths, 1992) and continental weathering models (Berner, 2006; Nardin et al., 2011) indicate that the Ordovician CO₂ atmospheric partial pressure (pCO₂) was equal to some 8–20 times its preindustrial atmospheric level (1 PAL = 280 ppm). Faintner Sun (Gough, 1981) somewhat reduces this apparent discrepancy between high CO₂ levels and glacial inception by suggesting that the solar constant was lower in the past. Nevertheless, it seems insufficient; the best hypotheses today consider a Late Ordovician drop in CO₂ to be the trigger of glaciation (Brenchley et al., 1994; Kump et al., 1999; Villas et al., 2002; Lenton et al., 2012). A glacial threshold of 14 to 8 PAL has been found with various models, from a conceptual energy balance model (EBM) (Crowley and Baum, 1991) to a more complex general circulation model (GCM) coupled to an ice-sheet model (Herrmann et al., 2003). Various data suggest that the lower CO₂ levels are probably the best estimates (Vandenbroucke et al., 2010; Pancost et al., 2013).

While the role of the pCO₂ has been fully investigated to explain the Ordovician glacial onset, much less has been done about the effect of ocean dynamics. Most of the modeling studies used atmospheric GCMs with slab mixed-layer oceans (Crowley et al., 1993; Crowley and Baum, 1995; Gibbs et al., 1997; Herrmann et al., 2003, 2004b). In these models, there are no currents and the poleward oceanic heat transport (OHT) is prescribed through a diffusion coefficient. Ocean dynamics are thus totally absent. Herrmann et al. (2004b) gave a precursory approach of this question by varying the diffusion coefficient in the slab model GENESIS. They found that the paleogeographical evolution and a drop in CO₂ were only preconditioning factors that were not sufficient to trigger continental glaciation during the Hirnantian. Additional changes, such as a sea level drop or a lowered OHT, were necessary to cause glaciation within a reasonable CO₂ interval (8–18 PAL). These results suggest that the ocean played a fundamental role in the Late Ordovician. But, as Herrmann et al. (2004a) noticed later, the non-dynamical ocean precluded any further investigation. Only two studies have so far investigated the role of ocean dynamics in the Late Ordovician climate changes. Poussart et al. (1999) coupled an energy/moisture balance model to an oceanic GCM and to a sea ice model. They demonstrated that the Ordovician geographical configuration may have led to critical changes in the OHT with an up to 42% increase in the southern OHT compared to present-day and a highly asymmetric OHT relative to the equator. These results point out a major limitation inherent to the slab models that are unable to explicitly resolve the OHT, which value and pattern are often prescribed based on present-day estimates. Poussart et al. (1999) also demonstrated an important sensitivity of the oceanic overturning circulation to changes in the ice–snow

albedo parameters, highlighting a strong coupling between atmospheric and oceanic components in the Ordovician climate. Herrmann et al. (2004a) forced an oceanic GCM with the GENESIS outputs from Herrmann et al. (2004b). They compared the results obtained for two Late Ordovician paleogeographical configurations. They noticed a significant impact of paleogeography on the oceanic overturning circulation and thus on the OHT. They moreover demonstrated a non-linear response of the OHT to the atmospheric forcing. Their uncoupled modeling procedure however did not allow for any feedback of ocean on global climate.

Here, we propose to extend previous studies by exploring the effects of changed radiative forcing due to CO₂ decrease on the Ordovician climate with the fully coupled ocean–atmosphere–sea ice GCM FOAM (Jacob, 1997) for two time intervals within the Middle and Late Ordovician. The climate model and the experimental setup are briefly described in Sect. 2. The modeling results are presented in Sect. 3. The main results are discussed in Sect. 4 and a conclusion is given in Sect. 5.

## 2 Model description and experimental setup

### 2.1 The coupled climate model FOAM

We used the Fast Ocean Atmosphere Model (FOAM) version 1.5, a 3-D GCM with coupled components for ocean, atmosphere, and sea ice (Jacob, 1997). The atmospheric module is a parallelized version of the National Center for Atmospheric Research's (NCAR) Community Climate Model 2 (CCM2) with the upgraded radiative and hydrologic physics from CCM3 version 3.2 (Kiehl et al., 1998). It was run at a R15 spectral resolution (4.5° × 7.5°) with 18 vertical levels. The ocean module, the Ocean Model version 3 (OM3), is a 24-level z coordinate ocean GCM providing a 1.4° × 2.8° resolution. The sea ice module uses the thermodynamic component of the CSM1.4 sea ice model, which is based on the Semtner 3-layer thermodynamic snow/ice model (Semtner, 1976). The coupled model, FOAM, is well designed for paleoclimate studies. It has no flux correction and its quick turnaround time allows for long millennium-scale integrations. FOAM has been widely used in paleoclimate studies (Poulsen and Jacob, 2004; Donnadieu et al., 2009; Zhang et al., 2010; Dera and Donnadieu, 2012; Lefebvre et al., 2012) including for the Ordovician (Nardin et al., 2011).

It is noteworthy that no land ice component was included. Correctly simulating deep-time ice sheets constitutes a whole study as such (see for example the review from Pollard, 2010) and it was therefore left for future work. Preliminary simulations conducted for the Late Ordovician (not shown here) however confirm that the climatic behavior demonstrated in this study persists in the presence of continental-sized ice sheets.
<figure><img src="image_71.png" /><figcaption>**Figure 1.** Continental configurations with the FOAM ocean grid based on the paleogeographical reconstructions from Blakey (2007), with the color code for the bathymetric and topographic categories used in this study. G: Gondwana, A: Avalonia, B: Baltica, L: Laurentia, S: Siberia. Ocean names are shown. Each time slice is replaced on the synthetic Ordovician glaciation chronology (inspired from Finnegan et al., 2011). Age limits are taken from the geological time scale v2014/02 (International Commission on Stratigraphy, 2014). Question marks are time intervals when hypothetical first ice sheets may have appeared according to Turner et al. (2012). Blue bars represent the two scenarios considered today for the duration of the Ordovician glaciation (see Sect. 2.2. for a full description).</figcaption></figure>

## 2.2 Experimental setup

Spjeldnaes (1962), working on floro-faunal data, was the first to postulate a glaciation during the Ordovician. Since then, numerous studies widely documented the Ordovician glacial sedimentary record in North Africa (e.g., Beuf et al., 1971; Denis et al., 2007; Ghienne et al., 2007), South Africa (e.g., Young et al., 2004) and South America (e.g., Díaz-Martínez and Grahn, 2007). Dating these deposits remains highly difficult (e.g., Díaz-Martínez and Grahn, 2007) and geochemical methods bring conflicting results (see for example Trotter et al., 2008 and Finnegan et al., 2011). As a consequence, the duration of the glaciation is still not well constrained. Two models are still debated: a short-lived (< 2 Myr) glaciation limited to the Late Ordovician Hirnantian (445–443 Ma) (Brenchley et al., 1994, 2003; Sutcliffe et al., 2000), or a long-lived glaciation (> 20 Myr) that extended from the Late Ordovician Katian to the Silurian Wenlock, with the Hirnantian only representing the glacial maximum (Saltzman and Young, 2005; Loi et al., 2010; Finnegan et al., 2011) (Fig. 1). Some authors even propose a Darriwilian ice age based on eustatic cyclicities (Turner et al., 2012).

In this study, simulations were carried out for two paleogeographical configurations reconstructed by Blakey (2007), one for 470 Ma (Middle Ordovician Dapingian) and one for 450 Ma (Late Ordovician Katian) (Fig. 1). These two reconstructions have been chosen to be respectively representative of the following: (i) the pre-glacial interval; (ii) the undisputable Late Ordovician glacial interval. The 450 Ma reconstruction is considered here as a good estimate for the Late Ordovician Hirnantian (445 Ma). The paleogeographic

changes are limited at that timescale and they are assumed to not critically impact global climate. Since previous studies have revealed the major role played by the topography to initiate glacial conditions during the Ordovician (Crowley and Baum, 1995), we defined five realistic topographic classes (Fig. 1). The ocean model bathymetry is a flat-bottom case. As Ordovician vegetation was limited to non-vascular plants (Steemans et al., 2009; Rubinstein et al., 2010) – the coverage of these can hardly be estimated – we imposed a bare soil (rocky desert) at every land grid point. No ice sheet was prescribed on the land. The solar luminosity was set 3.5% below present according to the model of solar physics provided by Gough (1981), and we used the “cold summer orbital configuration” (CSO). The CSO is defined with an obliquity of 22.5°, an eccentricity of 0.05 and a longitude of perihelion of 270°. Because most explanatory hypotheses require a major $pCO_2$ fall to initiate glacial conditions in the Middle–Late Ordovician (Brenchley et al., 1994; Kump et al., 1999; Villas et al., 2002; Lenton et al., 2012), we finally used five $pCO_2$ values representative for presumed Ordovician conditions: 4, 6, 8, 12 and 16 PAL. A lower value, 2 PAL, was also tested in order to better understand the climatic behavior although this $pCO_2$ clearly stands beyond the lower bound of the $pCO_2$ range defined for the Ordovician Period. For each simulation, the model was run until deep ocean equilibrium was reached, which means 2000 years of integration except for the 450 Ma experiment with the $pCO_2$ set to 6 PAL that required 1000 years more.
<figure><img src="image_72.png" /></figure>

**Figure 2.** (a) Mean annual, globally averaged surface air temperature for various *pCO*₂ values from 2 to 16 PAL, with the slab-version of FOAM and with the fully coupled ocean–atmosphere–sea ice version of FOAM at 470 and at 450 Ma. The sketch illustrates the explanations provided in the main text. Note the logarithmic scale on the *x* axis. (b) Same plot as (a), for the mean annual global sea ice cover expressed in M km⁻².

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

**Figure 3.** SST (color shading) and permanent sea ice cover (blue-white shading) simulated at 450 Ma with the fully coupled version of FOAM within the cold climatic mode (at 6 PAL) and within the warm climatic mode (at 8 PAL). The permanent sea ice cover is defined as the sum of all grid points whose at least 80 % of the surface is covered by sea ice during a year, allowing to discard temporary sea ice cover by icebergs or rapidly advancing and retreating ice fronts. The gray mask stands for continental masses.

# 3 Climate simulation results

## 3.1 Sensitivity to an atmospheric CO₂ decrease

Figure 2 presents the relationship between the global mean surface air temperature (SAT) and the *pCO*₂ for all the simulations discussed in this paper. For comparison with the previous studies undertaken with slab oceans (Crowley et al., 1993; Crowley and Baum, 1995; Gibbs et al., 1997; Herrmann et al., 2003, 2004b), we first ran FOAM with a 50 m

deep mixed-layer slab ocean (black dashed lines in Fig. 2). All boundary conditions were kept the same as for the runs with the ocean–atmosphere GCM (see Sect. 2.2) and we followed the usual conventions by setting the diffusion coefficient so that total OHT reaches its present value. As expected from previous studies, the Late Ordovician (450 Ma) climate system is sensitive to a *pCO*₂ decrease and the temperature–CO₂ relationship is linear in a log-axis graph (Fig. 2a). The model sensitivity, 3 °C per halving of CO₂, lies very close to the value 2.5 °C obtained by Herrmann et al. (2003). It is however noteworthy that, for an identical atmospheric forcing value, climate is systematically warmer compared to GENESIS (Gibbs et al., 1997; Herrmann et al., 2004b). An immediate explanation lies in the prescribed solar constant. We chose a solar constant 3.5 % weaker than today following the parametric law from Gough (1981). Previous studies mostly used a 4.5 % depressed insolation compared to present-day (Crowley and Baum, 1995; Gibbs et al., 1997; Herrmann et al., 2003, 2004b), following the constant rate of decrease of 1 % every hundred million years proposed by Crowley and Baum (1992). Both values are within the range –3.5 to –5 % defined by Endal and Sofia (1981), but our weaker decrease of the solar constant logically results in higher temperatures. We must note that modeled absolute temperatures are also sensitive to ice albedo and cloud parameterizations and are thus highly model-dependent (see for example Braconnot et al., 2012). Proposing absolute climate estimates for the Ordovician remains highly challenging and lies far beyond the main purpose of our study. Aside that absolute temperature offset, FOAM produces a trend similar to results previously published when it is used as a slab.

Using the same (450 Ma) paleogeographical configuration, FOAM experiments with ocean dynamics lead to a much different climatic response to the CO₂ forcing. The temperature–CO₂ relationship is non-linear (Fig. 2a). Two
<figure><img src="image_74.png" /><figcaption>**Figure 4.** Mean annual, zonally averaged SAT simulated for each *pCO*₂ at 450 Ma with the fully coupled version of FOAM.</figcaption></figure>

climatic modes are simulated as forcing parameters slowly
change.

- For high CO₂ concentrations (≥ 8 PAL), the warm mode is marked by a hot climate (Figs. 2a and 3) and no sea ice in the Northern Hemisphere (NH) (Figs. 2b and 3). The model sensitivity is similar to the one obtained with the slab (3 °C per halving of CO₂) and lies within the range 2.1-4.7 °C defined with various up-to-date GCMs for present-day simulations (IPCC, 2013). The temperature decrease associated to a CO₂ drop is almost evenly distributed at all latitudes (Fig. 4).

- When atmospheric forcing is decreased from 8 to 4 PAL, the annual global mean SAT falls by 8.65 °C from 8 to 6 PAL, and by 7.6 °C from 6 to 4 PAL. The Ordovician climate abruptly enters a cool climatic mode characterized by a fast extension of the sea ice cover (Figs. 2, 3 and 5). A tipping point is reached, a slight *pCO₂* drop leading to a sharp cooling. The temperature decrease associated to a halving of CO₂ (from 8 to 4 PAL), -16 °C, lies far beyond the IPCC (2013) upper bound. The SAT drop is unevenly distributed against the latitudes. The minimum mean annual cooling is observed in the tropics (up to -5 °C from 8 to 6 PAL) whereas the maximum cooling (up to -27 °C from 8 to 6 PAL) occurs in both hemispheres at polar latitudes (Fig. 4). The amplification of the thermal response at high latitudes is due to the sea ice albedo positive feedback. The hemispheric sea ice cover (Table 1, Fig. 5b and c) suggests in particular a strong relationship between the NH sea ice extension and the global climate state. No ice is present in the NH within the warm mode (Fig. 3), and it is the expansion of the sea ice from the North Pole that marks the shift to the cold climatic mode. It is noteworthy here that the simulation of the sea ice in the model FOAM has been proven correct enough in the modern climate to be reasonably extended to paleoclimate studies. The reader is invited to

<figure><img src="image_75.png" /><figcaption>**Figure 5.** Zonally averaged sea ice cover fraction (expressed in percents) simulated for each *pCO₂* at 450 Ma with the fully coupled version of FOAM. DJF: December, January and February. JJA: June, July and August.</figcaption></figure>

refer to the study from Lefebvre et al. (2012, paragraph
14) for the whole analysis.

- From 4 to 2 PAL, the Earth's climate is still in cold mode but the SAT drop tends to slow down. Another threshold is crossed and the model sensitivity decreases back to a level characteristic of glacial climate (~5 °C per CO₂ halving, see Berner and Kothavala, 2001). The climate change does not differ so much from the cooling occurring within the warm climatic mode between 16 and 8 PAL for a corresponding halving of CO₂, though glacial conditions still provide a higher global sensitivity (Kothavala et al., 1999) and an enhanced polar cooling. Sea ice continues to grow almost with the same intensity in the Southern Hemisphere (SH) but its spread is very limited in the NH (Table 1), where the sea ice edge stabilizes near 40° N (Fig. 5a).

The simulations were replicated for the 470 Ma continen-
tal configuration. A similar evolution is observed. With the
slab model (Fig. 2), the temperature-CO₂ relationship is still
linear in the log-axis graph. For the coupled runs (Fig. 2 and
Table 1), two phases of stable climate are separated by a ma-
jor climatic instability within the interval 12–6 PAL. The tip-
ping point showing the shift from the warm to the cold cli-
matic mode is reached this time between 12 and 8 PAL and
is still associated to a sudden extension of the sea ice in the
NH. Three main conclusions follow from theses results: (i)
the climatic instability originates in ocean dynamics, as it is
not observed with the slab, (ii) NH sea ice cover plays a crit-
ical role in the shift to the cold climatic mode, (iii) reaching
the cold climatic mode requires a lower CO₂ level at 450 Ma
than at 470 Ma.
<table><caption>**Table 1.** Overview of the coupled simulations. The name of each experiment includes the age of the paleogeography used in input followed by the *pCO*₂ expressed in PAL. These abbreviations are used throughout the manuscript. Slab runs are not shown.</caption><thead><tr><th>Simulations</th><th><i>pCO</i><sub>2</sub> (ppm)</th><th>Surface air temperature (°C)</th><th>Global temperature of the Ocean (°C)</th><th>Mean annual global sea ice cover (10<sup>6</sup> km<sup>2</sup>)</th><th>Mean annual Southern Hemisphere sea ice cover (10<sup>6</sup> km<sup>2</sup>)</th><th>Mean annual Northern Hemisphere sea ice cover (10<sup>6</sup> km<sup>2</sup>)</th></tr></thead><tbody><tr><td>450 Ma-2X</td><td>560</td><td>-0.3</td><td>1.8</td><td>123.9</td><td>49.3</td><td>74.6</td></tr><tr><td>450 Ma-4X</td><td>1120</td><td>4.5</td><td>2.5</td><td>105.1</td><td>35.2</td><td>69.9</td></tr><tr><td>450 Ma-6X</td><td>1680</td><td>12.1</td><td>4.1</td><td>58.5</td><td>15.5</td><td>43.0</td></tr><tr><td>450 Ma-8X</td><td>2240</td><td>20.8</td><td>9.0</td><td>1.8</td><td>1.8</td><td>0.0</td></tr><tr><td>450 Ma-12X</td><td>3360</td><td>22.5</td><td>11.4</td><td>1.1</td><td>1.1</td><td>0.0</td></tr><tr><td>450 Ma-16X</td><td>4480</td><td>23.7</td><td>12.7</td><td>0.9</td><td>0.9</td><td>0.0</td></tr><tr><td>470 Ma-2X</td><td>560</td><td>4.7</td><td>2.3</td><td>102.3</td><td>28.3</td><td>74.0</td></tr><tr><td>470 Ma-4X</td><td>1120</td><td>8.6</td><td>3.1</td><td>85.4</td><td>19.2</td><td>66.2</td></tr><tr><td>470 Ma-6X</td><td>1680</td><td>11.8</td><td>3.9</td><td>63.1</td><td>10.9</td><td>52.2</td></tr><tr><td>470 Ma-8X</td><td>2240</td><td>16.7</td><td>6.7</td><td>26.0</td><td>3.2</td><td>22.8</td></tr><tr><td>470 Ma-12X</td><td>3360</td><td>22.5</td><td>11.1</td><td>0.3</td><td>0.3</td><td>0.0</td></tr><tr><td>470 Ma-16X</td><td>4480</td><td>23.5</td><td>12.0</td><td>0.4</td><td>0.4</td><td>0.0</td></tr></tbody></table>

<figure><img src="image_76.png" /><figcaption>**Figure 6.** Comparison of the SST simulated at 12 PAL with the fully coupled version of FOAM for the two continental configurations used in this study. The black mask stands for continental masses. Note that, for the plot presented in (c), the mask actually represents all the grid points where the SST difference could not be calculated, cumulating each point that is a continent either in the 450 or in the 470 Ma configuration.</figcaption></figure>

This type of non-linearity within the Earth climate system has been extensively studied in the past few years. We will in particular discuss our results in the light of the major studies from Rose and Marshall (2009) and Ferreira et al. (2011) in

Sect. 4, but we propose for now to continue our analysis by focusing on the response of ocean dynamics to the paleogeographical changes occurring between 470 and 450 Ma.

## 3.2 Sensitivity of ocean dynamics to the paleogeography

The shift from the warm to the cold climatic mode occurs between 12 and 8 PAL at 470 Ma, and between 8 and 6 PAL at 450 Ma. In order to constrain why the 450 Ma climatic system has a lower CO₂ threshold, we compare the climatic features at both time slices for the same *pCO*₂. We use the simulations preceding the tipping points (12 PAL) to identify the preconditioning factors for the climatic instability. These two simulations have a same global SAT (22.5 °C, Table 1) and a close global ocean temperature (11.1 °C at 470 Ma and 11.4 °C at 450 Ma, Table 1) that make the comparison reasonable.

Figure 6 displays the sea-surface temperatures (SST) for the 470 Ma-12X and the 450 Ma-12X experiments together with the SST difference between these two simulations. From 470 to 450 Ma, the SH cools and the NH warms (Fig. 6c). The temperature changes are on the order of a few degrees, locally up to 3 °C in both hemispheres. Assuming that the NH sea ice extension plays a major role in the abrupt transition from one climatic state to another through the sea ice albedo positive feedback (see Sect. 3.1), the NH warming observed in Fig. 6c appears as the key point accounting for the increasing stability from 470 to 450 Ma. It becomes more difficult to form sea ice and thus to trigger the albedo feedback in the NH through time, and a lower CO₂ concentration is required to enter the cool climatic mode at 450 Ma.

The origin of the temperature pattern observed in Fig. 6c is investigated through the study of the oceanic heat transport (Fig. 7). From 470 to 450 Ma, the NH is marked by
<figure><img src="image_77.png" /></figure>

**Figure 7.** Comparison of the total, advective and diffusive OHT simulated with the fully coupled version of FOAM at 12 PAL for the two continental configurations used in this study. The OHT is given in petawatts (1 PW = $10^{15}$ W). The present-day curve is the control experiment conducted at present-day conditions with the up-to-date IPSL (Institut Pierre Simon Laplace) CM4 GCM (see the main text for explanations).

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

**Figure 8.** Siberia's northward drift and the consequences on the longitudinal stream function. Eastward stream is positive (red), westward stream is negative (blue). The stream function is integrated over the whole water column. S: Siberia. Simulations were conducted with the fully coupled version of FOAM. It is noteworthy that the northward drift of Siberia from the Middle to the Late Ordovician is a robust feature of all available Ordovician paleogeographical reconstructions published so far, including the maps from Torsvik and Cocks (Cocks and Torsvik, 2007, Fig. 6a; Torsvik and Cocks, 2013: Fig. 2.12, 2.14 and 2.15), Golonka and Gaweda (2012, Figs. 6 and 7) and Scotese and McKerrow (1991, Figs. 4 and 5).

a significant increase in the northward OHT exceeding 100% at 20° N (from 0.4 to 0.8 PW respectively, Fig. 7a). In the SH, a slight decrease in the southward OHT is observed at all latitudes except for a narrow latitudinal band between 6 and 18° S (Fig. 7a). This new energy redistribution brings more heat to the NH and less to the SH. As a direct consequence, the NH warms from 470 to 450 Ma, and the SH cools.

Dividing the total OHT into its advective and diffusive components (Fig. 7b and c) reveals a striking feature within the NH advective heat transport between ~35 and ~65° N (Fig. 7b). The transport is negative, meaning a southward heat transport, which is a priori unexpected in the NH. The advective heat transport scales as ~$\Psi \times \Delta T \times \rho \times C_p$ where $\Psi$ is the strength of the overturning circulation, $\Delta T$ the temperature difference between the northward and

the southward limbs of the overturning circulation, $\rho$ the water density and $C_p$ the water specific heat. A negative heat transport implies a southward transport of surface warm waters. The southward water flow results from an intense wind-driven eastward circumpolar current that develops at these latitudes in the absence of landmasses (Fig. 8). As a consequence of the Ekman deflection, the water masses are deviated to the right of the main eastward stream and thus constitute a southward mass transport. From 470 to 450 Ma, the northward drift of Siberia (Fig. 8) inhibits the zonal current and the associated Ekman cell, explaining why the reverse OHT weakens through time (Fig. 7b). A similar reverse advective OHT is observed in present-day experiments in the SH due to the Antarctica Circumpolar Current. Figure 7b shows the results obtained for a present-day experiment conducted with a state-of-the-art ocean-atmosphere-vegetation coupled model: the IPSL CM4 GCM (Marti et al., 2010). A reverse heat transport is obtained near 40° S. From this point of view, the Ordovician quasi-oceanic NH keeps some common features with our present-day ~80% oceanic SH.

The global Northern Hemisphere OHT strengthening, from 470 to 450 Ma (Fig. 7a), is due to an increase of the transport in three main areas (Fig. 9), the location of which turns out to be constrained by changes in the overturning circulation (Fig. 10).

- Between 5 and 30° N, the advective OHT increases (Fig. 7b) in response to the intensifying surface subtropical cell that directly brings heat from the equator.

- Between 30 and ~42° N, the diffusive OHT increases (Fig. 7c) as a consequence of the strengthening of the subtropical cell. Much more heat is “deposited” at 30° N, which imposes a strong temperature gradient. Heat diffuses from the warm tropical latitudes where the heat is “deposited” to the cold waters at higher latitudes.

- Between ~42 and 60° N, the southward advective OHT weakens through time (Fig. 7b). The northward drift of Siberia inhibits the 40° N zonal current (Fig. 8): the
counter-clockwise cell associated to the zonal current acts over the entire water column at 470 Ma and not beyond 2500 m at 450 Ma. Hence less cold waters are brought from the high latitudes to the mid-latitudes resulting in an increased net northward heat transport. An additional process is the northward deflection of surface warm currents by Siberia at 450 Ma. This can be seen by comparing the isotherms 20 and 15 °C in Fig. 6a and b. Although it is restricted to upper levels oceanic waters, this second mechanism is significant because the first tens of meters represent an important heat reservoir of the ocean.

The most significant outstanding issue is to determine why the temperature drawdown occurring through the tipping point is sharper at 450 than at 470 Ma. Although the tipping point occurs at lower CO₂ at 450 Ma, the two paleogeographical configurations present almost the same global SAT and an equivalent sea ice cover at 6 PAL (Fig. 2 and Table 1). Paradoxically, this sharper response through the 450 Ma tipping point is due to the paleogeographically induced increase in the Northern Hemisphere OHT that previously inhibited the cooling. While the NH warms from 470 to 450 Ma (Fig. 6c), the NH equator-to-pole temperature profile flattens (not shown here). These changes are observed within the warm mode at 12 PAL, but they are critical since they constitute the preconditioning factors for the differential climatic sensitivity as climate shifts to the cold climatic mode. The flattened 450 Ma SST gradient delays the tipping point to a certain extent by inhibiting sea ice growth at polar latitudes, but it increases the global climate sensitivity to the radiative forcing as soon as NH sea ice begins to expand.

## 4 Discussion

The potential existence of multiple equilibria in the climate system has a long story in the literature. Lorenz (1968) discussed whether the physical laws governing climate are “transitive”, supporting a unique climate state, or “intransitive” leading to several sets of long-term statistics. Two main mechanisms supporting climate intransitivity have so far emerged. The first mechanism originates in the study of Stommel (1961) who showed that the thermohaline circulation may sustain two stable regimes depending on whether the driver of the circulation, the density difference, is dominated by the temperature or the haline gradient, leading respectively to an intense or a weak circulation. This concept motivated later a wealth of studies linking the last glacial abrupt events to these multiple thermohaline states, studies that are today well known as “fresh water hosing” experiments (see for example Kageyama et al., 2013). The second mechanism is the ice albedo positive feedback. It was independently introduced by Budyko (1969) and Sellers (1969). The authors developed simple one-dimension climate models in which surface temperature is computed by taking into

account the local radiative budget and the meridional heat transport. Two stable climate states are obtained for a same solar forcing depending on the initial conditions, describing a classical hysteresis loop: a warm state with a small sea ice cap or no sea ice at all (Eq3, Fig. 11a), and a cold completely ice-covered state (Eq1, Fig. 11a). A large sea ice cap of finite size extending to the mid-latitudes constitutes another analytical solution, but this third equilibrium turns out to be very unstable and thus not realizable (Eq2 in Fig. 11a; Rose and Marshall, 2009). Critical non-linearities are observed at both latitudinal extremes (Rose and Marshall, 2009 and references therein). At low latitudes, near the equator, the “large ice cap instability” rules. An sea ice margin reaching these latitudes expands up to the equator (“snowball” state Eq1, Fig. 11a). At very high latitudes, this is the “small ice cap instability” (SICI): a reduced sea ice cover will either completely disappear or spread until it forms a small ice cap restricted to the high latitudes (warm state, Eq3 in Fig. 11a). In the absence of the equilibrium Eq2, a unique sea ice cap of large (beyond the SICI domain) but finite (no snowball state) size is obtained for a given climatic forcing (Fig. 11b). A comprehensive review of the successive improvements associated to the energy balance models (EBMs) is provided by North et al. (1981).

Recently, Rose and Marshall (2009) managed to stabilize the equilibrium Eq2 by extending the previous EBMs through two key processes: (i) the effect by which sea ice insulates the ocean surface from the atmosphere, (ii) the wind-driven structure of the OHT. A primary maximum at low-latitudes accounts for the heat transport by the subtropical gyre and a secondary maximum at high latitudes for the OHT generated by the subpolar gyre, with a local minimum between (Fig. 11c). Contrary to classical EBMs in which the OHT describes a large parabolic curve from the equator to the pole (Fig. 11b), the OHT pattern from Rose and Marshall (2009) shows a local minimum at mid-latitudes where the sea ice edge can stably rest. This OHT structure is closer to the observations (see for example Garnier et al., 2000). Within the improved EBM of Rose and Marshall (2009), the small ice cap instability does not exist anymore: sea ice

<figure><img src="image_79.png" /><figcaption>**Figure 9.** Difference in total, advective and diffusive OHT between the 450 Ma-12X and the 470 Ma-12X simulations.</figcaption></figure>
<figure><img src="image_80.png" /><figcaption>**Figure 10.** Comparison of the global meridional overturning circulation simulated at 12 PAL with the fully coupled version of FOAM for the two continental configurations used in this study. Red shading and solid contours indicate a clockwise circulation, blue shading and dashed contours indicate a counter-clockwise circulation.</figcaption></figure>

caps of various size can persist at polar latitudes. Because the equilibrium Eq2 is now stable, two sea ice caps of large but finite size can be sustained for the same climatic forcing (e.g., for F5 in Fig. 11c). The position of the sea ice edge in the equilibrium Eq2 is captured within the position of the mid-latitude OHT minimum and thus very stable when external forcing is varied.

The question however still remained if the mid-latitude equilibrium Eq2 would hold in a more complex climate model with many degrees of freedom and a strong internal variability. Ferreira et al. (2011) demonstrated multiple equilibria using an up-to-date ocean GCM (the MITgcm) on two idealized configurations: a purely aqua planet, and an ocean planet with a strip of land extending from one pole to the other. The three states from Rose and Marshall (2009) are obtained: (i) a warm, ice-free state, (ii) a cold state with sea ice extending down to mid-latitudes, (iii) and the completely ice-covered “snowball” state. The authors highlight the consistency of their results with the study of Rose and Marshall (2009). In their complex GCM, the multiple equilibria owe their existence to the latitudinal structure of the OHT. The authors moreover extend the previous results from Rose and Marshall (2009) by demonstrating that no OHT minimum is required. A weak high-latitude OHT within the sub-polar gyre is sufficient to stabilize a finite sea ice cap at mid-latitudes, in that critical weak-OHT zone situated between

the subtropical and the subpolar gyres that they widely refer to as the “OHT convergence”.

Here, we have investigated the Ordovician climate sensitivity to a varying radiative forcing. The hallmark of our results is the existence of a global climatic instability in the NH. A tipping point separates two climatic modes: a warm mode with no sea ice in the NH, and a cold mode with sea ice extending down to the mid-latitudes (Fig. 5a). In the light of the results from Rose and Marshall (2009) and Ferreira et al. (2011), these climatic modes respectively correspond to the classical warm mode obtained with an EBM and to the third equilibrium previously described, respectively Eq3 and Eq2 in Fig. 11. The simulations marked by a weak *pCO₂* stabilize in Eq2, whereas the experiments conducted with high *CO₂* levels lead to the warm state of Eq3.

It is noteworthy that the latitudinal position at which the sea ice edge stabilizes within the cold climatic mode effectively lies immediately poleward of the Northern Hemisphere OHT convergence (Fig. 12, ~40° N), which is a key prediction from Rose and Marshall (2009) and Ferreira et al. (2011). According to these authors, the strong low-latitude OHT encounters here a weaker transport and heat is “deposited”, thus preventing the sea ice edge from moving further equatorward.

The fact that the ice edge is strongly enslaved to the OHT convergence (Fig. 12; Rose and Marshall, 2009) explains why the model returns to its initial sensitivity below 4 PAL (see Sect. 3.1). Ice cannot spread further and stabilizes. The sea ice albedo feedback is geometrically stopped by the subtropical high-OHT barrier and cooling is consequently slowed-down.

The comparison of Fig. 11b and c permits to explain why tipping points were not previously observed for the Ordovician. Similarly to the classical EBMs, the slab models use a purely diffusive transport. The resulting OHT describes the same large parabolic curve from the equator to the pole (Fig. 12). There is no OHT minimum or more largely no OHT convergence where the sea ice edge could rest to make the mid-latitude equilibrium stable. Climate cannot switch between stable states which inhibits the climatic instability.

In addition, the similarity between EBMs and slab models explains the results from Gibbs et al. (1997) who observed with the slab model GENESIS a runaway icehouse that they interpreted as resulting from the small ice cap instability. As it was previously explained, this mechanism is not observed anymore when the OHT structure is taken into account, but it is a classical feature of the EBMs and by extension of all models that use a diffusive OHT. We emphasize that the runaway icehouse proposed by Gibbs et al. (1997) and the climatic instability suggested in this study rely on much different mechanisms.

The absence of ocean dynamics in the slab models makes their response to an atmospheric forcing much different from the ocean-atmosphere GCM when simulating time periods of particular land-sea configuration as it is the case for the
<figure><img src="image_81.png" /></figure>

**Figure 11.** (a) Graph showing the latitude of the sea ice front (x axis) vs. the climatic forcing (y axis) for a classical Budyko–Sellers type EBM. Dark arrows indicate the pathway described by the model through a whole cycle of forcing decrease and increase. Red arrows account for major climatic instabilities. Red circles represent the sea ice edge latitudes corresponding to each analytical solution for a given climatic forcing F1. Dashed lines represent regions that are correct analytical solutions, but are unable to be realized within the model due to instability. SICI: small ice cap instability; LICI: large ice cap instability. (b) represents the same plot as (a), while focusing on the domain shaded in grey in (a). The forcing is decreased and then increased without reaching the snowball equilibrium Eq1. Several forcing values (F2, F3 and F4) are represented, each one leading to a single sea ice cap of finite size (red circles). The OHT pattern corresponding to the purely diffusive transport characteristic of classical EBMs is overlayed in grey (right y axis). (c) Same plot as (b) for the extended EBM from Rose and Marshall (2009) (see the main text for explanations) in which the OHT pattern is modified (grey curve). Note the two sea ice caps of finite size corresponding to the forcing level F5. (a) is adapted from North et al. (1981), (b) and (c) are modified after Rose and Marshall (2009).

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

**Figure 12.** Comparison of the total OHT at 450 Ma and at 8 PAL for the slab version of FOAM and for the fully coupled version of FOAM. The OHT pattern is relatively stable when CO₂ is varied. Color bars represent the maximum (boreal winter) Northern Hemisphere sea ice extent simulated within the cold climatic mode at 450 Ma with the fully coupled model.

Ordovician. An efficient way to overcome that apparent deficiency would be to impose the structure of the latitudinal OHT, as Rose and Marshall (2009) did in their extended EBM. The decisive answer is brought by Ferreira et al. (2011), who tested this hypothesis within the MIT-gcm used as a slab. With the simple diffusive coefficient, the mid-latitude equilibrium is not stable, as for classical EBMs. However, by imposing the OHT deduced from their coupled runs to the slab, they managed to reproduce exactly the same climatic response as with their complex GCM. These results confirm that the differential behavior of ocean–atmosphere and slab–atmosphere GCMs lies within ocean dynamics through the OHT.

The Ordovician Period turned out to bring a favorable paleogeographical context to propose the sudden shift between different climate states on a realistic configuration. The stabilization of a sea ice edge at mid-latitudes requires a wind-driven latitudinal OHT structure, a condition that can be reached only in an oceanic context. It is notably the case of the Ordovician 95% oceanic NH which ocean zonal circulation is not hampered by continental masses. If continents were largely distributed at all latitudes, then continental edges would strongly alter the pattern originally imposed by the winds and by the Ekman pumping mechanism. The OHT may thus be more even at all latitudes, as it is the case in our study in the SH due to the supercontinent Gondwana (Figs. 1 and 12). No OHT convergence would be obtained and the OHT pattern would become very similar to the one imposed by a simple diffusion coefficient. From this point of view, it is logical that the global instability lies in our study in the oceanic NH, and it would anyway be inconceivable to obtain a similar behavior in the SH.
## 5 Conclusions

We have examined the response of the Ordovician climate to a decreasing $pCO_2$ for two continental configurations at 470 and at 450 Ma with the coupled ocean, atmosphere and sea ice GCM FOAM. We have found the existence of tipping points in the Ordovician climate owing to the particular paleogeography with a quasi-oceanic Northern Hemisphere. As a result, the climate may have abruptly shifted from a warm climatic mode to a cold climatic mode. Warm SST and no sea ice in the NH characterize the warm mode, whereas much colder SST and a large sea ice extension to the mid-latitudes mark the cool mode. We also note that subtle changes occurring in the continental configuration between 470 and 450 Ma make the latter less sensitive to the radiative forcing. Indeed, the shift from the warm to the cold climatic mode takes place between 12 and 8 PAL at 470 Ma, and between 8 and 6 PAL at 450 Ma. One potential mechanism relies on the changes occurring in the oceanic heat transport between the two simulated time periods. The paleogeographical evolution promotes a larger northward oceanic heat transport at 450 Ma which in turn induces warmer Northern Hemisphere SST. Spread of sea ice becomes more difficult and a lower $CO_2$ level is required to enter the cold climatic mode. The fact that the tipping point is reached for a higher $CO_2$ level at 470 Ma whereas climate was warm at that time (e.g., Trotter et al., 2008), implies that the $pCO_2$ must have been high during the Middle Ordovician, thus precluding climate to reach the climatic instability.

The tipping point is due to a climatic instability occurring in the quasi-oceanic Northern Hemisphere. The instability relies on the sudden shift from a warm climate to a climatic equilibrium characterized by a sea ice front resting at mid-latitudes. This particular climate state has been widely discussed in the literature and proven to require a wind-driven latitudinal structure of the OHT (Rose and Marshall, 2009; Ferreira et al., 2011). The classical EBMs and slab models use a purely diffusive OHT and are by construction unable to stabilize an ice edge at the mid-latitudes (Rose and Marshall, 2009; Ferreira et al., 2011), explaining why that behavior was never noticed in previous studies about the Ordovician climate (e.g., Gibbs et al., 1997; Herrmann et al., 2004b). Published studies performed with idealized land-sea configurations however demonstrate that this limitation may be overcome by imposing a latitudinally varied diffusion coefficient reproducing the OHT structure. Such extended EBMs and slab models successfully reproduce the behavior of the more complex ocean-atmosphere GCMs (Rose and Marshall, 2009; Ferreira et al., 2011).

The Late Ordovician Hirnantian is characterized by the rapid growth of large ice sheets (e.g., Ghienne et al., 2007; Loi et al., 2010; Finnegan et al., 2011). Recent geochemical data suggest an abrupt tropical SST cooling at that time (Trotter et al., 2008; Finnegan et al., 2011). Several $CO_2$ sinks have been proposed to explain cooling but no consen-

sus has been found so far (Brenchley et al., 1994; Kump et al., 1999; Villas et al., 2002; Lenton et al., 2012). In this context, the tipping point becomes a very interesting working hypothesis. It requires only a slight $CO_2$ decrease and at the same time induces a sharp global temperature drawdown, in the same order of magnitude as the Hirnantian cooling observed by Trotter et al. (2008) and Finnegan et al. (2011). The climatic instability demonstrated in this study therefore represents a step forward in understanding the Ordovician glacial events and provides in particular a minimal hypothesis for the Hirnantian abrupt cooling and associated glacial pulse. Our results also shed a new light on the $CO_2$ sinks envisaged to reduce the $pCO_2$ at that time. We now suggest that a moderate sink (e.g., Young et al., 2009; Nardin et al., 2011) may have been sufficient, thus making relevant some mechanisms that have been neglected so far because of their moderate impact on the $pCO_2$.

*Acknowledgements.* We thank Mark Williams and one anonymous referee for their constructive remarks that have greatly improved this paper. We thank Pierre Sepulchre for providing the IPSL CM4 GCM run. We are also grateful to David G. Ferreira for his constructive and helpful comments. This work was granted access to the HPC resources of TGCC under the allocation 2014-012212 made by GENCI.

Edited by: A. Haywood

## References

* Berner, R. A.: GEOCARBSULF: a combined model for Phanerozoic atmospheric $O_2$ and $CO_2$, Geochim. Cosmochim. Ac., 70, 5653–5664, doi:10.1016/j.gca.2005.11.032, 2006.

* Berner, R. A. and Kothavala, Z.: GEOCARB III: a revised model of atmospheric $CO_2$ over Phanerozoic time, Am. J. Sci., 301, 182–204, doi:10.2475/ajs.301.2.182, 2001.

* Beuf, S., Biju-Duval, B., De Charpal, O., Rognon, P., Gariel, O., and Bennacef, A.: Les grès du Paléozoïque Inférieur au Sahara, Technip, Paris, 464 pp., 1971.

* Blakey, R. C.: Carboniferous-Permian paleogeography of the assembly of Pangea, in: Fifteenth International Congress on Carboniferous and Permian Stratigraphy, edited by: Wong, T. E., Utrecht, the Netherlands, pages 443–465, Royal Netherlands Academy of Arts and Sciences, Amsterdam, 2007.

* Braconnot, P., Harrison, S. P., Kageyama, M., Bartlein, P. J., Masson-Delmotte, V., Abe-Ouchi, A., Otto-Bliesner, B., and Zhao, Y.: Evaluation of climate models using palaeoclimatic data, Nat. Clim. Change, 2, 417–424, doi:10.1038/NCLIMATE1456, 2012.

* Brenchley, P. J., Marshall, J. D., Carden, G., Robertson, D., Long, D., Meidla, T., Hints, L., and Anderson, T. F.: Bathy-metric and isotopic evidence for a short-lived Late Ordovician glaciation in a greenhouse period, Geology, 22, 295–298, doi:10.1130/0091-7613(1994)022<0295:BAIEFA>2.3.CO;2, 1994.

* Brenchley, P. J., Carden, G. A., Hints, L., Kaljo, D., Marshall, J. D., Martma, T., Meidla, T., and Nölvak, J.: High-
resolution stable isotope stratigraphy of Upper Ordovician sequences: constraints on the timing of bioevents and environmental changes associated with mass extinction and glaciation, Geol. Soc. Am. Bull., 115, 89–104, doi:10.1130/0016-7606(2003)115<0089:HRSISO>2.0.CO;2, 2003.

Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619, doi:10.1111/j.2153-3490.1969.tb00466.x, 1969.

Cocks, L. R. M. and Torsvik, T. H.: Siberia, the wandering northern terrane, and its changing geography through the Palaeozoic, Earth-Sci. Rev., 82, 29–74, doi:10.1016/j.earscirev.2007.02.001, 2007.

Crowley, T. J. and Baum, S. K.: Toward reconciliation of Late Ordovician (~440 Ma) glaciation with very high CO₂ levels, J. Geophys. Res.-Atmos., 96, 22597–22610, doi:10.1029/91JD02449, 1991.

Crowley, T. J. and Baum, S. K.: Modeling late Palaeozoic glaciation, Geology, 2, 507–510, doi:10.1130/0091-7613(1992)020<0507:MLPG>2.3.CO;2, 1992.

Crowley, T. J. and Baum, S. K.: Reconciling Late Ordovician (440 Ma) glaciation with very high (14X) CO₂ levels, J. Geophys. Res.-Atmos., 100, 1093–1101, doi:10.1029/94JD02521, 1995.

Crowley, T. J., Baum, S. K., and Kim, K. Y.: General circulation model sensitivity experiments with pole-centered supercontinents, J. Geophys. Res.-Atmos., 98, 8793–8800, doi:10.1029/93JD00122, 1993.

Denis, M., Buoncristiani, J.-F., Konate, M., Ghienne, J.-F., and Guiraud, M.: Hirnantian glacial and deglacial record in SW Djado Basin (NE Niger), Geodin. Acta, 20, 177–195, doi:10.3166/ga.20.177-195, 2007.

Dera, G. and Donnadieu, Y.: Modeling evidences for global warming, Arctic seawater freshening, and sluggish oceanic circulation during the Early Toarcian anoxic event, Palaeoceanography, 27, PA2211, doi:10.1029/2012PA002283, 2012.

Díaz-Martínez, E. and Grahn, Y.: Early Silurian glaciation along the western margin of Gondwana (Peru, Bolivia and northern Argentina): Palaeogeographic and geodynamic setting, Palaeogeogr. Palaeocl., 245, 62–81, doi:10.1016/j.palaeo.2006.02.018, 2007.

Donnadieu, Y., Goddéris, Y., and Bouttes, N.: Exploring the climatic impact of the continental vegetation on the Mezosoic atmospheric CO₂ and climate history, Clim. Past, 5, 85–96, doi:10.5194/cp-5-85-2009, 2009.

Endal, A. S. and Sofia, S.: Rotation in solar-type stars, I – Evolutionary models for the spin-down of the sun, Astrophys. J., 243, 625–640, doi:10.1086/158628, 1981.

Ferreira, D., Marshall, J., and Rose, B.: Climate determinism revisited: multiple equilibria in a complex climate model, J. Climate, 24, 992–1012, doi:10.1175/2010JCLI3580.1, 2011.

Finnegan, S., Bergmann, K., Eiler, J. M., Jones, D. S., Fike, D. A., Eisenman, I., Hughes, N. C., Tripati, A. K., and Fischer, W. W.: The magnitude and duration of Late Ordovician-Early Silurian glaciation, Science, 331, 903–906, doi:10.1126/science.1200803, 2011.

Garnier, E., Barnier, B., Siefridt, L., and Béranger, K.: Investigating the 15 years air-sea flux climatology from the ECMWF re-analysis project as a surface boundary condition for ocean

models, Int. J. Climatol., 20, 1653–1673, doi:10.1002/1097-0088(20001130)20:14<1653::AID-JOC575>3.0.CO;2-G, 2000.

Ghienne, J.-F., Le Heron, D., Moreau, J., Denis, M., and Deynoux, M.: The Late Ordovician glacial sedimentary system of the North Gondwana platform, in: Glacial sedimentary processes and products, Special Publication, edited by: Hambrey, M., Christoffersen, P., Glasser, N., Janssen, P., Hubbard, B., and Siegert, M., 39, 295–319, International Association of Sedimentologists, Blackwells, Oxford, 2007.

Gibbs, M. T., Barron, E. J., and Kump, L. R.: An atmospheric $pCO_2$ threshold for glaciation in the Late Ordovician, Geology, 25, 447–450, doi:10.1130/0091-7613(1997)025<0447:AAPCTF>2.3.CO;2, 1997.

Golonka, J. and Gaweda, A.: Plate tectonic evolution of the southern margin of Laurussia in the Paleozoic, in: Tectonics – Recent advances, edited by: Sharkov, E., 261–282, InTech, 2012.

Gough, D. O.: Solar interior structure and luminosity variations, Sol. Phys., 74, 21–34, doi:10.1007/BF00151270, 1981.

Herrmann, A. D., Patzkowsky, M. E., and Pollard, D.: Obliquity forcing with 8–12 times preindustrial levels of atmospheric $pCO_2$ during the Late Ordovician glaciation, Geology, 31, 485–488, doi:10.1130/0091-7613(2003)031<0485:OFWTPL>2.0.CO;2, 2003.

Herrmann, A. D., Haupt, B. J., Patzkowsky, M. E., Seidov, D., and Slingerland, R. L.: Response of Late Ordovician palaeoceanography to changes in sea level, continental drift, and atmospheric $pCO_2$: potential causes for long-term cooling and glaciation, Palaeogeogr. Palaeocl., 210, 385–401, doi:10.1016/j.palaeo.2004.02.034, 2004a.

Herrmann, A. D., Patzkowsky, M. E., and Pollard, D.: The impact of paleogeography, $pCO_2$, poleward ocean heat transport and sea level change on global cooling during the Late Ordovician, Palaeogeogr. Palaeocl., 206, 59–74, doi:10.1016/j.palaeo.2003.12.019, 2004b.

International Commission on Stratigraphy: International Chronostatigraphic Chart v2014/02, available at: www.stratigraphy.org, 2014.

IPCC: Climate Change 2013: The physical science basis, in: Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 1535 pp., 2013.

Jacob, R. L.: Low frequency variability in a simulated atmosphere ocean system, Ph.D. thesis, University of Wisconsin-Madison, 1997.

Kageyama, M., Merkel, U., Otto-Bliesner, B., Prange, M., Abe-Ouchi, A., Lohmann, G., Ohgaito, R., Roche, D. M., Singarayer, J., Swingedouw, D., and Zhang, X.: Climatic impacts of fresh water hosing under Last Glacial Maximum conditions: a multi-model study, Clim. Past, 9, 935–953, doi:10.5194/cp-9-935-2013, 2013.

Kiehl, J. T., Hack, J. J., Bonan, G. B., Boville, B. A., Williamson, D. L., and Rasch, P. J.: The national center for Atmospheric Research community climate model: CCM3, J. Climate, 11, 1131–1149, doi:10.1175/1520-0442(1998)011<1131:TNCFAR>2.0.CO;2, 1998.
Kothavala, Z., Oglesby, R. J., and Saltzman, B.: Sensitivity of equilibrium surface temperature of CCM3 to systematic changes in atmospheric CO₂, Geophys. Res. Lett., 26, 209–212, doi:10.1029/1998GL900275, 1999.

Kump, L. R., Arthur, M. A., Patzkowsky, M. E., Gibbs, M. T., Pinkus, D. S., and Sheehan, P. M.: A weathering hypothesis for glaciation at high atmospheric pCO₂ during the Late Ordovician, Palaeogeogr. Palaeocl., 152, 173–187, doi:10.1016/S0031-0182(99)00046-2, 1999.

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

Lenton, T. M., Crouch, M., Johnson, M., Pires, N., and Dolan, L.: First plants cooled the Ordovician, Nat. Geosci., 5, 86–89, doi:10.1038/ngeo1390, 2012.

Loi, A., Ghienne, J.-F., Dabard, M. P., Paris, F., Botquelen, A., Christ, N., Elaouad-Debbaj, Z., Gorini, A., Vidal, M., Videt, B., and Destombes, J.: The Late Ordovician glacio-eustatic record from a high-latitude storm-dominated shelf succession: the Bou Ingarf section (Anti-Atlas, Southern Morocco), Palaeogeogr. Palaeocl., 296, 332–358, doi:10.1016/j.palaeo.2010.01.018, 2010.

Lorenz, E. N.: Climatic determinism, Meteor. Mon., 8, 1–3, 1968.

Marti, O., Braconnot, P., Dufresne, J. L., Bellier, J., Benshila, R., Bony, S., Brockmann, P., Cadule, P., Caubel, A., Codron, F., Noblet, N., Denvil, S., Fairhead, L., Fichefet, T., Foujols, M. A., Friedlingstein, P., Goosse, H., Grandpeix, J. Y., Guilyardi, E., Hourdin, F., Idelkadi, A., Kageyama, M., Krinner, G., Lévy, C., Madec, G., Mignot, J., Musat, I., Swingedouw, D., and Talandier, C.: Key features of the IPSL ocean atmosphere model and its sensitivity to atmospheric resolution, Clim. Dynam., 34, 1–26, doi:10.1007/s00382-009-0640-6, 2010.

Nardin, E., Goddéris, Y., Donnadieu, Y., Le Hir, G., Blakey, R. C., Pucéat, E., and Aretz, M.: Modeling the early Paleozoic long-term climatic trend, Geol. Soc. Am. Bull., 123, 1181–1192, doi:10.1130/B30364.1, 2011.

North, G. R., Cahalan, R. F., and Coakley, J. A.: Energy balance climate models, Rev. Geophys., 19, 91–121, 80R1502, doi:10.1029/RG019i001p00091, 1981.

Pancost, R. D., Freeman, K. H., Herrmann, A. D., Patzkowsky, M. E., Ainsaar, L., and Martma, T.: Reconstructing Late Ordovician carbon cycle variations, Geochim. Cosmochim. Ac., 105, 433–454, doi:10.1016/j.gca.2012.11.033, 2013.

Pollard, D.: A retrospective look at coupled ice sheet–climate modeling, Climate Change, 100, 173–194, doi:10.1007/s10584-010-9830-9, 2010.

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

Poussart, P. F., Weaver, A. J., and Barnes, C. R.: Late Ordovician glaciation under high atmospheric CO₂: a coupled model analysis, Paleoceanography, 14, 542–558, doi:10.1029/1999PA900021, 1999.

Rose, B. E. J. and Marshall, J.: Ocean heat transport, sea ice, and multiple climate states: insights from energy balance models, J. Atmos. Sci., 66, 2828–2843, doi:10.1175/2009JAS3039.1, 2009.

Rubinstein, C. V., Gerrienne, P., de la Puente, G. S., Astini, R. A., and Steemans, P.: Early Middle Ordovician evidence for land plants in Argentina (eastern Gondwana), New Phytol., 188, 365–369, doi:10.1111/j.1469-8137.2010.03433.x, 2010.

Saltzman, M. R. and Young, S. A.: Long-lived glaciation in the Late Ordovician? Isotopic and sequence-stratigraphic evidence from western Laurentia, Geology, 33, 109–112, doi:10.1130/G21219.1, 2005.

Scotese, C. R. and McKerrow, W. S.: Ordovician plate tectonic reconstructions, in: Advances in Ordovician geology, edited by: Barnes, C. R. and Williams, S. H., 90–9, 271–282, Geological Survey of Canada, 1991.

Sellers, W. D.: A global climatic model based on the energy balance of the earth-atmosphere system, J. Appl. Meteorol., 8, 392–400, doi:10.1175/1520-0450(1969)008<0392:AGCMBO>2.0.CO;2, 1969.

Semtner Jr, A. J.: A model for the thermodynamic growth of sea ice in numerical investigations of climate, J. Phys. Oceanogr., 6, 379–389, doi:10.1175/1520-0485(1976)006<0379:AMFTTG>2.0.CO;2, 1976.

Servais, T., Owen, A. W., Harper, D. A., Kröger, B., and Munnecke, A.: The Great Ordovician Biodiversification Event (GOBE): the palaeoecological dimension, Palaeogeogr. Palaeocl., 294, 99–119, doi:10.1016/j.palaeo.2010.05.031, 2010.

Sheehan, P. M.: The late Ordovician mass extinction, Annu. Rev. Earth Pl. Sc., 29, 331–364, doi:10.1146/annurev.earth.29.1.331, 2001.

Spjeldnaes, N.: Ordovician climatic zones, Norsk Geol. Tidsskr., 41, 45–77, 1962.

Steemans, P., Le Herissé, A., Melvin, J., Miller, M. A., Paris, F., Verniers, J., and Wellman, C. H.: Origin and radiation of the earliest vascular land plants, Science, 324, 353–353, doi:10.1126/science.1169659, 2009.

Stommel, H.: Thermohaline convection with two stable regimes of flow, Tellus, 13, 224–230, doi:10.1111/j.2153-3490.1961.tb00079.x, 1961.

Sutcliffe, O. E., Dowdeswell, J. A., Whittington, R. J., Theron, J. N., and Craig, J.: Calibrating the Late Ordovician glaciation and mass extinction by the eccentricity cycles of Earth's orbit, Geology, 28, 967–970, doi:10.1130/0091-7613(2000)28<967:CTLOGA>2.0.CO;2, 2000.

Torsvik, T. H. and Cocks, L. R. M.: New global palaeogeographical reconstructions for the Early Palaeozoic and their generation, in: Early Palaeozoic biogeography and palaeogeography, edited by: Harper, D. A. T. and Servais, T., 38, 5–24, Geological Society of London, Memoirs, 2013.

Trotter, J. A., Williams, I. S., Barnes, C. R., Lécuyer, C., and Nicoll, R. S.: Did cooling oceans trigger Ordovician biodiversification? Evidence from conodont thermometry, Science, 321, 550–554, doi:10.1126/science.1155814, 2008.

Turner, B. R., Armstrong, H. A., Wilson, C. R., and Makhlouf, I. M.: 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, doi:10.1016/j.sedgeo.2012.01.002, 2012.

Vandenbroucke, T. R., Armstrong, H. A., Williams, M., Paris, F., Zalasiewicz, J. A., Sabbe, K., Nölvak, J., Challands, T. J., Verniers, J., and Servais, T.: Polar front shift and atmospheric CO₂ during the glacial maximum of the Early Pa-
leozoic Icehouse. P. Natl. Acad. Sci., 107, 14983–14986, doi:10.1073/pnas.1003220107, 2010.

Villas, E., Vennin, E., Álvaro, J. J., Hammann, W., Herrera, Z. A., and Piovano, E. L.: The late Ordovician carbonate sedimentation as a major triggering factor of the Hirnantian glaciation, B. Soc. Geol. Fr., 173, 569–578, doi:10.2113/173.6.569, 2002.

Yapp, C. J. and Poths, H.: Ancient atmospheric CO₂ pres-
sures inferred from natural goethites, Nature, 355, 342–344,
doi:10.1038/355342a0, 1992.

Young, G. M., Minter, W., and Theron, J. N.: Geochemistry and palaeogeography of upper Ordovician glaciogenic sedimentary rocks in the Table Mountain Group, South Africa, Palaeogeogr. Palaeocl., 214, 323–345, doi:10.1016/j.palaeo.2004.07.029, 2004.

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

Zhang, Z.-S., Yan, Q., and Wang, H.-J.: Has the Drake passage played an essential role in the Cenozoic cooling?, Atmos. Ocean. Sci. Lett., 3, 288–292, 2010.
# The role of palaeogeography in the Phanerozoic history of atmospheric CO₂ and climate

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

CrossMark

Yves Goddéris a,*, Yannick Donnadieu b, Guillaume Le Hir c, Vincent Lefebvre a, Elise Nardin a

## ARTICLE INFO

Article history:

Received 27 June 2013

Accepted 9 November 2013

Available online 16 November 2013

Keywords:

Carbon cycle

Climate

Palaeogeography

Phanerozoic

Weathering

Modelling

## ABSTRACT

The role of the palaeogeography on the geological evolution of the global carbon cycle has been suspected since the development of the first global geochemical models in the early 80s. The palaeogeography has been rapidly recognized as a key factor controlling the long-term evolution of the atmospheric CO₂ through its capability of modulating the efficiency of the silicate weathering. First the role of the latitudinal position of the continents has been emphasized: an averaged low latitudinal position promotes the CO₂ consumption by silicate weathering, and is theoretically associated to low CO₂ periods. With the increase of model complexity and the explicit consideration of the hydrological cycle, the importance of the continentality factor has been recognized: periods of supercontinent assembly coincide with high pCO₂ values due to the development of arid conditions which weaken the silicate weathering efficiency. These fundamental feedbacks between climate, carbon cycle and tectonic have been discovered by pioneer modelling studies and opened new views in the understanding of the history of Earth's climate. Today, some of the key features of the Phanerozoic climate can be explained by: (1) continental drift; (2) small continental blocks moving to tropical belts; and (3) modulation of the climate sensitivity to CO₂ by palaeogeography changes. Those results emphasize the need for a careful process-based modelling of the water cycle and climate response to the continental drift.

© 2013 Published by Elsevier B.V.

## Contents

<table><tr><td>1. Introduction</td><td>122</td></tr><tr><td>2. Pioneer modelling studies of the impact of palaeogeography on the global carbon cycle</td><td>123</td></tr><tr><td>3. Exploring the role of palaeogeography on the water and carbon cycle</td><td>123</td></tr><tr><td>4. Coupling climate and carbon cycle</td><td>124</td></tr><tr><td>4.1. GEOCARB and similar models</td><td>124</td></tr><tr><td>4.2. Sensitivity of the GEOCARB family of models to palaeogeography</td><td>125</td></tr><tr><td>4.3. The GEOCLIM model and the palaeogeography</td><td>125</td></tr><tr><td>4.3.1. Previous investigations with GEOCLIM</td><td>125</td></tr><tr><td>4.3.2. Exploring the impact of palaeogeography on the Phanerozoic CO₂ level</td><td>128</td></tr><tr><td>4.3.3. Role of small continents</td><td>131</td></tr><tr><td>4.3.4. Climate sensitivity and palaeogeographical configuration</td><td>132</td></tr><tr><td>4.3.5. GEOCLIM output and proxy data</td><td>133</td></tr><tr><td>5. Conclusions</td><td>134</td></tr><tr><td>Acknowledgements</td><td>136</td></tr><tr><td>Appendix A</td><td>136</td></tr><tr><td>References</td><td>137</td></tr></table>

## 1. Introduction

Palaeogeography and Earth climate are tightly linked. As noticed by Hay et al. (1990), the existence of a link between the latitudinal

<footnote>a Géosciences-Environnement Toulouse, CNRS-Observatoire Midi-Pyrénées, Toulouse, France</footnote>

<footnote>b Laboratoire des Sciences du Climat et de l'Environnement, CNRS-CEA, Gif-sur-Yvette, France</footnote>

<footnote>c Institut de Physique du Globe de Paris, Paris, France</footnote>

<footnote>* Corresponding author. Tel.: +33 5 61 33 26 15; fax: +33 5 61 33 25 60. E-mail address: yves.godderis@getobs-mip.fr (Y. Goddéris).</footnote>
distribution of land mass and the continental climate has been first postulated by Lyell (1830), almost a century before the publication of the continental drift theory (Wegener, 1912). Assuming a constant land surface, Lyell (1830) proposes that filling the Arctic region with emerged continents at the expense of equatorial territories would result in a severe polar glaciation and in a drop of the global temperature. Then comes the continental drift theory (Wegener, 1912) and the prediction of Köppen and Wegener (1924) that if all continents are located in the equatorial area for a given past period, continental and shelfal markers of climate will all indicate warm conditions, making a strong case for a link between palaeogeography and climate.

Since the end of the seventies, an abundant literature deals with the link between continental configuration and climate. Land mass distribution modifies the regional and global climate (Gyllenhaal et al., 1991), through changes in the oceanic circulation (Kennett, 1977; von der Heydt and Dijkstra, 2006; Lagabrielle et al., 2009 and reference therein; Lefebvre et al., 2012; Zhang et al., 2012), albedo feedbacks (Hay et al., 1990; Horton et al., 2010), and direct impact on the precipitation/evaporation balance (Barron et al., 1989; Otto-Bliesner, 1995). In addition to geological features, climate modelling studies explore the role of the continental configuration on pCO₂ thresholds needed to explain ice/green-house successions (Barron and Washington, 1982; Crowley and Baum, 1991, 1995; Barron et al., 1993; Gibbs et al., 1997; Hyde et al., 1999; Herrmann et al., 2004; Donnadieu et al., 2006; Craggs et al., 2012; Horton et al., 2012; Spicer et al., 2008). However all these studies share a common characteristic: they prescribe atmospheric greenhouse gas levels. Specifically, CO₂ level is a boundary condition and is thus assumed to be an external forcing of the climate system, whatever the continental configuration.

But things are more complex. Since weathering of Ca–Mg silicates of crustal rocks consumes atmospheric CO₂ at a rate depending on the Earth's climate, CO₂ levels are necessarily depending on the palaeogeography at long-term scale. Basically, the rock mineral dis-solution directly depends: (1) on the availability of water, thus on the continental area located within humid climatic belts; and (2) on the temperature, thus on the latitudinal distribution of the continents. Surprisingly, the interplay between palaeogeography and global carbon cycle has not been often studied in the recent years, despite the fact that it may be potentially a first order forcing of climate at the geological timescale. Instead of working on the geography effect, most of the efforts have been concentrated on the role of land plants (Berner, 2004), and of mountain building and physical erosion (Raymo et al., 1988; Raymo, 1991).

In this contribution, we propose an overview of the studies devoted to the impact of palaeogeography on the long-term carbon cycle (>>10⁵ yrs). We show that continental precipitations and runoff are highly sensitive to the palaeogeography. We demonstrate that the con- tinental configuration is a first order factor controlling the atmospheric CO₂ level and the Earth's climate evolution.

## 2. Pioneer modelling studies of the impact of palaeogeography on the global carbon cycle

According to the Earth thermostat theory (Walker et al., 1981), for a given solid Earth CO₂ degassing, the total CO₂ consumption by continen- tal weathering should be the same, whatever the location of continents. This is a direct consequence of the short response time of carbon and alkalinity in the ocean–atmosphere system, respectively 200 and 3 kyr (François and Goddéris, 1998) and of the negative feedback exerted by climate on weathering (Walker et al., 1981). The CO₂ consumption by silicate weathering is always tracking the CO₂ supply by solid Earth degassing, to avoid too large CO₂ fluctuations and repeated climatic catas- trophe (Berner and Caldeira, 1997).

Marshall et al. (1988) first proposed that the palaeogeography might be a key parameter of the evolution of the geological carbon cycle (Table 1). They found that atmospheric CO₂ must be at high levels

when most of the continents are located around the poles. Conversely, CO₂ must be low and global climate cold when continents are located around the equator. They argue that this is a direct consequence of a promoted CO₂ consumption by continental silicate rock dissolution when continents are located in the warm equatorial belt. Because temperatures are higher at the equator than at the pole, Marshall et al. (1988) note that the weatherability (i.e., the susceptibility to weathering) is higher when continents are close to the equator. As a consequence, the same total weathering flux will be reached at a much lower CO₂ level compared to a world where continents are located at the pole. Although their study was conceptual with generic continental configurations, Marshall et al. (1988) had already proposed that the late Precambrian glaciations might have been promoted by an equatorial location of the continents.

Worsley and Kidder (1991) further explore the Marshall et al. (1988) hypothesis through additional qualitative considerations. Similarly to Marshall et al., they suggest that a ring world (all continents located along the equator) should be characterized by lower atmospheric CO₂ (and hence colder global conditions) than a cap world (all continents around the poles). Provocatively, Worsley and Kidder (1991) even pro- pose that the averaged latitude of the continents alone specify CO₂ levels.

Marshall and collaborators were the first to couple a process-based climate model with a mathematical description of the continental weathering. The climate model used was quite simple: an energy-balanced model (EBM). This model calculates the Earth surface temper- ature as a function of latitude accounting for the vertical budget of energy and for the meridional heat transport. Similar models, coupling a description of the carbon cycle with a zonal energy-balanced model, have been used to explore the Phanerozoic CO₂ history (François and Walker, 1992; Goddéris and François, 1995; Goddéris and Joachimski, 2004; Tajika, 2007). However, it is important to note that if EBM models are reliable in terms of temperature predictions, they lack a physical description of the water cycle and, as such, are not able to calculate accurately the runoff, a key parameter of continental weathering (Dessert et al., 2003; Oliva et al., 2003). To overcome this limitation, a new generation of models has been developed.

## 3. Exploring the role of palaeogeography on the water and carbon cycle

Continental silicate weathering and associated CO₂ sink are heavily dependent on continental runoff, which is the difference between rainfall and evapotranspiration. As such, the relationship between palaeogeographical setting and hydrological cycle must be un- derstood, and this was beyond the abilities of the energy-balanced models previously used. The first study exploring the role of the continent distribution on the hydrologic cycle with an appropriate model was per- formed in 1989 by Barron and co-authors. They note that the water cycle is an essential part of the geochemical cycles (weathering, erosion, sedi- mentation), and that too much attention has been paid to temperature reconstructions at the expense of water cycle reconstructions. They explore the role of palaeogeography on the water cycle with the NCAR Community Climate Model (spatial resolution of 7.5° long × 4.5° lat). Contrary to energy-balanced models, the CCM includes a calculation of evaporation and condensation, and cloud cover. Such 3D climate models (also called General Circulation Models or GCM) explicitly simulate the water cycle, calculating the evaporation and rainfall on a 2-D map of the world, and the three-dimensional transport of water vapour in the atmosphere. Their main finding is that continental positions can strongly alter the hydrologic cycle (rainfall and runoff) if their positions affect the supply of moisture to the atmosphere. This is particularly the case in a simulation assuming a continuous belt of continents between 15°S and 15°N. Compared to the present day, this ring world is 8 °C globally warm- er, but global rainfall is maintained at its present day rate because the area of intense evaporation along the equator is lacking oceanic surfaces. Because of this intense continental evaporation, continental runoff
<table><caption>**Table 1**
A synthesis of numerical models exploring the role of palaeogeography on the geological carbon cycle.</caption><thead><tr><th>Model</th><th>Palaeogeography</th><th>Mechanistic calculation of the air temperature</th><th>Mechanistic calculation of the water cycle component</th><th>Spatial resolution</th><th>Parametric calculation of the air temperature</th><th>Parametric calculation of the continental runoff</th></tr></thead><tbody><tr><td>Marshall et al. (1988)</td><td>Conceptual</td><td>X</td><td>X</td><td>Zonal</td><td></td><td></td></tr><tr><td>Barron et al. (1989)</td><td>Conceptual and realistic</td><td>X</td><td>X</td><td>7.5°long × 4.5° lat (atmosphere) slab ocean</td><td></td><td></td></tr><tr><td>Worsley and Kidder (1991)</td><td>Conceptual</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Otto-Bliesner (1995)</td><td>Realistic</td><td>X</td><td>X</td><td>11.6°lat × 11.25°long (atmosphere) slab ocean</td><td></td><td></td></tr><tr><td>Gibbs et al. (1997)</td><td>Realistic</td><td>X</td><td>X</td><td>4.5°lat × 7.5°long (atmosphere) slab ocean</td><td></td><td></td></tr><tr><td>Berner (2004) and all GEOCARB versions</td><td></td><td></td><td></td><td>Global, no spatial resolution</td><td>X</td><td>X</td></tr><tr><td>Donnadieu et al. (2004)</td><td>Realistic</td><td>X</td><td>X</td><td>10°lat × 50° long (atmosphere) zonal dynamic ocean</td><td></td><td></td></tr><tr><td>Donnadieu et al. (2006, 2009), Goddéris et al. (2008a, 2008b), Nardin et al. (2011) GEOCLIM model</td><td>Realistic</td><td>X</td><td>X</td><td>4.5°lat × 7.5°long (atmosphere) slab ocean</td><td></td><td></td></tr></tbody></table>

decreases by a factor of 2 in this extreme ring configuration, with respect to present day runoff. Despite lower continental rainfall, a world with all the continents located beyond 50° latitude (polar location, cap world) displays a much higher runoff rate. These results contradict the finding of Worsley and Kidder (1991). If weathering is firstly dependent on runoff, then a cap world would be more weatherable than a ring world, thus characterized by a lower CO₂ level.

Otto-Bliesner (1995) further investigates the palaeogeographical forcing of runoff by running the NCAR 3D atmospheric general circulation model (i.e., AGCM) of the climate for 14 time slices spanning the Phanerozoic. The palaeogeographical reconstructions are taken from Scotese (www.scotese.com). Aiming at isolating the palaeogeographical effect on the global climate with a focus on the water cycle, Otto-Bliesner (1995) has fixed the CO₂ level at 280 ppmv, the solar constant at its present day value, and has assumed flat continents (no mountains) covered by rocky deserts. She finds that the relationship between runoff and palaeogeography is critical but more complex than the Marshall et al. (1988) and Worsley and Kidder (1991) simple scheme. The size, and not only the location of continents, must be considered. As long as continents are small and widely dispersed, precipitations can reach easily the inland areas. As a consequence, rainfall is not a limiting factor for runoff. Rather evaporation becomes the limiting factor, especially if the continents are located in the tropics, as already noticed by Barron et al. (1989). In this specific configuration (small continents and low latitude averaged location), runoff stays low, being constrained by the intense evaporation and responds weakly to global warming induced by a given CO₂ rise. Conversely, if small continents are moving poleward, the associated decreasing temperature induces a decrease in evaporation, possibly promoting runoff, and hence CO₂ consumption by weathering. Once cratons are aggregated into large continents (such as the Pangea), inland rainfall is strongly reduced and becomes limiting for runoff, instead of evaporation. Runoff is kept at low levels but its response to global warming events becomes stronger, forced by the rise in rainfall. Otto-Bliesner (1995) evidences the complex interplay between latitudinal location of the continents and their size. She emphasizes that the sensitivity of runoff to global warming might be very small for particular continental configuration, namely for small continents located at low latitudes.

In 1991, Bluth and Kump paved the way to a spatially-resolved calculation of the past weathering rates by defining lithological maps for the geological past, with a resolution of 2°long × 2°lat. Using their Quaternary map, they first calculated the weathering rates and fluxes during the last glacial maximum (Gibbs and Kump, 1994). In this study, the authors emphasize the importance of accounting for the spatial distribution in runoff when reconstructing past weathering rates. Going then into the distant past, Gibbs et al. (1999) explicitly calculate

the spatial distribution of the weathering rates of the continents over the last 250 Myr, using a 3D general circulation model for the atmosphere (GENESIS v1.02). The solar constant is kept at present day value and two sets of simulation are performed: one assuming a constant CO₂ level of 340 ppmv, and one at high CO₂ taken from simulations of the GEOCARB geochemical model (Berner, 1991). The continental runoff and temperature maps generated by the climate model for each time slice are used to calculate weathering rate for each of the 5 main lithological types (granite, basalt, carbonate, shale and sandstone). The fundamental conclusion of their work is that the overarching effect of continental collision and accretion is increased continentality (first identified by Kutzbach and Gallimore in 1989 for Pangean climates), which inhibits runoff and hence CO₂ consumption by silicate weathering. Times of super continent assembly should correspond to periods of high pCO₂, keeping all other parameters constant, such as the solid Earth degassing rate. Overall, chemical erosion is limited by the supply of water, which depends heavily on the palaeogeographical setting (super continent versus dispersed continents).

## 4. Coupling climate and carbon cycle

### 4.1. GEOCARB and similar models

All the above models (Marshall et al., 1988; Worsley and Kidder, 1991; Otto-Bliesner, 1995; Gibbs et al., 1999) give insights into the role of palaeogeography on silicate weathering without computing the CO₂ levels. This has been achieved by another type of model, the most cited being the GEOCARB model (versions I, II, III, and GEOCARBSULF; Berner, 2004, 2006). The models of the GEOCARB family compute the Phanerozoic atmospheric CO₂ through the solving of a set of equations, the most important being the carbon and alkalinity balance at the million year timescale, the carbon isotopic balance, the expression of silicate weathering as a function of environmental parameters such as the vegetation cover, physical erosion, and climate. The GEOCARB models being zero-dimensional, they only calculate a mean global temperature and runoff, without resolving the spatial distribution of the continental climate. In details, mean global runoff and temperature are calculated through the solving of two equations.

$$T_{glob} = T_0 + \Gamma \log(pCO_2) - W_s \frac{t}{570} + geog \quad (1)$$

$$R_{glob} = R_0 [1 + Y(T_{glob} - T_0)] f_D \quad (2)$$
In Eq. (1), $T_0$ is the present day global air temperature, $\Gamma$ is a parameter fixing the sensitivity of temperature to CO$_2$ variations, pCO$_2$ is the atmospheric CO$_2$ partial pressure in PAL (Pre-industrial Atmospheric Level), $W_s$ is the rate of temperature rise with time related to the increase in the solar constant. In this equation, two parameters are functions of the continental configuration. Firstly the factor *geog* is the correction applied to the calculation of $T_{glob}$ in response to continental drift. This factor is taken from Otto-Bliesner (1995) and is equal to the mean air continental temperature anomaly (relative to the present day control run) for each geological time slice from the simulations described above (pCO$_2$ = 1 PAL, no mountains, present day solar constant). Second, $\Gamma$ is the sensitivity of the air temperature to a CO$_2$ change for a given palaeogeography. Two values are adopted for the Phanerozoic: 3.3 for a priori warm climate modes (before 340 Ma, and between 260 and 40 Ma, corresponding to a global warming of 2.3 °C for a CO$_2$ doubling), and 4 for other periods assumed to be cold (2.8 °C warming for a CO$_2$ doubling). Those values come from GCM simulation at various CO$_2$ levels for the present day (assumed to be representative of the Phanerozoic cold modes) and the mid-Cretaceous (a model for warm modes) configurations (Hay et al., 1999; Kothavala et al., 1999). Recently, Park and Royer (2011) have revisited the value of the $\Gamma$ parameter, increasing it strongly during the cold phase of the Phanerozoic ($\Gamma$=8.7 to 11.5, corresponding to a 6 to 8 °C warming for a CO$_2$ doubling). Whatever the absolute values, it is important to note that the values of *geog* and $\Gamma$ are fixed prior to the GEOCARB simulations, and do not evolve as a function of the calculated atmospheric CO$_2$. They are thus forcing functions and may introduce some inconsistency in the model: for instance low calculated CO$_2$ may be coeval with low climate sensitivity $\Gamma$.

In Eq. (2), $R_0$ is the present day runoff, Y is a factor fixing the sensitivity of runoff to air temperature, and $f_D$ is an additional correction applied to runoff in response to continental drift. Y is fixed to 0.025 during warm climate modes (before 340 Ma and between 260 and 40 Ma, corresponding to a 2.5% increase in global runoff for a global mean temperature rise of 1 °C, Kothavala et al., 1999), and to 0.045 for other periods corresponding to a 4.5% increase in global runoff for a global mean temperature rise of 1 °C (Hay et al., 1999), a value close to the 4% estimated for the present day world by Labat et al. (2004) from river runoff time series. $f_D$ is taken from the Otto-Bliesner (1995) simulations and is equal to the runoff calculated for each timeslice assuming 1 PAL of CO$_2$, flat continents and a present-day solar constant divided by the present day runoff.

These two equations are critical. They define how these key parameters of the silicate weathering rate (temperature and runoff) evolve through Phanerozoic times in response to CO$_2$ and tectonic changes. Being zero-dimensional models, GEOCARB and GEOCARBSULF cannot explicitly account for the role of palaeogeography, despite its potential key role (Marshall et al., 1988; Worsley and Kidder, 1991; Otto-Bliesner, 1995). To overcome this critical limitation, the GEOCARB family of models includes several parameters derived from climate evolution studies (such as Otto-Bliesner, 1995), in such a way that the role palaeogeography is somehow included ($\Gamma$, *geog*, $f_D$, and Y to some extent). The benefit of the use of these parameters to calculate temperature and runoff instead of physically resolving the role of palaeogeography is the gain in time-computation. But the risk of introducing inconsistencies, or overestimating or underestimating the role of palaeogeography cannot be neglected.

## 4.2. Sensitivity of the GEOCARB family of models to palaeogeography

Most of the other published models calculating past CO$_2$ levels use the same formalism for the impact of palaeogeography on climate (Wallmann, 2001, 2004) or simply neglect the potential effect of continental drift (Bergman et al., 2004).

What is the contribution of palaeogeography to the Phanerozoic CO$_2$ signal calculated by the GEOCARB type of models? To answer this question, we built a simplified version of the GEOCARB model

(see Appendix A). When running the reference simulation for 21 time slices spanning the Phanerozoic, our simplified model version reproduces the general features of the GEOCARB family of models (Fig. 1): high CO$_2$ levels in the early Paleozoic, then a rapid decrease starting at 400 Ma. Low levels are predicted between 350 and 260 Ma. Then CO$_2$ rises again up to 10 times the present day values around 100 Ma. Finally, the late Cretaceous and Cenozoic are characterized by an overall CO$_2$ decrease towards the present day low level. We refer to this Phanerozoic trend as the “double-peak” feature in the following. Slight differences can be observed between this CO$_2$ curve and the latest GEOCARBSULF output, particularly higher levels around 300 and 100 Ma (Fig. 1). Indeed, we do not account for the role of the sulfur cycle, and silicate weathering laws and apparent activation energies for silicate dissolution are slightly different. Moreover, the temporal evolution of the continental area is taken from reconstructions by Blakey for the Paleozoic, Besse and Fluteau and Sewall for the Mesozoic and Cenozoic and differs thus slightly from the forcing used in GEOCARB. Finally, the time resolution of our simulations is much lower than Berner (2006). Anyway, our objective is not to reproduce precisely the GEOCARBSULF output, but rather to explore the contribution of the palaeogeographical setting to the “double-peak” evolution of atmospheric CO$_2$.

Then, we run our version of GEOCARB by fixing successively all parameters to their present day value except the parameters depending on the palaeogeography. At the end of this procedure (GLOBAL#4 simulation), only *geog*, $f_D$, $\Gamma$, and Y are allowed to fluctuate over the Phanerozoic. Solid Earth degassing rate is held constant and the organic carbon subcycle is assumed at steady state, which means that it does not impact the CO$_2$ level. A linear rise of the solar constant with time is still assumed. By doing so, we find that the double peak feature is mainly due to: (1) the impact of land plant evolution on continental weathering; (2) the variable volcanism; and (3) the fluctuations in the organic carbon cycle. When accounting only for the palaeogeographical parameters and the solar constant evolution, the major features of the GEOCARB CO$_2$ curve disappeared. The CO$_2$ displays now a general declining trend from about 6 PAL in the Cambrian to 1 PAL at present day, mainly forced by the rise in the solar constant. A sharp decrease around 460 Ma is predicted. This decline is related to the prediction of very high continental temperatures at 458 Ma in the simulations of Otto-Bliesner (1995), increasing sharply the *geog* factor, promoting silicate weathering, and thus forcing the predicted CO$_2$ to decrease to maintain the balance between silicate weathering and volcanic degassing. This warm climate is the result of the mean low latitude of the continental blocks during this period (Otto-Bliesner, 1995), and has been confirmed by more recent simulations (Nardin et al., 2011).

Corresponding calculated air temperatures are almost constant over the whole Phanerozoic, close to 14 °C (Fig. 2, GLOBAL#4 simulation), while runoff is predicted to be always larger than at present day, except around 250 Ma where it is very close to the present day value (Fig. 3, GLOBAL#4 simulation).

Overall, the sensitivity of the GEOCARB family of models to palaeogeography is low. The double-peak feature, which is the most visible feature of the CO$_2$ curve, is thus weakly dependent on the palaeogeographical setting.

## 4.3. The GEOCLIM model and the palaeogeography

### 4.3.1. Previous investigations with GEOCLIM

A new step in the carbon-climate modelling is achieved with Donnadieu et al. (2004) who introduce a new method to calculate the atmospheric CO$_2$ in the specific context of a supercontinent breakup at the end of the Proterozoic. The model they used, the first seed of the GEOCLIM model, couples a global carbon-alkalinity cycle to a climate model of intermediate complexity with a coarse spatial resolution (50°long × 10°lat). Continental weathering rates are resolved spatially, allowing to account explicitly for the continental configuration. Donnadieu et al. (2004) show that the Rodinia super continent,
<figure><img src="image_84.png" /><figcaption>**Fig. 1.** Evolution of atmospheric CO₂ calculated by zero-dimensional models of the GEOCARB type. The curve labelled GEOCARBSULF is taken from Berner (2006). GLOBAL is the simplified version of the GEOCARB model built for the present purpose. The GLOBAL label reminds that this model includes a zero-dimensional description of the global climate (only global air temperature and runoff are calculated). GLOBAL#1 is the same simulation where the solid Earth degassing is held constant. GLOBAL#2 is the same as GLOBAL#1, but with the organic carbon subcycle always assumed at steady-state. GLOBAL#3 further assumes a constant impact of land plants on weathering (same as today). Finally GLOBAL#4 is the same than GLOBAL#3 but lithology and erosion rates are kept at their present day values. The only forcing functions of GLOBAL#4 are thus the rising solar constant and the palaeogeographical changes.</figcaption></figure>

although located close to the equator, was characterized by low runoff, low weathering rates and hence high CO₂ (1800 ppmv). This is the result of the intense continentality, a finding in agreement with Gibbs et al. (1999). The subsequent breakup of the super continent into small pieces along the equator promotes the weatherability. The steady-state CO₂ level required for weathering to balance the solid Earth degassing equals about 500 ppmv in the dispersed continental configuration, corresponding to the onset of much colder climatic conditions compared to the super continent configuration (8 °C global cooling). When the weathering of freshly erupted basalts is accounted for, the breakup of the Rodinia might even be the ultimate cause for

the onset of the Neoproterozoic snowball glaciations (Donnadieu et al., 2004).

The Donnadieu et al. (2004) model was a first 'low-resolution' synthesis of the ideas developed by Marshall et al. (1988) and Gibbs et al. (1999). Evolving from this first model, GEOCLIM is a numerical model of the next generation able to calculate the past geological CO₂ levels (Fig. 4). It couples an 11-box model of the carbon-alkalinity-phosphorus-oxygen cycles (Goddéris and Joachimski, 2004) to the 3-dimensional climate model FOAM (Jacob, 1997). At each time step of the calculation and for any CO₂ level, the output of the FOAM general circulation model is used to estimate the coeval air temperature and

<figure><img src="image_85.png" /><figcaption>**Fig. 2.** Globally averaged surface temperature calculated in the various GEOCARB-type simulations described in Fig. 1. The curve GLOBAL#4 stands for the temperature change induced by the rising solar constant and the palaeogeographical changes solely.</figcaption></figure>
<figure><img src="image_86.png" /><figcaption>**Fig. 3.** Globally averaged continental runoff (normalized to present day) calculated in the various GEOCARB-type simulations described in Fig. 1. The curve GLOBAL#4 stands for the runoff change induced by the rising solar constant and the palaeogeographical changes solely.</figcaption></figure>

continental runoff with a spatial resolution of 7.5°long × 4.5°lat. While GEOCARB is a zero-dimensional model (only able to calculate global averages) in which corrective factors have been introduced from several offline GCM simulations (Otto-Bliesner, 1995), GEOCLIM explicitly resolve the impact of palaeogeography coevally with the calculation of the atmospheric CO₂ level. It is important to note that the two methods (GEOCARB vs GEOCLIM) are fundamentally different. In GEOCLIM, temperature, runoff, and more generally the hydrologic cycle, are resolved physically by the GCM (conservation of mass and energy in the atmospheric system) in GEOCLIM while GEOCARB relies on parametric laws. Furthermore, in GEOCARB, Y and geog are fixed assuming a priori

the timing of warm and cold climatic modes. In GEOCLIM, the sensitivity of runoff to temperature and of temperature to CO₂ is a function of the continental configuration, and fluctuates over the Phanerozoic. Furthermore, GEOCLIM allows the discrimination between continental and global temperatures.

In the following discussion, the continental runoff and its dependency on the continental configuration will be a key parameter. The mean global runoff calculated by the FOAM GCM for the present day configuration equals 330 mm/yr, in close agreement with the observations (313 mm/yr, UNH/GRDC database). Furthermore, the FOAM GCM is able to reproduce the observed spatial distribution of the continental

<figure><img src="image_87.png" /><figcaption>**Fig. 4.** Schematic description of the GEOCLIM model. A box model describing the geological cycles of carbon and alkalinity is coupled to the 3D FOAM climate model. Temperature and runoff calculated by the climate model are used to estimate the continental weathering fluxes, including the CO₂ consumption by silicate weathering at each time step with a spatial resolution of 7.5° long × 4.5° lat. For each continental configuration spanning the whole Phanerozoic, the GEOCLIM model is run until a steady-state is reached. Each continental configuration can thus be associated to a calculated CO₂ level, depending on other parameters such as the solid Earth degassing rate.</figcaption></figure>
runoff, the main discrepancy between the model and the data being a slight underestimation of the runoff for the northern part of South America by the model (Fig. 5). Anyway, in terms of runoff predictions, the GEOCLIM model appears to be quite robust.

Several studies have been undergone with the GEOCLIM model, keeping all parameters constant but the continental configuration and the slowly rising solar constant. Nardin et al. (2011) demonstrate that the drifting of the Baltica, Laurentia, and Siberia across the equatorial belt promotes continental weathering. As a consequence, CO₂ goes down since the Middle Ordovician from about 5600 ppmv to 2800 ppmv to maintain the balance between silicate weathering and the solid Earth degassing. About 70% of this CO₂ drop is fully related to the drifting of continents, confirming the key role played by palaeogeography. The remaining 30% are due to enhanced young volcanic rock weathering (Nardin et al., 2011). The lowest CO₂ level is reached at the end of the Ordovician/beginning of the Silurian. Then CO₂ rises again in excess of 6000 ppmv in the Early Devonian because of the accretion of the Baltica and Laurentia. The continental configuration thus promotes the onset of colder climatic conditions from the Middle Ordovician to the Early Silurian, in agreement with many geological indicators (Nardin et al., 2011).

Still using the GEOCLIM model, Donnadieu et al. (2006) and Goddéris et al. (2008a) have shown that a <10° northward drift of the Pangea supercontinent from the Carnian to the Rhetian induces a drastic decrease in atmospheric CO₂ from about 3000 ppmv down to 800 ppmv in the late Triassic. This cooling was related to the specific shape of Pangea and does not involve any continental break up process. The continental surface located in the tropical humid zone is almost doubled when Pangea drifts 10° north. Silicate weathering is promoted and CO₂ goes down to maintain the carbon balance. This is in agreement with sedimentary δ¹⁸O data which predicts a global cooling of 3.7 °C,

closely matching the 3.3 °C calculated by GEOCLIM. When accounting for the continental drift, what matters is not only the latitudinal location of the continents, but also their size and our ability to reconstruct accurately the past water cycle.

An improvement of the GEOCLIM model is the coupling with a global dynamic vegetation model. Vegetation, which responds to the continental drift, can change the albedo of the continental surfaces, the cycling of water through enhanced evapotranspiration, and the roughness of the continents. All these factors may change the temperature and precipitation patterns, and hence the weathering rates and the calculated CO₂ level. Le Hir et al. (2011) and Donnadieu et al. (2009) have both demonstrated that the key climatic effect of land plants with land plants is to modify the surface albedo. For instance, in the Devonian, the colonization of continents by plants decreases the albedo, resulting in a global warming trend. This global warming promotes CO₂ consumption by continental weathering, and CO₂ is forced to decline to maintain the carbon cycle balance (Le Hir et al., 2011).

### 4.3.2. Exploring the impact of palaeogeography on the Phanerozoic CO₂ level

GEOCLIM is run for the 21 time slices covering the Phanerozoic. The model is run until steady-state is reached, meaning that continental silicate weathering matches the prescribed solid Earth degassing. We calculate the corresponding atmospheric CO₂ and climate for each of the 21 time slices. The simulations account for dynamic vegetation by coupling the FOAM GCM to the LPJ biospheric model (Sitch et al., 2003; Donnadieu et al., 2009). The FOAM GCM is used in mixed-layer mode, i.e., the atmospheric model is linked to a 50-metre mixed-layer ocean, which parameterizes heat transport through diffusion. Ocean surface temperatures are thus calculated, but not the oceanic circulation. A

<figure><img src="image_88.png" /><figcaption>**Fig. 5.** (A) Present day observed distribution of the mean annual continental runoff (mm/yr) (UNH/GRDC: Composite runoff fields V1, 0.5°lat × 0.5°long). (B) Present day distribution of the mean annual runoff (mm/yr) calculated by the FOAM GCM under 280 ppmv CO₂.</figcaption></figure>
first set of simulations is run assuming the palaeogeography and the increase in the solar constant as the unique forcing functions (GEOCLIMtec simulations). The second set of simulations accounts additionally for the role of land plants on weathering, for a variable degassing, and for imbalances in the organic carbon subcycles, all these forcing functions coming from the GEOCARB model and added in GEOCLIM (GEOCLIM simulations).

Accounting only for palaeogeography and the linear rise in the solar constant, the GEOCLIM CO₂ curve (GEOCLIMtec, Fig. 6) differs considerably from the CO₂ history produced by the GEOCARB family of models. A huge drop in CO₂ is predicted during the Ordovician and an ample CO₂ peak is calculated between 280–220 Ma. CO₂ levels are much higher before the mid-Cretaceous. Very low CO₂ levels dominate after 130 Ma (Fig. 6). The palaeogeographical imprint on the CO₂ history is major, contrasting with the GEOCARB-type simulation in which the palaeogeography has only a second order impact on the CO₂ level. High CO₂ levels at the beginning of the Phanerozoic are forced by the polar location of most of the continents. As identified by Nardin et al. (2011), atmospheric CO₂ then decreases due to the northward drift of Laurentia, Baltica and Siberia. This period of change is followed by a period of relatively constant CO₂ levels, suggesting that the palaeogeographical configuration has played a little role in the onset of the Permo-Carboniferous glacial episode. The end Permian to middle Triassic period is characterized by very high CO₂ levels, promoted by the assembly of Pangea resulting in a dry continental climate. CO₂ abruptly declines in the latest Triassic in response to a northward drift of Pangea (Donnadieu et al., 2006; Goddéris et al., 2008b). A striking feature is the persistent low CO₂ levels predicted for the Jurassic and Cretaceous. If palaeogeography is the only forcing factor, the post-Triassic dispersed continental configuration should favour the onset of a globally cold climatic mode since the Jurassic. The end-Cretaceous and beginning of the Cenozoic are even marked by CO₂ levels below the present day value.

Part of the GEOCLIM CO₂ signal is forced by the response of runoff to palaeogeography. The amplitude of change in runoff over the Phanerozoic is larger in GEOCLIM (both for GEOCLIM and GEOCLIMtec

simulations) and in the simulations performed by Gibbs et al. (1999) compared to the simulations of Otto-Bliesner (1995) (Fig. 7). Furthermore, in Otto-Bliesner (1995), the runoff is generally higher than at present day, except for the Permo-Triassic period where it is close to the present day value. This explains partly why the CO₂ levels calculated by GEOCARB (GLOBAL#4 run, only solar constant and palaeogeography changing) are much lower than in the GEOCLIM run (GEOCLIM run) (Fig. 8). Consistently with Gibbs et al. (1999), the GEOCLIM runoff is conversely generally lower than at present day, at least before the mid-Cretaceous. Consequently, the GEOCLIM CO₂ is higher, because weathering is partly inhibited. The dry peak of the Permo-Triassic is also much more ample in more recent GCM simulations than in Otto-Bliesner (1995), explaining the existence of a calculated CO₂ peak (global runoff is below 60% of the present day value).

It is not easy to identify a unique cause for the calculated CO₂ features by GEOCLIM. Considering the GEOCLIMtec simulation, where continental drift is the sole forcing function, a plot of the calculated atmospheric CO₂ versus the total continental surface allows us to identify at least three groups of points (Fig. 9A). The first group has a standard behaviour (“standard” pool): atmospheric CO₂ decreases gently with the continental surface. More continental surface implies more silicate weathering, and hence less CO₂. Then the second group displays a range of high CO₂ values despite a roughly constant continental surface. This is the case for the 260–220 Ma time window and for the early Phanerozoic. These points define the “warm” Phanerozoic pool. Finally, a last group displays very low CO₂ levels whatever the total continental size. This group is identified here as the “cold” pool and covers the end of the Phanerozoic, starting from the middle Cretaceous.

The existence of a “warm” pool is linked to the aridity developing during the assembly of the Pangea super continent (Fig. 9B), but not only. This “warm” pool is also the consequence of a weak fraction of emerged continents in the warm and humid equatorial zone (10°S–10°N), where modelled weathering fluxes are maximized (Fig. 9C). This is the case during the Pangea phase and also during the early Phanerozoic. Consequently, the global weathering rates are inhibited and CO₂ rises. The existence

<figure><img src="image_89.png" /><figcaption>**Fig. 6.** Calculated evolution of the atmospheric CO₂ level over the Phanerozoic. The curve GLOBAL#4 stands for the CO₂ evolution induced by the rising solar constant and the palaeogeographical changes solely in a GEOCARB-type simulation. Similarly, the GEOCLIMtec curve is the output of the GEOCLIM model when the only forcings are the solar constant and the palaeogeographical changes. Accounting for the spatial distribution of the rainfall and temperature (and hence of weathering) in GEOCLIM impacts strongly on the calculated CO₂ history, suggesting that the palaeogeographical setting is a first order forcing function of the Phanerozoic climatic evolution. Red boxes figurate periods where the palaeogeography acts as a warming forcing function, promoting high CO₂ levels. Conversely, blue boxes stand for period where the continental configuration acts as a cooling forcing function. The ultimate causes for these trends are written in each box.</figcaption></figure>
<figure><img src="image_90.png" /><figcaption>**Fig. 7.** Time evolution of the global runoff calculated by 3D climate models. The Gibbs et al. (1999) ctrl simulation has been performed assuming a constant 340 ppm of CO₂ and moving continents. The Gibbs et al. (1999) high simulation was performed assuming variable CO₂ taken from GEOCARB. The Otto-Bliesner (1995) simulation was performed assuming 280 ppm of CO₂. Finally, GEOCLIM and GEOCLIMtec are respectively the output of the full GEOCLIM simulation, and of a GEOCLIM simulation where the only forcings are the rising solar constant and the palaeogeography.</figcaption></figure>

of a “cold” pool since the middle Cretaceous is linked to the ability of the Earth system to sustain an intense runoff over large continental surfaces (Fig. 9B). This is true even if the continental surface in the equatorial area stays generally below the mean Phanerozoic value. Runoff is no more dependent on the emerged continental size and stays at a high value. Since about 90 Ma, the continental configuration (highly dispersed continents, with only small polar continents) maintains high runoff and hence high CO₂ consumption by silicate weathering, forcing the Earth system into a cold state.

Consistently with Otto-Bliesner (1995), we found a dependence of runoff on the averaged latitudinal location of the continents (Pearson

correlation coefficient 0.60 in the GEOCLIMtec simulation) (Fig. 10A). The highest runoff rates coincide with the highest averaged latitudinal location. At first glance, this result is in agreement with previous calculations performed with generic continental configurations in which an equatorial “ring” world (extending from 17°S to 17°N) was characterized by low global runoff (Hay et al., 1990). Indeed, large continental surfaces located along the equator may inhibit the evaporation in the tropical convergence zone, decreasing the global runoff. However, in our simulations, the increase in global runoff with the mean latitude of the continents can be explained by simple “geometric” consideration. Whatever the considered period, runoff is always minimal between

<figure><img src="image_91.png" /><figcaption>**Fig. 8.** Phanerozoic CO₂ history calculated by the zero-dimensional GEOCARB-type model and by the three-dimensional GEOCLIM model. GLOBAL#4 and GEOCLIMtec are simulations where the only forcings are the rising solar constant and the palaeogeography.</figcaption></figure>
<figure><img src="image_92.png" /><figcaption>**Fig. 9.** (A) Atmospheric CO₂ versus the total emerged continental surface in the GEOCLIMtec simulation. The black numbers stand for the age of each dot, in million years. (B) Global continental runoff normalized to its present day value versus the total emerged continental surface in the GEOCLIMtec simulation. (C) Atmospheric CO₂ versus the percentage of the total continental surface located in the 10°S–10°N latitude band in the GEOCLIMtec simulation. The orange envelope represents the “warm” pool, the blue envelop the “cold” pool and the green envelop the “standard” pool. Each of these pools is exclusively forced by the palaeogeographical forcing.</figcaption></figure>

roughly 15° and 35° latitude, because this zone coincides with the sub-sidence of the Hadley cell. The width of this low rainfall/high evaporation area is somehow dependent on the continental configuration, but not critically. As a consequence, when the continental fraction increases in the 15°–35° zone, we might expect that the global runoff will decrease simply because the continental area exposed to very low rainfall increases. This “geometric” effect is illustrated on Fig. 10B, where the global runoff is negatively correlated to the continental fraction in the 15°–35° latitude band (Pearson correlation coefficient — 0.84). The positive correlation between the averaged latitudinal position of the continents and the runoff arises from the “fortuité” existence of a negative correlation between the averaged latitudinal location of the continents and the continental fraction between 15°–35° latitude (Fig. 10C).

It is worth mentioning that both model (GEOCLIM and GEOCARB) have a rather different behaviour with respect to the Permo-Triassic dry event. Although less pronounced than in GEOCLIMtec, the runoff calculated in GLOBAL#4 simulation also displays a decrease around 250 Ma (both GEOCLIM and GEOCARB-style simulations being performed with the palaeogeographical forcing only). Despite this runoff decrease by about 30% from the Carboniferous to the Permo-Triassic, calculated atmospheric CO₂ in GLOBAL#4 does not display any rise in CO₂ at the same time, while GEOCLIMtec peaks at almost 20 times the present day level with a 40% decrease in runoff (Fig. 6). The reason for these distinctive behaviours is rooted in the parametric expression of silicate weathering as a function of runoff. In GEOCLIMtec, silicate weathering is linearly dependent on runoff, as suggested by recent complications of field studies of silicate weathering (Dessert et al., 2003; Oliva et al., 2003). In GLOBAL#4, silicate weathering linearly increases with the runoff to the power 0.65, to account for dilution effects (Berner, 1994, 2004). This dilution effect is introduced to fit the

database of Dunne (1978) for Kenyan rivers and of Peters (1984) for rivers of the United States. As a result, the response of weathering to runoff changes is seriously damped in GLOBAL#4, explaining why there is no CO₂ peak during the Permo-Triassic low runoff interval. This illustrates that the superposition of numerous parametric laws of various origin to describe the climate/carbon system, may finally hide the role of the palaeogeography on the atmospheric CO₂ level. This result also stresses that runoff and its impact on the CO₂ consumption by rock weathering must be modelled accurately. The use of a GCM to simulate the water cycle in GEOCLIM is a major improvement compared to zero-dimensional models, but we are still dependent on the reliability of parametric laws to describe silicate weathering. Because these parametric laws are fitting functions of the field data (for instance, HCO₃⁻ fluxes at the watershed outlet versus mean annual runoff, Dessert et al., 2003), the most complete databases must be used. Both the Oliva et al. (2003) and Dessert et al. (2003) databases cover a wide variety of climatic environments. But the spreading of their data points, particularly important for granitic environments, still questions the reliability of parametric laws. This is an inherent limitation of the method.

### 4.3.3. Role of small continents

The drastic decrease in atmospheric CO₂ in the Rhetian stage of the Triassic is mainly linked to the Northward drift of Pangea (Goddéris et al., 2008a). But there is an additional contributor to this sharp cooling trend. South China is located in the humid tropical area at the same time. In the simulation where the continental configuration is the only forcing (GEOCLIMtec), silicate weathering in South China is accounting for 17% of the global CO₂ consumption by weathering (Fig. 11). If CO₂ consumption by silicate weathering in South China during the Rhetian is removed from the simulation, atmospheric CO₂ would rise from
<figure><img src="image_93.png" /><figcaption>**Fig. 10.** (A) Global runoff as a function of the averaged latitudinal positions of the continents. (B) Global runoff as a function of the continental fraction within the latitude band extending from 15° to 35° latitude, in both hemisphere. (C) Averaged continental latitude of the continents as a function the continental fraction within the latitude band extending from 15° to 35° latitude, in both hemisphere. Runoff output from the GEOCLIMtec simulation. The black numbers stand for the age of each dot, in million years.</figcaption></figure>

3.15 PAL (880 ppmv) to 5 PAL (1400 ppmv). Because of the logarithmic dependence of temperature on CO₂, the subsequent global warming of the Rhetian stage is significant (3 °C continental warming). This sensitivity test emphasized the potential key role played by small continents on the global carbon cycle, especially when they are located in the humid tropical belt (Oglesby and Park, 1989). This result somehow contradicts Otto-Bliesner (1995), who found that small continental blocks located in the tropical area might be characterized by low runoff, because evapotranspiration exceeds rainfall. If this was the case, small

continental blocks should have only a minor effect on the long-term climate evolution. Conversely, our simulations show that small tropical continents might contribute significantly to the global CO₂ sink by weathering.

### 4.3.4. Climate sensitivity and palaeogeographical configuration

In the recent years, geological climate reconstructions have been used to assess the climate sensitivity to atmospheric CO₂ changes (Berner and Kothavala, 2001; Lunt et al., 2010; Park and Royer, 2011).

<figure><img src="image_94.png" /><figcaption>**Fig. 11.** CO₂ consumption by silicate weathering (10¹⁰ mol CO₂/yr) calculated for the Rhetian stage of the Triassic in the GEOCLIMtec simulation (pCO₂ = 5040 ppmv). South China crosses the equatorial humid belt and consumes 17% of the total CO₂ consumption.</figcaption></figure>
We show here that climate sensitivity is substantially dependent on the palaeogeographical setting. Sensitivity is maximum during the Early Cambrian and during the Early to Middle Triassic period, exceeding 4 °C for a CO₂ doubling (Fig. 12). Except for this two particular geological period, the sensitivity is more less constant oscillating around 2.8 °C with a minimum at 2.4 °C at 90 Ma. What make the Early Cambrian and the Early to Middle Triassic so sensitive to a CO₂ doubling is not straightforward at first. It is important to note that these calculations of the climatic sensitivity account for fast and slow climatic feedbacks, excepting the role of the cryosphere, which is not modelled.

Climate sensitivity seems to respond partly to fluctuations in the continental surface as a function of latitude (Fig. 12). It roughly follows the evolution of the continental fraction between 60° and 30° latitude. Sensitivity increases when mid-latitudinal continental fraction rises, as it is the case in the Early Cambrian and during the Pangean period. But in details, things are trickier. The rapid increase in climate sensitivity from about 2.8 °C to more than 4 °C around 240 Ma is correlated to a strong decrease in the polar continental fraction. Also during the end Mesozoic and the Cenozoic, climate sensitivity decreases while the mid-latitude continental fraction stays roughly constant.

Runoff dependence on global temperature is also modulated by the palaeogeographical evolution. During most of the Paleozoic, from the Ordovician to the Mississippian, the runoff sensitivity to temperature is around 1.5 cm yr⁻¹ K⁻¹ (Fig. 13). Then it rapidly decreases with the assembly of Pangea, reaching 0.6 cm yr⁻¹ K⁻¹. The Mesozoic is then characterized by a slow increase in runoff sensibility. Finally a major rise is predicted in the Miocene, where sensitivity goes back to Paleozoic values.

### 4.3.5. GEOCLIM output and proxy data

The objective of this contribution is to emphasize the key role played by the palaeogeographical configuration on the carbon cycle and on the climate evolution. At this stage, it is not easy to compare the GEOCLIM outputs with proxy data, because GEOCLIM has been designed to model mechanistically the role of palaeogeography, and does not aim at being a model of the complete carbon cycle. The other processes, including the role of the organic carbon cycle, of the vegetation on weathering, are either not modelled (GEOCLIMtec simulation), or

modelled using simple parametric laws (GEOCLIM simulation). Nevertheless, some conclusions can be drawn from the comparison between proxies and model output.

Regarding the CO₂ level, the agreement between the GEOCLIM outputs and the CO₂ proxies is acceptable only for the last 200 Ma of the Earth's history, considering the large uncertainties inherent to the reconstructions (Fig. 14A). Between 200 and 400 Ma, GEOCLIM tends to overestimate the past CO₂ levels. This overestimation is particularly striking in the Triassic, where the supercontinental configuration forces GEOCLIM to very high CO₂ values. These high values are not recorded in the CO₂ proxies. Then, the 250–350 Ma period is characterized by very low proxy-based CO₂ values, around 500 ppmv, while GEOCLIM predicts values in excess of 1000 ppmv.

We also compare the mean annual sea surface temperature (SST) calculated by GEOCLIM within the 30°S–30°N latitudinal band to the temperature reconstructed using the tropical δ¹⁸O data compiled by Prokoph et al. (2008) (Fig. 14B). To translate δ¹⁸O into SST, we follow the method described in Veizer et al. (2000). δ¹⁸O data have been detrended, using a root mean square linear fit. The detrending is applied to avoid unrealistically high temperatures mainly before 100 Ma. Then a running average has been calculated (20 Myr time window), including the calculation of the range containing 68% of the data points. Finally, we set the seawater δ¹⁸O to −0.7 permil in the Cenomanian (a period where continental ice is assumed to be absent), and to 1 permil for the coldest stage of the Phanerozoic (most positive δ¹⁸O detrended value). Because this represents a significant handling of the data, we express the result in terms of relative changes to the present day instead of absolute temperatures. Furthermore, the detrending procedure may hide real long-term climatic trends in the oxygen record. This is particularly critical for the last 100 Ma of Earth history, where the δ¹⁸O signal provides a realistic temperature signal even without detrending (Donnadieu et al., 2009). For that reason, we also display on Fig. 14B the temperature calculated for the last 100 Ma without detrending of the oxygen isotopic data.

The spreading of the temperature values is quite large. The most striking feature is the existence of a warm peak in the Late-Permian/Early Triassic, in agreement with the GEOCLIM output (Goddéris et al., 2008a). A recent study has evidenced the

<figure><img src="image_95.png" /><figcaption>**Fig. 12.** Mean global temperature change for a CO₂ doubling as a function of age, calculated by the GEOCLIM model (solid line, scale on the right). The dotted blue line stands for the climate sensitivity prescribed in the GEOCARB model. Superimposed is the fraction of continents in various latitude bands (scale on the left).</figcaption></figure>
<figure><img src="image_96.png" /><figcaption>**Fig. 13.** Runoff change in cm/yr for an increase of 1 °C in the global temperature, calculated by the GEOCLIM model using a three-dimensional calculation of the water budget (solid line). The dotted blue line is the same calculated by the zero-dimensional GEOCARB model, which includes a parametric description of the water cycle.</figcaption></figure>

persistence of very high temperatures in the tropical ocean over the whole Early Triassic, in excess of 35 °C (from conodont apatite δ¹⁸O) (Sun et al., 2012). As noted by the authors themselves, high temperatures for the 5 Myr of the Early Triassic require strong and persistent greenhouse conditions. These high temperature conditions may have lasted until the end of the Carnian stage, as inferred from brachiopod δ¹⁸O (Korte et al., 2005). In the GEOCLIM simulations as well as in the Gibbs et al.'s (1999) study, these long lasting warm conditions can be sustained by the developing aridity, limiting the efficiency of continental weathering and promoting high CO₂ levels for millions of years. Conversely, the major temperature decrease occurring after 100 Ma and predicted by the model is not observed in the Prokoph et al. (2008) database when the oxygen isotopic data are detrended. But a cooling of about 5 °C in the tropical SST is reconstructed between 100 and 0 Ma when the oxygen isotopic data are not detrended, in close agreement with the model output (Donnadieu et al., 2009).

Several conclusions can be drawn. (1) When considering the early Triassic warm event, there is a clear disagreement between CO₂ proxy data and temperature data that we already note in a previous contribution (Goddéris et al., 2008a). A warm event seems to be associated to low CO₂ values. A coupled carbon-climate model like GEOCLIM, and even GEOCARB, cannot calculate at the same time high temperatures and low CO₂ levels. This apparent discrepancy between proxy-based temperatures and proxy-based CO₂ levels might be linked to the pretty large uncertainties inherent to the reconstruction methods, but this is beyond the scope of this study. (2) It appears that the palaeogeographical forcing is not responsible for the Permo-Carboniferous cool mode. (3) Reconstructing δ¹⁸O-based temperature curves for the whole Phanerozoic is not an easy task. However, such composite regional or latitudinal curves would greatly help at validating spatially resolved numerical models. The various available databases are not really self-consistent. They mix phosphate and carbonate signals, the mathematical relationship between δ¹⁸O and temperature fluctuates from one study to the other, and the data reflects the temperatures at various depths, depending on the organisms. Furthermore, the role of diagenesis is still not fully understood. Most of the geochemists now rely preferentially on phosphatic data for periods older than the Jurassic (Trotter et al., 2008;

Joachimski et al., 2009), because phosphatic materials are much less sensitive to diagenesis than pre-Jurassic carbonates. But the question of the phosphate-water fractionation equation is still debated (Longinelli, 2013; Pucéat et al., 2013) and the question of the seawater δ¹⁸O value is still critical. (4) Non-quantitative data should not be neglected when using spatially resolved models. Indeed, the GEOCLIM predictions of a very warm early Triassic are in agreement with sedimentological evidence for the deposition of redbeds and massive evaporites. The calculated sharp cooling during the latest Triassic coincides with major climatic changes impacting various regions of the world. Pedogenic and clay mineral data (smectite-rich redbeds) indicate a switch from arid/semi-arid conditions to year-round humid environments, a switch reproduced in our simulations (see Goddéris et al., 2008a, for references and more details). (5) The aim of the modelling efforts described in this work is not really to improve the results in terms of past atmospheric CO₂, but rather to improve the method by introducing physics in the numerical reconstructions of the past climate and carbon cycle. The GEOCARB model paved the way of the quantitative reconstructions of the Phanerozoic carbon cycle. The method used was mostly parametric. Now, recent studies relying on complex coupled climate-carbon numerical models have shown that the palaeogeographical factors have been underestimated. This review demonstrates that we now understand how and to what extent the palaeogeography really impacts the past carbon cycle and climate. This is a good starting point for revisiting the other processes at play, such as the links between weathering and land plants, between weathering and relief, among other parameters.

## 5. Conclusions

Two modelling pathways have been followed in the recent years for reconstructing the geological atmospheric CO₂ content. Zero-dimensional models (the so-called GEOCARB model family in this contribution) have been developed with the objective of accounting for as much processes as possible, such as the role of land plants or physical erosion on the weathering fluxes. Such models also include a simple parametric description of the climate response to the CO₂ content, the solar constant rise or the continental configuration on a global basis.
<figure><img src="image_97.png" /><figcaption>**Fig. 14.** (A) Comparison between the atmospheric CO₂ level calculated by the GEOCLIM model (dashed lines, the upper curve is the GEOCLIM simulation, the lower curve is the GEOCLIMtec simulation) and the CO₂ levels from various proxies (20 Myr moving averages). (B) The tropical sea surface temperature (averaged 30°S–30°N latitude) calculated by the GEOCLIM model (dashed lines, the upper curve is the GEOCLIM simulation, the lower curve is the GEOCLIMtec simulation, scale on the left). Also shown is the temperature anomaly for the same zonal area estimated from sedimentary carbonate δ¹⁸O (scale on the right). See text for more details.</figcaption></figure>

Those models gather information from other model runs or from field data and include those as parametric laws. Mixing several parametric descriptions from various origins is always at risk of losing the inherent physical behaviour of a complex system such as the Earth system. Nevertheless, these zero-dimensional models have the advantage of being integrative (many processes are considered), and fast in terms of computation time. They also generate past CO₂ reconstruction that can be compared to available data because they aim at accounting for all the factors controlling the evolution of the global carbon cycle.

Since the end of the 80s, physically-based models have also been developed to describe specifically the impact of the continental configuration of the global carbon cycle and climate geological evolution. Such models all include a spatially-resolved description of the climate system. This is the minimum requirement to consider physically the role of the palaeogeography on the global carbon cycle. Those models, one

to three-dimensional, proceed from a different philosophy than the zero-dimensional method. The objective is no more to calculate a time evolution of the atmospheric CO₂ than can be compared readily to data, but to explore the sensitivity of the Earth system to one parameter, here the palaeogeographical configuration. These studies come to four major findings:

(1) The palaeogeography is a key controlling factor of the geological evolution of the atmospheric CO₂. By modulating runoff and air temperature above the continental surface, it strongly controls the CO₂ consumption by silicate rock weathering. Supercontinent and polar configurations favour the onset of very high CO₂ pressure in the atmosphere. Conversely, since about the mid-Cretaceous, the palaeogeography promotes the onset of low CO₂ and cold climatic conditions. The role of the palaeogeography has been largely underestimated in zero-dimensional models.
(2) Some of the key features of the Phanerozoic climate can be explained by the modulation of the carbon cycle by continental drift. This is the case for the early Paleozoic (Cambrian to early Silurian) global cooling, the warm and arid early Triassic, and the cold trend of the Cenozoic.

(3) Small continental blocks may contribute significantly to the total CO₂ consumption by weathering when entering the sub-tropical warm and humid zone.

(4) Climate sensitivity is substantially dependent on the continental configuration. As such, though the past might be the key to estimate the amplitude of the future climatic changes, it is important to understand the processes linked to the palaeogeography that have modulated the Earth climatic sensitivity in the past.

Instead of using parametric relationships to estimate global runoff and temperature, the models described in this review (Otto-Bliesner, 1995; Gibbs et al., 1999; Donnadieu et al., 2004, 2006; Goddérис et al., 2008b; Nardin et al., 2011; Lefebvre et al., 2013) resolve spatially the climatic conditions using process-based climate models. Such models implement the energy budget of the Earth atmosphere and oceans and a physical description of the water cycle. Those models are not yet able to propose a CO₂ history that can be compared to data. Indeed, despite a major improvement in the physical basis of the climate compared to parametric models, one major constraint remains the mathematical expression for CO₂ consumption by continental weathering. In Gibbs et al. (1999) and in the GEOCLIM model, continental weathering is calculated on each grid cell assuming runoff (and temperature) dependencies. But weathering is a complex process. It also depends on physical erosion, vegetation cover, soil thickness, and lithology. Regarding the vegetation cover, a few studies try to explore physically the climatic (Donnadieu et al., 2009; Le Hir et al., 2011) and geochemical (Taylor et al., 2012) impacts of land plants within the framework of a three-dimensional climate and geochemical modelling. All these studies show the critical role played by the land plant distribution on the weathering fluxes, and hence on the global carbon cycle. Regarding the thickness of the soil cover, protecting the bedrock from weathering when it becomes too thick, Goddérис et al. (2008b) have shown, using the GEOCLIM model that this effect has modulated the CO₂ levels. But the role of large mountain ranges on the global carbon cycle is still to be included in those spatially-resolved numerical model. Finally, although Gibbs et al. (1999) suggest that spatial variations in lithology account only for small fluctuations in the total chemical erosion rates, this was not the conclusion of Taylor et al. (2011, 2012), Nardin et al. (2011) and Kent and Muttoni (2013). Indeed these authors found that the bulk of weathering may come from basaltic hotspots.

In summary, it appears that we are slowly moving from simple parametric representations of the geological Earth system towards more comprehensive and physically-based descriptions. What we have learnt in the recent years is that the overall picture we had of the carbon cycle and climate evolution may drastically change when more complex models are used. Further works are clearly needed to account for a better representation of the weathering rates for instance, but they should include, as far as possible, a description of the physics at play.

## Acknowledgements

Funding for this work has been provided by the CNRS-INSU Eclipse and Syster programs, and by the ACCRO-Earth (ANR-06-BLAN-0347), COLORS (ANR-09-JCJC-0105), TERRES (ANR-10-BLAN-0607) and Anox-Sea (ANR-12-BS06-0011) ANR projects. We thank Jeffrey Park and Lee Kump for their constructive and fruitful reviews of this contribution.

## Appendix A

Basically, the calculation of the Phanerozoic atmospheric CO₂ pressure by the GEOCARB family of models relies upon a set of equations. Although the complete GEOCARB model is more complex (Berner, 2004), we only show here the key equations. The objective is to emphasize the role played by the palaeogeography on the climate evolution of the Earth, and not to reconstruct strictly the GEOCARB model.

(1) The carbon balance: given the short residence time of carbon in the ocean/atmosphere system compared to the geological time-scale, the exospheric carbon cycle is assumed to evolve through a succession of steady states:

$$F_{vol} + F_{MOR} + F_{cw} + F_{ow} = F_{cd} + F_{od} \quad (A1)$$

where $F_{vol}$ is the supply of CO₂ by subaerial volcanism, $F_{MOR}$ is the degassing flux of CO₂ at mid-oceanic ridges, $F_{cw}$ is the carbon flux released by continental carbonate rock dissolution, and $F_{ow}$ is the carbon flux released to the river system by the oxidation of reduced sedimentary rocks. The sink side corresponds to the storage of carbon on the seafloor either by carbonate sedimentation ($F_{cd}$) or by organic carbon burial into the sediment ($F_{od}$).

(2) The alkalinity balance:

$$F_{swo} + F_{swy} + F_{cw} = F_{cd} \quad (A2)$$

$F_{swo}$ and $F_{swy}$ are respectively half of the alkalinity flux created during the weathering of continental old cratonic silicate rocks and of young volcanic rocks. These two terms can be written as follows:

$$F_{swo} = (1 - \sigma_y) F_{sw} \quad (A3)$$

$$F_{swy} = \sigma_y F_{sw} \quad (A4)$$

where $F_{sw}$ equals half of the total alkalinity flux created during the weathering of continental silicate rocks. $\sigma_y$ is the fraction of the silicate rock outcrop occupied by young volcanic rocks.

Altogether, these two Eqs. (A1) and (A2) can be written as follows, stating that all sources of exospheric carbon by solid Earth degassing is compensated by the continental silicate weathering and the net sink of organic carbon (organic carbon burial minus the oxidation of reduced sedimentary carbon exposed at the Earth surface).

$$F_{vol} + F_{MOR} = F_{swo} + F_{swy} + (F_{od} - F_{ow}) \quad (A5)$$

(3) The organic carbon disequilibrium ($F_{od} - F_{ow}$) can be constrained by considering the carbon isotopic balance.

$$\begin{aligned} \delta_{vol} F_{vol} + \delta_{MOR} F_{MOR} + \delta_{cw} F_{cw} + \delta_{ow} F_{ow} \\ = \delta_{oc} F_{cd} + (\delta_{oc} - \varepsilon_p) F_{od} \end{aligned} \quad (A6)$$

where $\delta_{vol}$ is the averaged volcanic gas $\delta^{13}$C signature, and $\delta_{MOR}$ is the mantle carbon isotopic composition. $\delta_{cw}$ and $\delta_{ow}$ are the averaged $\delta^{13}$C of continental carbonated and reduced sedimentary rocks. $\delta_{oc}$ is the measured seawater $\delta^{13}$C (Veizer et al., 1999) and $\varepsilon_p$ is the photosynthetic carbon isotopic fractionation (Hayes et al., 1999).

(4) The strontium isotopic budget can be used to constrain the relative importance of young volcanic rocks versus older silicate rocks (Berner, 2004). Given the residence time of strontium in seawater (1.5 to 2 Ma) and the time resolution