Articles | Volume 14, issue 3
Research article
03 May 2023
Research article |  | 03 May 2023

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

Georg Feulner, Mona Bukenberger, and Stefan Petri

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 (CO2) at which this Snowball bifurcation occurs changes through Earth's history, most notably because of the slowly increasing solar luminosity. Quantifying this critical CO2 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 CO2 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 CO2 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. North1990). 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 (Öpik1953). 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 (Faegre1972; Schneider and Gal-Chen1973; Held and Suarez1974; North1975a, b; Gal-Chen and Schneider1976; Ghil1976; Drazin and Griffel1977; Lindzen and Farrell1977; Cahalan and North1979; North and Coakley1979; North et al.1981, 1983; North1990; Huang and Bowman1992; Ikeda and Tajika1999; Shen and North1999; Rose and Marshall2009; Roe and Baker2010).

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. Gough1981). The time evolution of the critical CO2 concentration required to prevent global glaciation has typically been studied with radiative–convective models (RCMs; Ramanathan and Coakley1978) 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 =109 years ago), with evidence for habitable conditions on Earth despite a considerably fainter young Sun (Feulner2012; 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 Chen2013) and towards the end of the eon (Hoffman and Schrag2002).

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 (Feulner2012; 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 Toon2013; 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 CO2 concentrations.

In contrast to the Archean, there is geologic evidence for low-latitude glaciations during the Proterozoic (Tang and Chen2013; Hoffman and Schrag2002), 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 Sohl2000; 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 Abbot2012; Liu et al.2013; Feulner and Kienert2014; Feulner et al.2015; Liu et al.2017; Braun et al.2022; Hörner et al.2022), yielding critical CO2 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 Abbot2012). 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 Abbot2012). 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 Griffies1999; Hofmann and Morales Maqueda2006), with a horizontal resolution of 3.75×3.75 and 24 vertical levels; a dynamic–thermodynamic sea-ice model (Fichefet and Morales Maqueda1997) 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 Abbot2012). 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 Kienert2014) 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 Abbot2012; 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 S0=1361 W m−2 (Kopp and Lean2011) 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 CO2 concentrations bracketing the Snowball bifurcation (see Table 1). In addition, we have run model experiments at two fixed levels of CO2 (10 000 and 10 ppm) and decreasing solar luminosities of 1140, 1130, 1125, and 1120 W m−2 for 10 000 ppm, and 1334, 1329, and 1324 W m−2 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 CO2 concentrations. All simulations were integrated for at least 2000 model years after the last change in CO2 concentration to allow the ocean to approach equilibrium conditions.

Table 1Overview 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/S0 to the present-day value S0, and the atmospheric CO2 concentration of not-fully-glaciated and Snowball states.

Download Print Version | Download XLSX

3 Results

3.1 The Snowball bifurcation for aquaplanets through time

The critical CO2 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 CO2 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 CO2 concentration in Earth's atmosphere and to earlier modelling studies for individual time slices, the latter being differentiated by model type.

