Tracing the Snowball bifurcation of aquaplanets through time reveals a fundamental shift in critical-state dynamics

The instability with respect to global glaciation is a fundamental property of the climate system caused by the positive ice-albedo feedback. The atmospheric concentration of carbon dioxide (CO₂) at which this Snowball bifurcation occurs changes through Earth's history, most notably because of the slowly increasing solar luminosity. Quantifying this critical CO₂ concentration is not only interesting from a climate dynamics perspective but also constitutes an important prerequisite for understanding past Snowball Earth episodes, as well as the conditions for habitability on Earth and other planets. Earlier studies are limited to investigations with very simple climate models for Earth's entire history or studies of individual time slices carried out with a variety of more complex models and for different boundary conditions, making comparisons and the identification of secular changes difficult. Here, we use a coupled climate model of intermediate complexity to trace the Snowball bifurcation of an aquaplanet through Earth's history in one consistent model framework. We find that the critical CO₂ concentration decreased more or less logarithmically with increasing solar luminosity until about 1 billion years ago but dropped faster in more recent times. Furthermore, there was a fundamental shift in the dynamics of the critical state about 1.2 billion years ago (unrelated to the downturn in critical CO₂ values), driven by the interplay of wind-driven sea-ice dynamics and the surface energy balance: for critical states at low solar luminosities, the ice line lies in the Ferrel cell, stabilised by the poleward winds despite moderate meridional temperature gradients under strong greenhouse warming. For critical states at high solar luminosities, on the other hand, the ice line rests at the Hadley cell boundary, stabilised against the equatorward winds by steep meridional temperature gradients resulting from the increased solar energy input at lower latitudes and stronger Ekman transport in the ocean.

1 Introduction

Today, the climate of planet Earth is in a state which is neither too hot nor too cold for life, with the ice line being far from the Equator. The position of the ice line is ultimately determined by the planetary energy balance, particularly the incoming solar radiation and the concentration of greenhouse gases in the atmosphere, as well as the ice-albedo feedback and meridional heat transport in Earth's climate system (e.g. North et al., 1981).

One of the consequences of the positive ice-albedo feedback is the phenomenon that, over a range of boundary conditions, more than one state of Earth's climate can be stable for the same levels of incoming solar radiation and greenhouse gases. So even with the current solar luminosity and atmospheric composition, Earth might as well rest in a state that is a lot less welcoming to living beings – a fully glaciated Snowball state with surface temperatures multiple tens of degrees lower than today (e.g. North, 1990). There would be no liquid water at Earth's surface, which is at least a necessary condition for complex life as we know it.

For solar luminosities and/or greenhouse gas concentrations lower than today, a bifurcation point in phase space is reached at some point. This means that, for a certain set of parameters, the system ends up with only the Snowball state being stable. If the current climate cooled down, the ice line would descend slowly towards the Equator at first. But at some point, the system would tip; Earth would rapidly cool down, ending up in a global glaciation (Öpik, 1953). The Snowball instability crossed in this case is fundamental to the climate system and is ultimately caused by the positive ice-albedo feedback.

The first to quantify the Snowball instability in terms of the critical solar luminosities for given system parameters were Budyko (1969) and Sellers (1969), who analysed the Snowball bifurcation with one-dimensional energy balance models (EBMs; North et al., 1981). These EBMs are built upon the principle of energy conservation at a given latitude and consist of either a time-dependent differential equation for the temperature as a function of latitude or a time-independent equation assuming that the system is in equilibrium. Budyko and Sellers found that even a small decrease in solar luminosity would suffice for an Earth with today’s system parameters (in terms of the greenhouse effect) to end up in a state of complete glaciation.

Despite their simplicity, many authors have used EBMs – often incorporating various modifications and extensions – ever since the pioneering studies by Budyko (1969) and Sellers (1969), especially in order to understand how various factors influence the Snowball instability and the climate system's phase space properties (Faegre, 1972; Schneider and Gal-Chen, 1973; Held and Suarez, 1974; North, 1975a, b; Gal-Chen and Schneider, 1976; Ghil, 1976; Drazin and Griffel, 1977; Lindzen and Farrell, 1977; Cahalan and North, 1979; North and Coakley, 1979; North et al., 1981, 1983; North, 1990; Huang and Bowman, 1992; Ikeda and Tajika, 1999; Shen and North, 1999; Rose and Marshall, 2009; Roe and Baker, 2010).

In terms of Earth's long-term habitability, the Snowball bifurcation is particularly relevant in light of the fact that the solar luminosity was considerably lower in the past (e.g. Gough, 1981). The time evolution of the critical CO₂ concentration required to prevent global glaciation has typically been studied with radiative–convective models (RCMs; Ramanathan and Coakley, 1978) rather than EBMs. However, in RCMs, the ice-albedo feedback is not taken explicitly into account but has typically been considered by requiring a global mean surface temperature well above the freezing point of water. Examples for such investigations can be found in Kasting (1987) or von Paris et al. (2008), for example.

In addition to considerations of planetary habitability, open problems related to specific time periods in Earth's history are the second important reason why one should be interested in the Snowball instability. The most relevant geologic eons in this respect are the Archean (about 4 to 2.5 Ga – 1 Ga =10⁹ years ago), with evidence for habitable conditions on Earth despite a considerably fainter young Sun (Feulner, 2012; Charnay et al., 2020), and the Proterozoic (about 2.5 to 0.54 Ga), for which there is geological evidence for near-global glaciations both in the beginning (Tang and Chen, 2013) and towards the end of the eon (Hoffman and Schrag, 2002).

During the Archean, the luminosity of the young Sun was reduced by 20 %–25 % as compared to today. Results from EBMs and RCMs suggest that such a large reduction in the incoming solar radiation should have turned Earth into a Snowball, yet the geologic record suggests the presence of liquid water. This puzzling discrepancy, known as the faint young Sun paradox (Feulner, 2012; Charnay et al., 2020), has led to considerable interest in quantifying the Snowball bifurcation point on early Earth. For 4 decades, climate-modelling work on the faint young Sun paradox has been dominated by radiative–convective climate models that essentially neglect the ice-albedo feedback. More recently, however, a number of research groups have published first results from three-dimensional climate models investigating solutions for the faint young Sun paradox (Kienert et al., 2012, 2013; Charnay et al., 2013; Wolf and Toon, 2013; Kunze et al., 2014; Le Hir et al., 2014). However, all of these model studies are limited either by simplifications in their atmosphere or in their ocean and/or sea-ice components and use a variety of different assumptions and boundary conditions, leading to a large spread in simulated critical CO₂ concentrations.

