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 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/S_{\mathrm{0}}$ to the present-day value S₀, and the atmospheric CO₂ concentration of not-fully-glaciated and Snowball states.

Download Print Version | Download XLSX

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 1Snowball 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).

Download

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):

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

Download

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 3Global 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.

Download

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 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 $\mathrm{0.497}<f_{\mathrm{sea}\mathrm{ice}}<\mathrm{0.503}$ to define the transient Hadley phase, indicated by the horizontal lines.

Download

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⁻¹ indicated in the bottom right corner.

Download

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

Download

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.

Download

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 φ:

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:

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

Download

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.