Figure 1Snowball bifurcation (in terms of CO2 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 CO2. 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 CO2 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 CO2 concentration calculated from the combined forcing. All model results are scaled to a modern solar constant of 1361 W m−2 (Kopp and Lean2011). Estimates of past atmospheric CO2 levels are shown in grey; Phanerozoic CO2 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 CO2 concentration in our aquaplanet simulations falls from ∼105 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 CO2 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 CO2 radiative forcing (Huang and Bani Shahabadi2014). The downturn of the critical CO2 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 -15C of the critical states (see Fig. 2), there is greenhouse warming due to water vapour, which becomes significant at very low CO2 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 H2O rotational bands because of the lower temperature. (3) Finally, there is also a warming contribution from cloud radiative effects.

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

(1) p CO 2 , crit = a 1 exp a 2 a 3 - S S 0 a 4 .

Fitting this function to the points derived from the aquaplanet simulations yields the following parameters: a1=(0.0836±0.0721) ppm, a2=26.6±0.7, a3=1.0±0.00002, and a4=0.475±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 CO2, 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 (Feulner2011). Similarly, our simulations with Neoproterozoic continents (Feulner and Kienert2014) indicate slightly higher critical CO2 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 CO2 concentrations of the aquaplanet simulations for modern boundary conditions (Yang et al.2012a, b, c) and the Late Paleozoic ice age (Feulner2017). 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 CO2 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 Jacob2004; Voigt and Abbot2012; Rose2015).

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 Abbot2012).

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 CO2 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 Marotzke2010), 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 CO2 and a fixed present-day solar constant, the bifurcation point in our aquaplanet simulations is below a CO2 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 CO2 concentration of 1 ppm. Voigt and Marotzke (2010) used ECHAM5/MPI-OM, finding a Snowball state at 0.1 ppm of CO2 for present-day continents and a slightly higher value of the solar constant of 1367 W m−2, 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 CO2 concentrations from below 40 ppm in an AGCM simulation (Chandler and Sohl2000) to about 100–700 ppm in AOGCMs and OGCMs with sea-ice dynamics (Poulsen et al.2002; Voigt et al.2011; Voigt and Abbot2012; Liu et al.2013; Feulner and Kienert2014; 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 CO2 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 (Feulner2012; 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 Toon2013; Kunze et al.2014; Le Hir et al.2014; Teitler et al.2014). The critical CO2 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 CO2 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 CO2 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 -2C and average sea-ice fractions of about 33 %, whereas critical states at higher solar luminosities exhibit much lower temperatures of about -15C 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 3Global and annual mean sea-ice fraction in equilibrium climate states at 900 Ma for various concentrations of atmospheric CO2. 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 CO2 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 CO2 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 fsea ice. In agreement with the findings of Fig. 2, these states (defined by 0.497<fseaice<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 4Transient 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 0.497<fseaice<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 5Maps 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−1 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 6Zonal 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 7Schematic 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 φ:

(2) Q ( φ ) ( 1 - α ( φ ) ) = Q abs ( φ ) = σ ε ( φ ) T sfc 4 ( φ ) + div F ( φ ) .

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, Qabs 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), Tsfc 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 Tsfc,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 ΔTsfc between a given equilibrium state and a reference state, denoted by subscript 0, can then be calculated as follows:


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 CO2 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 CO2 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 CO2 concentration decreased more or less logarithmically (from ∼105 ppm 4 billion years ago to ∼102 ppm 700 million years ago), as one would expect from the logarithmic character of the CO2 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 CO2 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 -2C 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 -15C 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 CO2 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 CO2 concentrations or at the CO2 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 CO2 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 CO2 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 CO2 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 CO2 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. Schneider2006). 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.

Code and data availability

Model input and output files, as well as the scripts used to generate the figures, are available at the institutional repository of the Potsdam Institute for Climate Impact Research (, Feulner et al.2022). The model source code can be made available upon request.

Author contributions

GF designed the study. MB and GF performed and analysed the model simulations. SP provided technical assistance and compiled the data archive. GF prepared all figures and wrote the paper with input from all the co-authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors would like to thank Julius Eberhard, Alexey V. Eliseev, and Anna Feulner for help and discussions. We also thank the reviewers, Aiko Voigt and Yonggang Liu, for their constructive feedback, which helped to improve the paper. The European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg are gratefully acknowledged for supporting this project by providing resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research.

Financial support

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

Review statement

This paper was edited by Laurens Ganzeveld and reviewed by Yonggang Liu and Aiko Voigt.


Abbot, D. S., Voigt, A., and Koll, D.: The Jormungand global climate state and implications for Neoproterozoic glaciations, J. Geophys. Res., 116, D18103,, 2011. a

Abe, Y., Abe-Ouchi, A., Sleep, N. H., and Zahnle, K. J.: Habitable Zone Limits for Dry Planets, Astrobiology, 11, 443–460,, 2011. a

Bendtsen, J.: Climate sensitivity to changes in solar insolation in a simple coupled climate model, Clim. Dynam., 18, 595–609,, 2002. a

Bonanno, A., Schlattl, H., and Paternò, L.: The age of the Sun and the relativistic corrections in the EOS, Astron. Astrophys., 390, 1115–1118,, 2002. a

Braun, C., Hörner, J., Voigt, A., and Pinto, J. G.: Ice-free tropical waterbelt for Snowball Earth events questioned by uncertain clouds, Nat. Geosci., 15, 489–493,, 2022. a, b, c, d