In contrast to the Archean, there is geologic evidence for low-latitude glaciations during the Proterozoic (Tang and Chen, 2013; Hoffman and Schrag, 2002), so the central question becomes the following: what are the conditions required for Snowball events? This is in contrast to the question of what prevented Earth from freezing over completely. This problem has been investigated in a number of modelling studies, particularly for the Neoproterozoic (e.g. Chandler and Sohl, 2000; Lewis et al., 2003; Poulsen et al., 2002; Donnadieu et al., 2004; Lewis et al., 2007; Pierrehumbert et al., 2011; Voigt et al., 2011; Voigt and Abbot, 2012; Liu et al., 2013; Feulner and Kienert, 2014; Feulner et al., 2015; Liu et al., 2017; Braun et al., 2022; Hörner et al., 2022), yielding critical CO₂ concentrations from below 40 ppm to about 700 ppm depending on model type and boundary conditions.

In most climate-modelling studies, global glaciation is initiated once the ice margin moves to latitudes lower than about 30^∘. Some models, however, exhibit stable critical states with ice lines much closer to the Equator, often referred to as waterbelt or Jormungand states (e.g. Abbot et al., 2011). In these models, such states are stabilised by a number of mechanisms, most importantly the low albedo of bare sea ice and wind-driven meridional heat transport in the ocean towards the ice line (Voigt and Abbot, 2012). It should be noted, however, that the existence of stable states with tropical sea-ice margins often occurs in models lacking representation of sea-ice dynamics, which has been demonstrated to strongly destabilise such states (Voigt and Abbot, 2012). Recently, Braun et al. (2022) have shown that tropical waterbelt states are strongly influenced by cloud radiation and microphysics, which brings into question the stability of these states under Neoproterozoic boundary conditions.

Despite the relevance of the phenomenon for our understanding of past climate states and planetary habitability, there has – to our knowledge – not yet been an attempt to study the Snowball instability throughout Earth's history within one consistent framework and using a more complex climate model. At one end of the spectrum, earlier studies comprise conceptual investigations with very simple climate models like RCMs or EBMs. The latter, in particular, help to understand the principles of the instability and changes in phase space. At the other end of the spectrum, there are many investigations of single events in Earth's history with climate models of various complexities, ranging from models of intermediate complexity to atmosphere–ocean general circulation models (AOGCMs), and with a variety of different boundary conditions. These studies provide detailed insights about the time slices investigated, but the lack of uniform simulation design, model architecture, and boundary conditions makes it hard to compare them to each other or to study the evolution of the Snowball instability with time.

In this study, we want to bridge the gap between studies of the Snowball instability for single time slices with complex models and conceptual investigations of its time evolution. To this end, we use an Earth system model of intermediate complexity (EMIC; Claussen et al., 2002) in an aquaplanet configuration to scan for the Snowball bifurcation point for time slices spanning the last 4 billion years, thus quantifying the time evolution of the bifurcation and identifying a qualitative shift in critical-state dynamics. Although aquaplanet setups were used in the context of Snowball glaciations before (e.g. Pierrehumbert et al., 2011; Braun et al., 2022; Hörner et al., 2022), we uniquely focus on the long-term evolution of the bifurcation point.

This article is organised as follows. In Sect. 2, we describe the coupled climate model used to scan the Snowball bifurcation of aquaplanets through Earth history, as well as the boundary conditions and design of our numerical experiments. In Sect. 3, we present the Snowball bifurcation for aquaplanets through time and compare our results to earlier studies (Sect. 3.1). Furthermore, we describe global properties of the critical states for the different time slices (Sect. 3.2) and discuss the two different dynamical regimes for the critical state (Sect. 3.3). Finally, in Sect. 4, we summarise the major findings of our work in the context of earlier studies, discuss limitations, and outline potential future research.

2 Methods: coupled climate model simulations

2.1 Model description

Scanning for the Snowball bifurcation for more than a dozen time slices throughout Earth's history requires a relatively fast coupled climate model. We employ the Earth system model of intermediate complexity CLIMBER-3α (Montoya et al., 2005). CLIMBER-3α consists of a modified version of the ocean general circulation model (OGCM) MOM3 (Pacanowski and Griffies, 1999; Hofmann and Morales Maqueda, 2006), with a horizontal resolution of 3.75${}^{\circ }\times \mathrm{3.75}$^∘ and 24 vertical levels; a dynamic–thermodynamic sea-ice model (Fichefet and Morales Maqueda, 1997) at the same horizontal resolution, allowing for partially ice-covered grid cells; and a fast statistical–dynamical atmosphere model (Petoukhov et al., 2000) with a coarse horizontal resolution of 22.5^∘ in longitude and 7.5^∘ in latitude.

We emphasise that the sea-ice model explicitly takes into account sea-ice dynamics, a factor which has been found to be of crucial importance for the Snowball bifurcation (Lewis et al., 2003, 2007; Voigt and Abbot, 2012). The Snowball bifurcation also critically depends on cryosphere albedo (e.g. Yang et al., 2012a). Our model uses clear-sky albedo values (averaged over all wavelengths) of 0.50 and 0.40 for freezing and melting sea ice and 0.75 and 0.65 for cold and warm snow, respectively. Albedo values for ultraviolet+visible light are 0.30 larger than near-infrared albedos, and a partitioning of 60 % (ultraviolet+visible) and 40 % (near-infrared) is assumed. The effects of snow cover on sea ice are explicitly taken into account.

The main limitations of the model relate to its simplified atmosphere component (Petoukhov et al., 2000). Particularly relevant for this study are the coarse spatial resolution, the highly parameterised vertical structure, the simple two-layer cloud scheme with cloud fractions depending on humidity and vertical velocity, and the simplified description of large-scale circulation patterns, including the fixed annual-mean width of the Hadley and Ferrel cells. Note that the boundary between the Hadley cells moves with the thermal equator, with a corresponding but smaller shift in the boundaries between the Hadley and the Ferrel cells – see Petoukhov et al. (2000, Sect. 3.2). Thus the overall changes of the large-scale circulation with the seasonal cycle are represented in the model in principle.

