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

The Ordovician Period (485–443 Ma) is characterized by abundant evidence for continental-sized ice sheets. Modeling studies published so far require a sharp CO2 drawdown to initiate this glaciation. They mostly used nondynamic 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 CO2 and paleogeography. We conduct experiments for a wide range of CO2 (from 16 to 2 times the preindustrial atmospheric CO2 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-CO2 relationship is highly non-linear when ocean dynamics are taken into account. Two climatic modes are simulated as radiative forcing decreases. For high CO2 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 CO2 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 CO2 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 explanation to the sudden Late Ordovician Hirnantian glacial pulse that does not require any large CO2 drawdown.


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 Biodiversification 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.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Atmospheric CO 2 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 2 values.Climate proxies (Yapp and Poths, 1992) and continental weathering models (Berner, 2006;Nardin et al., 2011) indicate that the Ordovician CO 2 atmospheric partial pressure (pCO 2 ) was equal to some 8-20 times its preindustrial atmospheric level (1 PAL = 280 ppm).Fainter Sun (Gough, 1981) somewhat reduces this apparent discrepancy between high CO 2 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 2 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 2 levels are probably the best estimates (Vandenbroucke et al., 2010;Pancost et al., 2013).
While the role of the pCO 2 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., 2003Herrmann et al., , 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 2 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 2 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 2 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.

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 continentalsized ice sheets.(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 andFinnegan 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(Brenchley et al., , 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 paleogeographi-cal 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.3 Climate simulation results

Sensitivity to an atmospheric CO 2 decrease
Figure 2 presents the relationship between the global mean surface air temperature (SAT) and the pCO 2 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., 2003Herrmann et al., , 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 2 decrease and the temperature-CO 2 relationship is linear in a log-axis graph (Fig. 2a).The model sensitivity, 3 • C per halving of CO 2 , 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., 2003Herrmann et al., , 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 2 forcing.The temperature-CO 2 relationship is non-linear (Fig. 2a).Two -For high CO 2 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 2 ) and lies within the range 2.1-4.7 • C defined with various up-todate GCMs for present-day simulations (IPCC, 2013).
The temperature decrease associated to a CO 2 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  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 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 2 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 2 , 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 continental configuration.A similar evolution is observed.With the slab model (Fig. 2), the temperature-CO 2 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 major climatic instability within the interval 12-6 PAL.The tipping point showing the shift from the warm to the cold climatic 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 critical role in the shift to the cold climatic mode, (iii) reaching the cold climatic mode requires a lower CO 2 level at 450 Ma than at 470 Ma.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.

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 2 threshold, we compare the climatic features at both time slices for the same pCO 2 .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 2 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  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 ∼ × T × ρ × C p where is the strength of the overturning circulation, T the temperature difference between the northward and the southward limbs of the overturning circulation, ρ 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 winddriven 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 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 2 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 polar latitudes, but it increases the global climate sensitivity to the radiative forcing as soon as NH sea ice begins to expand.

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 winddriven structure of the OHT.A primary maximum at lowlatitudes 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 in 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 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 subpolar gyre is sufficient to stabilize a finite sea ice cap at midlatitudes, 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 2 stabilize in Eq2, whereas the experiments conducted with high CO 2 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 lowlatitude 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  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 MITgcm 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.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 midlatitudes.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 .

Figure 1 .
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).
Figure 2. (a) Mean annual, globally averaged surface air temperature for various pCO 2 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 −2 .

Figure 3 .
Figure3.SST (color shading) and permanent sea ice cover (bluewhite 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.

Figure 4 .
Figure 4. Mean annual, zonally averaged SAT simulated for each pCO 2 at 450 Ma with the fully coupled version of FOAM.

Figure 5 .
Figure 5. Zonally averaged sea ice cover fraction (expressed in percents) simulated for each pCO 2 at 450 Ma with the fully coupled version of FOAM.DJF: December, January and February.JJA: June, July and August.

Figure 6 .
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.

Figure 7 .
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 8 .
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) andScotese and McKerrow (1991, Figs. 4 and 5).

Figure 9 .
Figure 9. Difference in total, advective and diffusive OHT between the 450 Ma-12X and the 470 Ma-12X simulations.

Figure 10 .
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.
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 12 .
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 2 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.

www.clim-past.net/10/2053/2014/ Clim. Past, 10, 2053-2066, 2014Table 1 .
Overview of the coupled simulations.The name of each experiment includes the age of the paleogeography used in input followed by the pCO 2 expressed in PAL.These abbreviations are used throughout the manuscript.Slab runs are not shown.