Budyko, M. I.: The effect of solar radiation variations on the climate of the earth, Tellus, 21, 611–619,, 1969. a, b

Cahalan, R. F. and North, G. R.: A Stability Theorem for Energy-Balance Climate Models, J. Atmos. Sci., 36, 1178–1188,<1178:ASTFEB>2.0.CO;2, 1979. a

Chandler, M. A. and Sohl, L. E.: Climate forcings and the initiation of low-latitude ice sheets during the Neoproterozoic Varanger glacial interval, J. Geophys. Res., 105, 20737–20756,, 2000. a, b, c, d

Charnay, B., Forget, F., Wordsworth, R., Leconte, J., Millour, E., Codron, F., and Spiga, A.: Exploring the faint young Sun problem and the possible climates of the Archean Earth with a 3-D GCM, J. Geophys. Res., 118, 10414–10431,, 2013. a, b, c

Charnay, B., Wolf, E. T., Marty, B., and Forget, F.: Is the Faint Young Sun Problem for Earth Solved?, Space Sci. Rev., 216, 90,, 2020. a, b, c

Claussen, M., Mysak, L. A., Weaver, A. J., Crucifix, M., Fichefet, T., Loutre, M.-F., Weber, S. L., Alcamo, J., Alexeev, V. A., Berger, A., Calov, R., Ganopolski, A., Goosse, H., Lohmann, G., Lunkeit, F., Mokhov, I. I., Petoukhov, V., Stone, P., and Wang, Z.: Earth system models of intermediate complexity: closing the gap in the spectrum of climate system models, Clim. Dynam., 18, 579–586,, 2002. a

Donnadieu, Y., Ramstein, G., Fluteau, F., Roche, D., and Ganopolski, A.: The impact of atmospheric and oceanic heat transports on the sea-ice-albedo instability during the Neoproterozoic, Clim. Dynam., 22, 293–306,, 2004. a

Drazin, P. G. and Griffel, D. H.: On the Branching Structure of Diffusive Climatological Models., J. Atmos. Sci., 34, 1696–1706,<1696:OTBSOD>2.0.CO;2, 1977. a

Driese, S. G., Jirsa, M. A., Ren, M., Brantley, S. L., Sheldon, N. D., Parker, D., and Schmitz, M.: Neoarchean paleoweathering of tonalite and metabasalt: Implications for reconstructions of 2.69 Ga early terrestrial ecosystems and paleoatmospheric chemistry, Precambrian Res., 189, 1–17,, 2011. a

Faegre, A.: An Intransitive Model of the Earth-Atmosphere-Ocean System, J. Appl. Meteorol., 11, 4–6,<0004:AIMOTE>2.0.CO;2, 1972. a

Feulner, G.: Are the most recent estimates for Maunder Minimum solar irradiance in agreement with temperature reconstructions?, Geophys. Res. Lett., 38, L16706,, 2011. a

Feulner, G.: The Faint Young Sun Problem, Rev. Geophys., 50, RG2006,, 2012. a, b, c

Feulner, G.: Formation of most of our coal brought Earth close to global glaciation, P. Natl. Acad. Sci., 114, 11333–11337,, 2017. a, b, c

Feulner, G. and Kienert, H.: Climate simulations of Neoproterozoic snowball Earth events: Similar critical carbon dioxide levels for the Sturtian and Marinoan glaciations, Earth Planet. Sc. Lett., 404, 200–205,, 2014. a, b, c, d, e, f

Feulner, G., Hallmann, C., and Kienert, H.: Snowball cooling after algal rise, Nature Geosci., 8, 659–662,, 2015. a

Feulner, G., Bukenberger, M. S., and Petri, S.: Simulation data for tracing snowball bifurcation on an earth-like aquaplanet over 4 billion years, Scientific Data Set, V. 2022-07-27, GFZ Data Services [data set],, 2022. a

Fichefet, T. and Morales Maqueda, M. A.: Sensitivity of a global sea ice model to the treatment of ice thermodynamics and dynamics, J. Geophys. Res.-Oceans, 102, 12609–12646,, 1997. a