We note, however, that despite these limitations, the Snowball bifurcation points derived for Neoproterozoic time slices with our model (Feulner and Kienert, 2014) fall well within the range of those from state-of-the-art atmosphere–ocean general circulation models (AOGCMs; see also Fig. 1) and agree very well with models using similar cryosphere albedo values (Voigt and Abbot, 2012; Yang et al., 2012c; Liu et al., 2013). The impact of model limitations on our results will be discussed in Sect. 4.

2.2 Boundary conditions and design of numerical experiments

To facilitate comparison between the different time slices, we chose an aquaplanet configuration without any continents. In contrast to some other coupled-model simulations of aquaplanets, we do not place small islands at the poles; the poles are treated similarly to the North Pole in present-day simulations with ocean models using spherical grids, applying filtering in the polar regions to prevent numerical instability. The ocean topography was generated by randomly assigning an ocean depth to each grid cell using a Gaussian probability distribution with a mean depth of 3000 m and a variance of 450 m. For each grid cell, the resulting random depth value was then assigned to the corresponding vertical level of the ocean model’s grid. We chose to have an ocean with varying depths rather than a flat ocean floor in order to avoid potential numerical instabilities.

The Snowball bifurcation point is derived for a total of 18 time slices ranging from today to 3600 Myr ago (million years ago) – see Table 1. The solar constant was scaled based on the approximation formula from Gough (1981) assuming a present-day value of S₀=1361 W m⁻² (Kopp and Lean, 2011) and an age of the Sun of 4570 Myr (Bonanno et al., 2002). Orbital parameters were set to a circular orbit with obliquity of 23.5^∘. For each time slice, a number of equilibrium simulations were run for different CO₂ concentrations bracketing the Snowball bifurcation (see Table 1). In addition, we have run model experiments at two fixed levels of CO₂ (10 000 and 10 ppm) and decreasing solar luminosities of 1140, 1130, 1125, and 1120 W m⁻² for 10 000 ppm, and 1334, 1329, and 1324 W m⁻² for 10 ppm, respectively. For the lowest value of the solar constant in each of these cases, the model entered a Snowball state.

The total atmospheric pressure was kept constant at 1 bar in all simulations. Most simulations were initialised from a warm, ice-free state with idealised symmetric present-day ocean temperature and salinity fields. Note that critical-state characteristics might depend on initial conditions (e.g. Yang et al., 2012a). Our results are valid for trajectories starting from a warm, ice-free state; other initial conditions are not investigated here. In many cases, simulations pinpointing the Snowball bifurcation point were branched from runs with higher CO₂ concentrations. All simulations were integrated for at least 2000 model years after the last change in CO₂ concentration to allow the ocean to approach equilibrium conditions.

Table 1 Overview of simulation experiments with the age of each time slice, the value of the solar constant S used in the simulations, and the ratio $S/S_{\mathrm{0}}$ to the present-day value S₀, and the atmospheric CO₂ concentration of not-fully-glaciated and Snowball states.

3 Results

3.1 The Snowball bifurcation for aquaplanets through time

The critical CO₂ concentration for the Snowball bifurcation through Earth's history as derived from the aquaplanet simulations with CLIMBER-3α is shown in Fig. 1. As expected, there is excellent agreement between the bifurcation points derived at the two fixed CO₂ levels and the closest corresponding simulation derived at fixed values of the solar constant. In other words, there is no fundamental difference between scanning for the critical values in the horizontal and the vertical direction in the diagram. In the figure, the values are also compared to proxy estimates for the past CO₂ concentration in Earth's atmosphere and to earlier modelling studies for individual time slices, the latter being differentiated by model type.

Figure 1 Snowball bifurcation (in terms of CO₂ partial pressure assuming 1 bar total atmospheric pressure) as a function of solar luminosity for the aquaplanet simulations presented in this work (large blue circles). The open circles indicate the bifurcation points in solar luminosity for the additional model experiments at fixed levels of 10 000 and 10 ppm of CO₂. The red line shows a fit of the function defined in Eq. (1) to these bifurcation points – see the main text for details. The smaller black symbols denote earlier modelling studies with AOGCMs (squares), OGCMs coupled to a simplified atmosphere model (circles), and AGCMs coupled to a simplified ocean model (triangles). The model results are averages between the highest CO₂ concentration (or solar constant) of fully glaciated states and the lowest value of states with open water, with the range between these two values indicated by the error bars unless the latter are smaller than the symbols. Note that some model studies provide upper limits (indicated by arrows) rather than ranges. In cases where the contributions of other greenhouse gases (particularly methane) were included in the models, we show the effective CO₂ concentration calculated from the combined forcing. All model results are scaled to a modern solar constant of 1361 W m⁻² (Kopp and Lean, 2011). Estimates of past atmospheric CO₂ levels are shown in grey; Phanerozoic CO₂ values are taken from Foster et al. (2017), and the modern value is shown as a diamond. Global glaciations in Earth's history (light-blue shading) occurred during the Paleoproterozoic and Neoproterozoic eras. Proxy data: R95 – Rye et al. (1995), H04 – Hessler et al. (2004), S06 – Sheldon (2006), R10 – Rosing et al. (2010), D11 – Driese et al. (2011), KM15 – Kanzaki and Murakami (2015), L+20 – Lehmer et al. (2020), ST20 – Strauss and Tosca (2020). Model studies: CS00 – Chandler and Sohl (2000), L+03 – Lewis et al. (2003), P+02 – Poulsen et al. (2002), R+06 – Romanova et al. (2006), VM10 – Voigt and Marotzke (2010), P+11 – Pierrehumbert et al. (2011), V+11 – Voigt et al. (2011), K+12 – Kienert et al. (2012), VA12 – Voigt and Abbot (2012), Y+12a – Yang et al. (2012a), Y+12c – Yang et al. (2012c), C+13 – Charnay et al. (2013), L+13 – Liu et al. (2013), L+14 – Le Hir et al. (2014), WT13 – Wolf and Toon (2013), FK14 – Feulner and Kienert (2014), K+14 – Kunze et al. (2014), T+14 – Teitler et al. (2014), F17 – Feulner (2017), FS17 – Fiorella and Sheldon (2017), L+17 – Liu et al. (2017), B+22 – Braun et al. (2022), H+22 – Hörner et al. (2022).