Fiorella, R. P. and Sheldon, N. D.: Equable end Mesoproterozoic climate in the absence of high CO2, Geology, 45, 231–234,, 2017. a, b

Foster, G. L., Royer, D. L., and Lunt, D. J.: Future climate forcing potentially without precedent in the last 420 million years, Nat. Commun., 8, 14845,, 2017. a

Frierson, D. M. W., Lu, J., and Chen, G.: Width of the Hadley cell in simple and comprehensive general circulation models, Geophys. Res. Lett., 34, L18804,, 2007. a

Gal-Chen, T. and Schneider, S. H.: Energy balance climate modeling: Comparison of radiative and dynamic feedback mechanisms, Tellus, 28, 108–121,, 1976. a

Ghil, M.: Climate Stability for a Sellers-Type Model, J. Atmos. Sci., 33, 3–20,<0003:CSFAST>2.0.CO;2, 1976. a

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

Heinemann, M., Jungclaus, J. H., and Marotzke, J.: Warm Paleocene/Eocene climate as simulated in ECHAM5/MPI-OM, Clim. Past, 5, 785–802,, 2009. a

Held, I. M. and Suarez, M. J.: Simple albedo feedback models of the icecaps, Tellus, 26, 613–629,, 1974. a

Hessler, A. M., Lowe, D. R., Jones, R. L., and Bird, D. K.: A lower limit for atmospheric carbon dioxide levels 3.2 billion years ago, Nature, 428, 736–738,, 2004. a

Hoffman, P. F. and Schrag, D. P.: The snowball Earth hypothesis: testing the limits of global change, Terra Nova, 14, 129–155,, 2002. a, b

Hofmann, M. and Morales Maqueda, M. A.: Performance of a second-order moments advection scheme in an Ocean General Circulation Model, J. Geophys. Res., 111, C05006,, 2006. a

Hörner, J., Voigt, A., and Braun, C.: Snowball Earth initiation and the thermodynamics of sea ice, J. Adv. Model. Earth Syst., 14, e2021MS002734,, 2022. a, b, c

Huang, J. and Bowman, K. P.: The small ice cap instability in seasonal energy balance models, Clim. Dynam., 7, 205–215,, 1992. a

Huang, Y. and Bani Shahabadi, M.: Why logarithmic? A note on the dependence of radiative forcing on gas concentration, J. Geophys. Res.-Atmos., 119, 13683–13689,, 2014. a

Ikeda, T. and Tajika, E.: A study of the energy balance climate model with CO2-dependent outgoing radiation: implication for the glaciation during the Cenozoic, Geophys. Res. Lett., 26, 349–352,, 1999. a

Kanzaki, Y. and Murakami, T.: Estimates of atmospheric CO2 in the Neoarchean-Paleoproterozoic from paleosols, Geochim. Cosmochim. Acta, 159, 190–219,, 2015. a

Kasting, J. F.: Theoretical Constraints on Oxygen and Carbon Dioxide Concentrations in the Precambrian Atmosphere, Precambrian Res., 34, 205–229,, 1987. a

Kienert, H., Feulner, G., and Petoukhov, V.: Faint young Sun problem more severe due to ice-albedo feedback and higher rotation rate of the early Earth, Geophys. Res. Lett., 39, L23710,, 2012. a, b, c, d

Kienert, H., Feulner, G., and Petoukhov, V.: Albedo and heat transport in 3-D model simulations of the early Archean climate, Clim. Past, 9, 1841–1862,, 2013. a

Kopp, G. and Lean, J. L.: A new, lower value of total solar irradiance: Evidence and climate significance, Geophys. Res. Lett., 38, L01706,, 2011. a, b

Kunze, M., Godolt, M., Langematz, U., Grenfell, J. L., Hamann-Reinus, A., and Rauer, H.: Investigating the Early Earth Faint Young Sun Problem with a General Circulation Model, Planet. Space Sci., 98, 77–92,, 2014. a, b, c, d

Le Hir, G., Teitler, Y., Fluteau, F., Donnadieu, Y., and Philippot, P.: The faint young Sun problem revisited with a 3-D climate–carbon model – Part 1, Clim. Past, 10, 697–713,, 2014. a, b, c