With solar luminosity increasing over time, the critical CO₂ concentration in our aquaplanet simulations falls from ∼10⁵ ppm to below 1 ppm between 4 Ga and today. This decrease is more or less logarithmic from 4 to 1 Ga, with more steeply falling CO₂ levels in more recent times. The logarithmic decrease with linearly increasing insolation can be understood in terms of the well-known logarithmic behaviour of the CO₂ radiative forcing (Huang and Bani Shahabadi, 2014). The downturn of the critical CO₂ concentration for solar luminosities approaching the modern value can be attributed to the following three factors: (1) the steadily increasing incoming solar radiation (which shows a strong variation with latitude) leads to a more positive surface radiation balance, particularly over the open-ocean areas, despite a relatively weak greenhouse warming (which is more uniformly distributed), making sea-ice formation at low latitudes increasingly difficult. (2) Even at the very low global mean surface air temperatures of $\sim -\mathrm{15}{}^{\circ }$C of the critical states (see Fig. 2), there is greenhouse warming due to water vapour, which becomes significant at very low CO₂ levels. This is also facilitated by the fact that the maximum of the thermal-emission spectrum is shifted towards longer wavelengths and thus into the H₂O rotational bands because of the lower temperature. (3) Finally, there is also a warming contribution from cloud radiative effects.

The critical CO₂ concentrations as a function of the solar constant S derived from our aquaplanet simulations can be approximated by the following formula (S₀ is the present-day solar constant):

$$\begin{array}{}\text{(1)}& p{\mathrm{CO}}_{\mathrm{2},\mathrm{crit}}=a_{\mathrm{1}}\exp \left[a_{\mathrm{2}}{\left(a_{\mathrm{3}}-{\displaystyle \frac{S}{S_{\mathrm{0}}}}\right)}^{a_{\mathrm{4}}}\right].\end{array}$$

Fitting this function to the points derived from the aquaplanet simulations yields the following parameters: $a_{\mathrm{1}}=(\mathrm{0.0836}\pm \mathrm{0.0721})$ ppm, $a_{\mathrm{2}}=\mathrm{26.6}\pm \mathrm{0.7}$, $a_{\mathrm{3}}=\mathrm{1.0}\pm \mathrm{0.00002}$, and $a_{\mathrm{4}}=\mathrm{0.475}\pm \mathrm{0.051}$, providing a good approximation over the entire range of solar luminosities – see Fig. 1.

Before discussing our results in the context of previous model studies for several key time periods in Earth's history in Sect. 3.1.2, we will describe and discuss more general features which can be derived from this synthesis and earlier studies.

3.1.1 General observations on the Snowball bifurcation

The presence of continents and ice sheets makes global glaciations easier. For models of similar design, the presence of continents makes Earth more susceptible to glaciation as compared to an aquaplanet configuration, predominantly due to the higher surface albedo of land areas compared to open oceans and the lower water vapour content of the atmosphere (Poulsen et al., 2002; Voigt et al., 2011; Liu et al., 2013; Kunze et al., 2014). The aquaplanet simulation for a modern solar constant and 277 ppm of CO₂, for example, has a global and annual mean surface air temperature of 19.4 ^∘C compared to 15.1 ^∘C in a pre-industrial simulation using our model (Feulner, 2011). Similarly, our simulations with Neoproterozoic continents (Feulner and Kienert, 2014) indicate slightly higher critical CO₂ values than the aquaplanet simulation at similar solar luminosity, although the difference is small in this case due to the low albedo of bare land assumed in Feulner and Kienert (2014). Note that, in the extreme case of a fully land-covered planet, global glaciation becomes more difficult because the drier atmosphere leads to reduced cloud and snow cover and thus a lower albedo compared to the aquaplanet case (Abe et al., 2011).

It has also been shown already that the presence of tropical ice sheets shifts the Snowball bifurcation point to values ∼ 3–10 times higher (Liu et al., 2017) than the case without ice sheets (Liu et al., 2013). The combined effect of continents and polar ice sheets is also the most likely cause for the lower critical CO₂ concentrations of the aquaplanet simulations for modern boundary conditions (Yang et al., 2012a, b, c) and the Late Paleozoic ice age (Feulner, 2017). In addition, the simulations in Feulner (2017) were carried out for a glacial orbital configuration rather than for the circular orbit used in the aquaplanet simulations.

Meridional ocean heat transport makes global glaciation more difficult. It is evident from Fig. 1 that atmospheric general circulation models (AGCMs) without ocean heat transport consistently show bifurcation points at higher CO₂ levels than models with prescribed ocean heat transport or dynamic ocean models. This is in line with earlier findings showing that ocean heat transport towards the sea-ice edge, particularly by the wind-driven ocean circulation, makes Snowball initiation harder (Poulsen and Jacob, 2004; Voigt and Abbot, 2012; Rose, 2015).

Models without sea-ice dynamics are too stable. Studies carried out with AGCMs coupled to mixed-layer ocean models with prescribed ocean heat transport but without dynamic sea ice tend to predict lower values for the Snowball bifurcation point (see Fig. 1). Indeed, the fact that models without sea-ice dynamics are artificially stable with respect to the Snowball bifurcation has been noted before (Lewis et al., 2003, 2007; Voigt and Abbot, 2012).

The glaciation threshold depends on sea-ice and snow albedo. Even for models of the same design, there is considerable spread in the values for the Snowball bifurcation point for similar boundary conditions (see Fig. 1). Differences in cloud radiative forcing and simulated heat transport in the atmosphere and the oceans can contribute to this spread; however, the predominant causes are differences in sea-ice and snow albedo values and parameterisation (Pierrehumbert et al., 2011; Yang et al., 2012c). These have been identified as the cause of the difference between the bifurcation points derived with CCSM3 (Yang et al., 2012a) and CCSM4 (Yang et al., 2012c) for the same present-day boundary conditions, for example.

3.1.2 Comparison with earlier modelling studies for selected time periods

Modern Snowballs. The Snowball bifurcation point for modern boundary conditions has been quantified with an AGCM in terms of reduced CO₂ by Romanova et al. (2006) and with AOGCMs in terms of a reduced solar constant by Voigt and Marotzke (2010) and by Yang et al. (2012a, b, c). For pre-industrial greenhouse gas concentrations, the AOGCM studies put the bifurcation point in the ranges of 91 %–94 % (Voigt and Marotzke, 2010), 89.5 %–90 % (Yang et al., 2012a), and 91 %–92 % (Yang et al., 2012c) of the present-day solar constant, comparing well to about 91 % of today's solar constant in our simulations (see Fig. 1). For a reduction in CO₂ and a fixed present-day solar constant, the bifurcation point in our aquaplanet simulations is below a CO₂ concentration of 0.1 ppm. This agrees well with the AGCM study by Romanova et al. (2006), where their model is not in a Snowball state at their lowest CO₂ concentration of 1 ppm. Voigt and Marotzke (2010) used ECHAM5/MPI-OM, finding a Snowball state at 0.1 ppm of CO₂ for present-day continents and a slightly higher value of the solar constant of 1367 W m⁻², which is again in good agreement with our results given the cooling influence of continents and the generally higher susceptibility to global glaciation of their AOGCM.

Proterozoic. For the Neoproterozoic, earlier model studies indicate critical CO₂ concentrations from below 40 ppm in an AGCM simulation (Chandler and Sohl, 2000) to about 100–700 ppm in AOGCMs and OGCMs with sea-ice dynamics (Poulsen et al., 2002; Voigt et al., 2011; Voigt and Abbot, 2012; Liu et al., 2013; Feulner and Kienert, 2014; Liu et al., 2017). The reason for the higher value found by Lewis et al. (2003) remains unclear, but could be connected to the fact that surface winds are prescribed in their model. In their studies on modern Snowballs, Yang et al. (2012a, c) find critical CO₂ concentrations of about 20 ppm to 100 ppm, depending on the model version, for a solar constant representative of the Neoproterozoic, i.e. reduced by 6 % compared to its modern value. In our aquaplanet configuration, the corresponding value at 700 Ma is 95±5 ppm and is thus well within the typical range of the models with ocean and sea-ice dynamics. With the exception of Chandler and Sohl (2000) and Yang et al. (2012a), our value is somewhat lower than the one derived in other studies, as one would expect for a configuration without continents or ice sheets and for the snow and sea-ice albedo values employed in our model (see also Sects. 2.1 and 3.1.1). Furthermore, there is also excellent agreement with the bifurcation point quantified using an AOGCM for an earlier time slice at 1 Ga by Fiorella and Sheldon (2017), where our values for 900 and 1200 Ma are again compatible with the lower range of their estimate – see Fig. 1.

Archean. In many ways, the situation is most complicated for the Archean, where the Snowball bifurcation has been quantified in a number of model studies in the context of the faint young Sun paradox (Feulner, 2012; Charnay et al., 2020). So far, there are no AOGCM studies for the Archean. The bifurcation points from Kienert et al. (2012), quantified with an earlier version of the model used in this work, are higher than the ones for the aquaplanet configuration, mainly caused by the higher cryosphere albedo values used in Kienert et al. (2012) and the different boundary conditions, particularly the higher rotation rate. All other model studies use AGCMs with simplified ocean components and sea-ice models without dynamics (Charnay et al., 2013; Wolf and Toon, 2013; Kunze et al., 2014; Le Hir et al., 2014; Teitler et al., 2014). The critical CO₂ values from these studies are generally (and in many cases, significantly) below the values derived for the aquaplanets in this paper. While we cannot rule out a contribution from the simplified atmosphere component used in our model, the lack of full ocean and, in particular, sea-ice dynamics in the other studies is likely to be a significant factor in explaining this discrepancy. In the end, this question can only be fully answered by running AOGCM simulations with Archean boundary conditions, which is beyond the scope of the present work.

3.2 Global characteristics of aquaplanet critical states

In the following, we will have a more detailed look at the global characteristics of the climate states close to the Snowball bifurcation, i.e. the stationary states for all investigated solar luminosities with the lowest concentration of CO₂ in the atmosphere at which the system does not fall into a Snowball state.

Figure 2 (a) Global annual mean surface air temperature and (b) sea-ice fraction of critical states as a function of solar luminosity (blue circles). The right-hand vertical axis in (b) indicates the effective latitude of the critical ice line, calculated from the sea-ice fraction assuming full ice cover over two symmetric polar caps, which is a reasonable approximation for the aquaplanet critical states (see Fig. 5). The grey shading indicates the ±3σ ranges of the respective values. Note that the states with higher temperature and lower sea-ice fraction are also stable states at higher solar luminosities, as shown for the 900 Ma time slice in Fig. 3.

Let us first take a look at the global annual mean surface air temperature and the mean sea-ice fraction of the critical states for the different time slices shown in Fig. 2, where the averages are taken over the last 100 years of each simulation. Whereas the critical CO₂ concentrations as a function of solar luminosity would indicate a fairly smooth change over time (see Fig. 1), there is a marked shift in both the global annual mean surface air temperature and the mean sea-ice fraction; critical states at low solar luminosities consistently have global mean temperatures of about $-\mathrm{2}{}^{\circ }$C and average sea-ice fractions of about 33 %, whereas critical states at higher solar luminosities exhibit much lower temperatures of about $-\mathrm{15}{}^{\circ }$C and higher sea-ice fractions of 50 %. In addition, there is a slight cooling trend with increasing solar luminosities for the critical states at both lower and higher solar luminosities. Finally, the sequence of critical states with higher sea-ice fractions shows remarkably little scatter in this quantity – see Fig. 2b. Note that the shift in critical-state properties is not related to the downturn of the bifurcation limit at higher solar luminosities discussed in Sect. 3.1.

Figure 3 Global and annual mean sea-ice fraction in equilibrium climate states at 900 Ma for various concentrations of atmospheric CO₂. Grey shaded areas are as in Fig. 2. Note that the simulation at 230 ppm is on track to a Snowball state but has not reached equilibrium due to numerical instability.