Lehmer, O. R., Catling, D. C., Buick, R., Brownlee, D. E., and Newport, S.: Atmospheric CO2 levels from 2.7 billion years ago inferred from micrometeorite oxidation, Sci. Adv., 6, eaay4644,, 2020. a

Lewis, J. P., Weaver, A. J., Johnston, S. T., and Eby, M.: Neoproterozoic “snowball Earth”: Dynamic sea ice over a quiescent ocean, Paleoceanography, 18, 1092,, 2003. a, b, c, d, e

Lewis, J. P., Weaver, A. J., and Eby, M.: Snowball versus slushball Earth: Dynamic versus nondynamic sea ice?, J. Geophys. Res., 112, C11014,, 2007. a, b, c

Lindzen, R. S. and Farrell, B.: Some Realistic Modifications of Simple Climate Models, J. Atmos. Sci., 34, 1487–1501,<1487:SRMOSC>2.0.CO;2, 1977. a

Liu, Y., Peltier, W. R., Yang, J., and Vettoretti, G.: The initiation of Neoproterozoic “snowball” climates in CCSM3: the influence of paleocontinental configuration, Clim. Past, 9, 2555–2577,, 2013. a, b, c, d, e, f, g

Liu, Y., Peltier, W. R., Yang, J., Vettoretti, G., and Wang, Y.: Strong effects of tropical ice-sheet coverage and thickness on the hard snowball Earth bifurcation point, Clim. Dyn., 48, 3459–3474,, 2017. a, b, c, d

Montoya, M., Griesel, A., Levermann, A., Mignot, J., Hofmann, M., Ganopolski, A., and Rahmstorf, S.: The earth system model of intermediate complexity CLIMBER-3α. Part 1: description and performance for present-day conditions, Clim. Dyn., 25, 237–263,, 2005. a

North, G. R.: Analytical Solution to a Simple Climate Model with Diffusive Heat Transport., J. Atmos. Sci., 32, 1301–1307,<1301:ASTASC>2.0.CO;2, 1975a. a

North, G. R.: Theory of Energy-Balance Climate Models, J. Atmos. Sci., 32, 2033–2043,<2033:TOEBCM>2.0.CO;2, 1975b. a

North, G. R.: Multiple solutions in energy balance climate models, Global Planet. Change, 2, 225–235,, 1990. a, b

North, G. R. and Coakley, James A., J.: Differences between Seasonal and Mean Annual Energy Balance Model Calculations of Climate and Climate Sensitivity., J. Atmos. Sci., 36, 1189–1204,<1189:DBSAMA>2.0.CO;2, 1979. a

North, G. R., Cahalan, R. F., and Coakley, James A., J.: Energy Balance Climate Models, Rev. Geophys. Space Phys., 19, 91,, 1981. a, b, c

North, G. R., Short, D. A., and Mengel, J. G.: Simple energy balance model resolving the seasons and the continents – Application to the astronomical theory of the ice ages, J. Geophys. Res., 88, 6576–6586,, 1983. a

Öpik, E. J.: On the causes of palaeoclimatic variations and of the ice ages in particular, J. Glaciol., 2, 213–218,, 1953. a

Pacanowski, R. C. and Griffies, S. M.: The MOM-3 manual, Tech. Rep. 4, NOAA/Geophyical Fluid Dynamics Laboratory, Princeton, NJ, USA, (last access: 29 April 2023), 1999. a

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

Pierrehumbert, R. T., Abbot, D. S., Voigt, A., and Koll, D.: Climate of the Neoproterozoic, Annu. Rev. Earth Pl. Sc, 39, 417–460,, 2011. a, b, c, d

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

Poulsen, C. J., Jacob, R. L., Pierrehumbert, R. T., and Huynh, T. T.: Testing paleogeographic controls on a Neoproterozoic snowball Earth, Geophys. Res. Lett., 29, 1515,, 2002. a, b, c, d

Ramanathan, V. and Coakley, J. A., J.: Climate Modeling Through Radiative-Convective Models (Paper 8R0533), Rev. Geophys. Space Phys., 16, 465,, 1978. a

Roe, G. H. and Baker, M. B.: Notes on a Catastrophe: A Feedback Analysis of Snowball Earth, J. Climate, 23, 4694–4703,, 2010. a