The values of 33 % and 50 % for the global annual mean sea-ice fraction found for the critical states can be translated into latitudes of the sea-ice margin of 42 and 30^∘, respectively, if one assumes full ice cover over two symmetric polar caps, which is a reasonable assumption for aquaplanets. Thus, for lower solar luminosities, the sea-ice margin rests within the Ferrel cell of the atmosphere's large-scale circulation, whereas for higher solar luminosities, it corresponds to the Hadley cell boundary (which is fixed at 30^∘ in our simplified atmosphere model). We will therefore refer to the critical states at lower solar luminosity as Ferrel states and to the ones at higher solar luminosity as Hadley states for simplicity.

In Fig. 2, we focus on the critical states only, i.e. the coldest not fully glaciated state at each solar luminosity. While we did not find any Hadley states at low solar luminosities, the question remains whether there could be Ferrel states at higher solar luminosities. This can be investigated by looking at equilibrium states for a range of CO₂ concentrations for one particular time slice at higher solar luminosity. The result for the 900 Ma time slice is shown in Fig. 3, where one can see that, even for higher solar luminosities, Ferrel states can indeed be stable states at higher atmospheric CO₂ levels, as already indicated by the grey shading in Fig. 2.

On the other hand, we can also look for transient Hadley states in simulations for Ferrel state time slices where the system falls into the Snowball state. In these simulations, the unstable Hadley state can clearly be identified as a time period with close to 50 % global and annual mean sea-ice cover f_sea ice. In agreement with the findings of Fig. 2, these states (defined by $\mathrm{0.497}<f_{\mathrm{sea}\mathrm{ice}}<\mathrm{0.503}$) become more stable with increasing solar luminosity. While the system spends only 70 years in the transient Hadley state at 3600 Ma, for example, this length more than quintuples to 396 years at 1350 Ma (Fig. 4).

Figure 4 Transient Hadley states with a global mean sea-ice cover of about 50 % in simulations on track to a Snowball. The transient Hadley phase can be clearly seen as a plateau in all simulations. To capture the entry into and the exit from this state for all time slices, we use a range of $\mathrm{0.497}<f_{\mathrm{sea}\mathrm{ice}}<\mathrm{0.503}$ to define the transient Hadley phase, indicated by the horizontal lines.

In summary, we find a marked shift in the characteristics of the critical states in the time period from 1350 to 1200 Ma. While a Ferrel state with the sea-ice margin in the middle of the Ferrel cell is stable for all luminosities, a Hadley state with the ice line at the Hadley cell boundary can only be stable for solar luminosities above about 90 % of the modern solar constant. The fact that the critical state at 1800 Ma is a Hadley state already, whereas the critical states at 1650, 1500, and 1350 Ma are Ferrel states, could again be due to the sensitivity of the system with respect to small changes, particularly in the transition period.

Figure 5 Maps of annual mean surface air temperatures (a, c) and sea-ice fractions (b, d) for the critical states at 3000 Ma (a, b) and 900 Ma (c, d). Only one-quarter of the globe is displayed in each map, since the aquaplanet simulations exhibit a high degree of symmetry. Surface wind velocity vectors are shown in the right-hand panels, with the length scaling for a wind speed of 10 m s⁻¹ indicated in the bottom right corner.

3.3 Dynamics of the Snowball bifurcation of aquaplanets

In order to understand the different regimes of the critical states, we will have a more detailed look at typical Ferrel and Hadley states. To this end, we select the critical states at 3000 Ma (a Ferrel state) and 900 Ma (a Hadley state) and investigate the spatial patterns of surface air temperature, sea-ice fraction, and surface winds (Fig. 5). Maps of annual mean surface air temperatures show that the Ferrel state is considerably warmer at all latitudes (Fig. 5a, c). The annual mean sea-ice distributions exhibit marked differences as well. The Hadley state at 900 Ma has a very well defined, sharp sea-ice margin, whereas the Ferrel state at 3000 Ma is characterised by a much more fuzzy transition between open ocean and full ice cover (Fig. 5b, d). In both critical states, the wind fields show the usual pattern of the large-scale atmospheric circulation, with a high degree of symmetry about the Equator, as one would expect for aquaplanet boundary conditions.

Figure 6 Zonal means of (a) annual mean surface air temperatures, (b) sea-ice fractions, and (c) meridional wind speed for the Ferrel state at 3000 Ma (solid red lines) and the Hadley state at 900 Ma (dashed red line). The vertical lines indicate the sea-ice margin for the Ferrel state (dashed blue line) and for the Hadley state (solid blue line).

These findings are even more evident from the zonal averages of annual mean surface air temperature, sea-ice fraction, and the meridional component of the surface wind for the two critical states shown in Fig. 6. In particular, we would like to highlight (1) the steeper temperature gradient across the sea-ice margin in the Hadley state, (2) the position of the sea-ice margin in the Hadley state close to the point where the meridional surface winds change their direction, and (3) the latitude of the sea-ice margin in the Ferrel state close to the poleward maximum of the meridional surface wind field.

Figure 7 Schematic illustration of the differences between the critical states (a) at lower solar luminosities (Ferrel states) and (b) at higher solar luminosities (Hadley states) showing sea-ice thickness (cyan – not to scale) on the ocean (blue) and the large-scale atmospheric wind patterns (arrows); see the text for the discussion. The vertical cyan lines indicate the effective ice line latitudes calculated assuming symmetric and full ice cover, as in Fig. 2.

Indeed, the specific locations of the sea-ice boundaries close to the centre of the Ferrel cell and at the edge of the Hadley cell immediately suggest that they are influenced by large-scale atmospheric circulation patterns, as illustrated in Fig. 7. The Ferrel states are stabilised by poleward winds pushing the sea-ice margin away from the Equator and transporting relatively warm air towards the ice edge. The Hadley states, on the other hand, are the coldest possible not-fully-glaciated states because any further cooling would bring sea ice into the Hadley cell, where it would be pushed towards the Equator by the trade winds, which would also transport cold air towards the sea-ice margin. Since, in the case of the Hadley states, the wind fields are a destabilising factor, these states have to be protected from global glaciation by relatively high surface temperatures within the Hadley cells, causing any sea ice which enters the Hadley cells to melt rapidly.

In the following, we will investigate this hypothesis in more detail, particularly with respect to the important question of why the Hadley states can be stable only at higher solar luminosities. Answering this question requires, among other things, a closer look at the evolution of the energy balance of the Hadley states. In particular, we would like to understand how and why the surface temperature distribution of the Hadley states changes with increasing solar luminosity. To this end, we follow Heinemann et al. (2009), Voigt et al. (2011), and Liu et al. (2013) in applying the following simple equation for the annual mean equilibrium energy balance at each latitude φ:

$$\begin{array}{}\text{(2)}& Q\left(\mathit{\phi }\right)(\mathrm{1}-\mathit{\alpha }(\mathit{\phi }\left)\right)=Q_{\mathrm{abs}}\left(\mathit{\phi }\right)=\mathit{\sigma }\mathit{\epsilon }\left(\mathit{\phi }\right)T_{\mathrm{sfc}}^{\mathrm{4}}\left(\mathit{\phi }\right)+\mathrm{div}F\left(\mathit{\phi }\right).\end{array}$$

This is done in order to disentangle the contributions of changes in absorbed solar radiation, greenhouse warming, and meridional heat transport to the surface temperature evolution of the Hadley states. In this equation, Q is the incoming solar radiation flux at latitude φ, α is the top-of-atmosphere albedo, Q_abs is the absorbed solar radiation flux (directly diagnosed from model output), σ is the Stefan–Boltzmann constant, ε is the effective emissivity (diagnosed from the ratio of top-of-atmosphere to surface upward long-wave fluxes), T_sfc is the surface temperature, and divF is the divergence of the total meridional heat transport (diagnosed from the net top-of-atmosphere radiation balance). For a given equilibrium simulation, the surface temperature T_sfc,0 can then be derived by simply solving Eq. (2) for the surface temperature, showing good agreement with the surface temperature directly diagnosed from model output (not shown).

The different contributions to surface temperature changes ΔT_sfc between a given equilibrium state and a reference state, denoted by subscript 0, can then be calculated as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{sfc},\mathrm{solar}}\left(\mathit{\phi }\right)& ={\left({\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }{\mathit{\epsilon }}_{\mathrm{0}}\left(\mathit{\phi }\right)}}\left(Q_{\mathrm{abs}}\right(\mathit{\phi })-\mathrm{div}{F}_{\mathrm{0}}(\mathit{\phi }\left)\right)\right)}^{\mathrm{1}/\mathrm{4}}\\ & -T_{\mathrm{sfc},\mathrm{0}}\left(\mathit{\phi }\right),\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{sfc},\mathrm{greenhouse}}\left(\mathit{\phi }\right)& ={\left({\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }\mathit{\epsilon }\left(\mathit{\phi }\right)}}\left(Q_{\mathrm{abs},\mathrm{0}}\right(\mathit{\phi })-\mathrm{div}{F}_{\mathrm{0}}(\mathit{\phi }\left)\right)\right)}^{\mathrm{1}/\mathrm{4}}\\ & -T_{\mathrm{sfc},\mathrm{0}}\left(\mathit{\phi }\right),\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{sfc},\mathrm{transport}}\left(\mathit{\phi }\right)& ={\left({\displaystyle \frac{\mathrm{1}}{\mathit{\sigma }{\mathit{\epsilon }}_{\mathrm{0}}\left(\mathit{\phi }\right)}}\left(Q_{\mathrm{abs},\mathrm{0}}\right(\mathit{\phi })-\mathrm{div}F(\mathit{\phi }\left)\right)\right)}^{\mathrm{1}/\mathrm{4}}\\ & -T_{\mathrm{sfc},\mathrm{0}}\left(\mathit{\phi }\right).\end{array}\end{array}$$

The results of this exercise are presented in Fig. 8, where we show the zonally averaged surface temperature differences between all Hadley critical states and the 900 Ma Hadley state together with the contributions of the three factors.

Figure 8 (a) Total difference of zonally averaged surface temperatures as diagnosed using a one-dimensional energy balance equation (see text) between the Hadley states from 1800 to 0 Ma and the one at 900 Ma. The other panels show the contributions of changes in absorbed solar radiation (b), greenhouse warming (c), and the combined atmospheric and oceanic meridional heat transport (d). The different Hadley states are indicated by increasingly darker shades of grey for increasing solar luminosities. Note the different scale of the vertical axis in panel (a).

It is evident from Fig. 8a that there is a distinct geographic pattern in the surface temperature changes of the Hadley states over time. With increasing solar luminosity (i.e. going from 1800 to 0 Ma), there is a gradual cooling in areas covered by sea ice and a (less pronounced) warming in ice-free regions. These trends are driven by the different spatial characteristics of the absorbed solar radiation and the warming due to the greenhouse effect (already noted in earlier work, e.g. Yang et al., 2012a). While the radiative fluxes due to greenhouse gases are distributed relatively uniformly with latitude, the absorbed solar energy shows a marked maximum at the Equator and drops off towards the poles, due to both the spatial distribution of the incoming solar radiation and the higher albedos of the ice-covered middle and high latitudes. Going from 1800 to 0 Ma, solar luminosity increases, while the CO₂ concentrations of the Hadley critical states decrease, leading to the evolution of their respective contributions to the surface temperature trends shown in Fig. 8b and c. These changes are only partially compensated for by adjustments of the meridional heat transport – see Fig. 8d.

The surface temperature evolution of the Hadley states described above also provides a qualitative explanation of why Hadley states cannot be stable at low solar luminosities. Going back in time, the decreasing solar radiation is compensated for by an increasing greenhouse warming to prevent global glaciation, but the different spatial patterns of these factors lead to progressively shallower surface temperature gradients across the Hadley cell boundary. This effect is further enhanced by a weakening of the meridional heat transport with decreasing solar luminosity, driven by a slowdown of the trade winds leading to weaker Ekman transport in the ocean. At some point, the surface energy budget within the Hadley cell is insufficient to melt sea ice before it is pushed towards the Equator by the trade winds, triggering global glaciation due to the ice-albedo feedback. This is also in line with the observed shortening of the transient Hadley phase with decreasing solar luminosity in simulations on track to a Snowball state, as shown in Fig. 4.