Romanova, V., Lohmann, G., and Grosfeld, K.: Effect of land albedo, CO2, orography, and oceanic heat transport on extreme climates, Clim. Past, 2, 31–42,, 2006. a, b, c

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

Rose, B. E. J. and Marshall, J.: Ocean Heat Transport, Sea Ice, and Multiple Climate States: Insights from Energy Balance Models, J. Atmos. Sci., 66, 2828,, 2009. a

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

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

Schneider, S. H. and Gal-Chen, T.: Numerical experiments in climate stability, J. Geophys. Res., 78, 6182–6194,, 1973. a

Schneider, T.: The General Circulation of the Atmosphere, Annu. Rev. Earth Planet. Sci., 34, 655–688,, 2006.  a

Sellers, W. D.: A Global Climatic Model Based on the Energy Balance of the Earth-Atmosphere System, J. Appl. Meteorol., 8, 392–400,<0392:AGCMBO>2.0.CO;2, 1969. a, b

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

Shen, S. and North, G.: A simple proof of the slope stability theorem for energy balance climate models, Canadian Applied Mathematics Quarterly, 7, 203–215, 1999. a

Strauss, J. V. and Tosca, N. J.: Mineralogical constraints on Neoproterozoic pCO2 and marine carbonate chemistry, Geology, 48, 599–603,, 2020. a

Tang, H. and Chen, Y.: Global glaciations and atmospheric change at ca. 2.3 Ga, Geosci. Front., 4, 583–596,, 2013. a, b

Teitler, Y., Le Hir, G., Fluteau, F., Philippot, P., and Donnadieu, Y.: Investigating the Paleoproterozoic glaciations with 3-D climate modeling, Earth Planet. Sci. Lett., 395, 71–80,, 2014. a, b

Voigt, A. and Abbot, D. S.: Sea-ice dynamics strongly promote Snowball Earth initiation and destabilize tropical sea-ice margins, Clim. Past, 8, 2079–2092,, 2012. a, b, c, d, e, f, g, h, i

Voigt, A. and Marotzke, J.: The transition from the present-day climate to a modern Snowball Earth, Clim. Dynam., 35, 887–905,, 2010. a, b, c, d

Voigt, A., Abbot, D. S., Pierrehumbert, R. T., and Marotzke, J.: Initiation of a Marinoan Snowball Earth in a state-of-the-art atmosphere-ocean general circulation model, Clim. Past, 7, 249–263,, 2011. a, b, c, d, e

von Paris, P., Rauer, H., Lee Grenfell, J., Patzer, B., Hedelt, P., Stracke, B., Trautmann, T., and Schreier, F.: Warming the early earth – CO2 reconsidered, Planet. Space Sci., 56, 1244–1259,, 2008. a

Wolf, E. T. and Toon, O. B.: Hospitable Archean Climates Simulated by a General Circulation Model, Astrobiology, 13, 656–673,, 2013. a, b, c

Yang, J., Peltier, W. R., and Hu, Y.: The Initiation of Modern “Soft Snowball” and “Hard Snowball” Climates in CCSM3. Part I: The Influences of Solar Luminosity, CO2 Concentration, and the Sea Ice/Snow Albedo Parameterization, J. Climate, 25, 2711–2736,, 2012a. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Yang, J., Peltier, W. R., and Hu, Y.: The Initiation of Modern “Soft Snowball“ and “Hard Snowball” Climates in CCSM3. Part II: Climate Dynamic Feedbacks, J. Climate, 25, 2737–2754,, 2012b. a, b, c

Yang, J., Peltier, W. R., and Hu, Y.: The initiation of modern soft and hard Snowball Earth climates in CCSM4, Clim. Past, 8, 907–918,, 2012c. a, b, c, d, e, f, g, h, i, j, k

Short summary
One limit of planetary habitability is defined by the threshold of global glaciation. If Earth cools, growing ice cover makes it brighter, leading to further cooling, since more sunlight is reflected, eventually leading to global ice cover (Snowball Earth). We study how much carbon dioxide is needed to prevent global glaciation in Earth's history given the slow increase in the Sun's brightness. We find an unexpected change in the characteristics of climate states close to the Snowball limit.
Final-revised paper