4 Discussion and conclusions

In this paper, we have investigated the Snowball bifurcation point (in terms of atmospheric CO₂ concentration) of aquaplanets under the steady increase of solar luminosity over Earth's history in one consistent model framework. We find that, until about 1 billion years ago, the critical CO₂ concentration decreased more or less logarithmically (from ∼10⁵ ppm 4 billion years ago to ∼10² ppm 700 million years ago), as one would expect from the logarithmic character of the CO₂ forcing. In more recent times, critical values decreased more strongly, mainly due to the increasing low-latitude insolation making sea-ice formation more difficult and due to the greenhouse effect of water vapour. We have also put these new simulation results in context by presenting a comprehensive synthesis of proxy CO₂ estimates and findings from earlier modelling studies.

The second important new conclusion from our work is of a regime shift in critical-state properties about 1.2 billion years ago. While the coldest not-fully-glaciated climate states at earlier times (and thus lower solar luminosities) are characterised by relatively high global mean temperatures of about $-\mathrm{2}{}^{\circ }$C and a sea-ice margin close to the centre of the Ferrel cell (Ferrel states), critical states at later times (and thus higher solar luminosities) exhibit much lower global mean temperatures of roughly $-\mathrm{15}{}^{\circ }$C and a sea-ice margin at the Hadley cell boundary (Hadley states). These states result from the interplay of the surface energy balance and atmospheric dynamics; in the Ferrel states, the sea-ice margin is stabilised by the poleward push of the meridional winds, whereas in the Hadley states, the sea-ice margin has to be maintained by a steep temperature gradient across the Hadley cell boundary to prevent sea-ice being pushed towards the Equator by the trade winds within the Hadley cells. The ultimate cause for the two distinct critical-state regimes are the different spatial distributions of solar radiation and greenhouse gas forcings. With increasing solar luminosity and decreasing CO₂ concentrations, this difference in the spatial distributions will result in ever-steeper temperature gradients, making Hadley states stable at some point.

While our investigation of the glaciation threshold through time in one consistent three-dimensional modelling framework and the regime shift from Ferrel to Hadley states is novel, aspects of our work tie in well with earlier findings. The importance of the Hadley circulation for global glaciations, for example, has been noted already by Bendtsen (2002) based on simulations with a simple coupled model. Furthermore, in their studies of modern Snowballs with CCSM3 and CCSM4, Yang et al. (2012a, b, c) find two different critical positions for the ice line depending on whether they reduce the solar constant at modern CO₂ concentrations or at the CO₂ concentrations at 94 % of the present-day solar constant.

We would like to point out that both the values for the Snowball bifurcation points for the different time slices and the solar luminosity at which the regime shift from Ferrel to Hadley states occurs will depend on model physics (e.g. the snow and cloud schemes) and model parameters (e.g. the cryosphere albedos) and will thus most likely differ between models. Furthermore, we emphasise that aspects of our study are certainly affected by model limitations like the coarse spatial resolution, the simple cloud scheme, or possible deviations of the radiative transfer at very low or very high CO₂ concentrations. Most importantly, however, the Hadley cells have a prescribed annual mean width of 30^∘ in our simplified atmosphere model, whereas it is well known that they become narrower in colder climate states (e.g. Frierson et al., 2007).

Despite these limitations, we consider the main finding of our paper, the regime shift from Ferrel states at lower solar luminosities to Hadley states at higher solar luminosities, to be robust. Indeed, there are indications in studies with more complex models supporting this conclusion. Most importantly and as mentioned above, Yang et al. (2012a) and Yang et al. (2012c) investigate modern Snowballs with CCSM3 and CCSM4, respectively, and find that “it is likely that there are no stable states between ∼ 40 % and 60 % sea ice coverage during the initiation of the Snowball Earth” (Yang et al., 2012a, p. 2719) for a present-day continental configuration. Note that their values for the critical global sea-ice fractions are somewhat higher than ours (33 % and 50 %; see Sect. 3.2), which could partly be due to the differences in boundary conditions but most likely reflects the fact that, in cold climate states, the Hadley cells will be narrower than the fixed annual-mean width of 30^∘ used in our simplified atmosphere model (see above).

Moreover, Yang et al. (2012a) also report different critical states in CCSM3 depending on the path to global glaciation. For a reduction of the solar constant at present-day CO₂ levels, their model enters a Snowball state beyond a global sea-ice fraction of ∼40 %, corresponding to a Ferrel state in our simulations. For a reduction of the atmospheric CO₂ concentration at 94 % of the present-day solar constant (mimicking Neoproterozoic boundary conditions) on the other hand, the critical sea-ice fraction is ∼60 %, similar to a Hadley state in our simulations. Yang et al. (2012a) attribute this difference in the paths to global glaciation to the dissimilar spatial characteristics of solar radiation and CO₂ forcing, in line with our findings above. Looking at Fig. 1, this would put the regime shift from Ferrel to Hadley states in Yang et al. (2012a) somewhere between ∼1.2 and ∼0.7 billion years ago, which is in excellent agreement with our findings despite the different boundary conditions. CCSM4 is more sensitive with respect to global glaciation than CCSM3 (Yang et al., 2012c), and in this model, the critical state is a Hadley state for both paths to global glaciation. This would put the transition from Ferrel to Hadley states to times earlier than ∼1 billion years ago, which is again in principal agreement with our results.

One conclusion from our work is that the Snowball bifurcation (and thus one of the most important limits to planetary habitability) is determined by the interplay of the energy balance and internal dynamics, implying the need for knowledge of certain system parameters and for three-dimensional models in order to be able to assess planetary habitability. Note that our results cannot be easily generalised to other planets for the following three reasons: first, the idealised aquaplanet configuration is a rather special case; second, our simulations have been performed for one particular orbital configuration; and third, the results will likely depend on the planet's rotation rate, since the structure of the Hadley circulation, for example, changes significantly with the rotation rate (e.g. Schneider, 2006). Sampling the parameter space more completely remains to be done in future work. Finally, our results highlight once again the crucial importance of sea-ice dynamics for investigations of the Snowball bifurcation point, for example in the context of the faint young Sun paradox or the Paleoproterozoic and Neoproterozoic glaciations.