2.2.1 Sampling continental configurations

For each of the two Neoproterozoic snowball-Earth onsets – the Sturtian onset (taken here as 720 Ma) and the Marinoan (here 650 Ma) – four different continental configurations are applied as boundary conditions: (1) two reconstructions and (2) two idealised configurations. Figure 2 shows the distribution of continents and the bathymetry in the different configurations along with the zonally averaged land fraction. The idea behind this sampling is (1) to obtain an estimate of the uncertainty induced by using different realistic palaeogeographies and (2) to roughly cover the largest possible uncertainty of global-mean temperatures caused by the distribution of continents.

Figure 2Distributions of the exposed continental surface, ocean bathymetry, and zonally averaged land fraction for the different continental configurations. Age specifications (720 and 635 Ma) refer to the time of the reconstructions, L+13 denotes reconstructions based on Li et al. (2013), and M+21 denotes those based on Merdith et al. (2021). Land fractions are shown as simple zonal averages (hatched) and zonal averages weighted by the actual area of the latitudinal bands (grey).

The reconstructions are based on published datasets from Li et al. (2013) and Merdith et al. (2021) dated to 720 Ma for the Sturtian and 635 Ma for the Marinoan snowball onset. (Note, however, that we perform Marinoan simulations with the solar constant at 650 Ma.) These reconstructions rely upon narrowly constrained palaeomagnetic datasets derived from cratonic lithologies as well as plate kinematics and tectonostratigraphic correlations. The more refined methods applied in the Merdith et al. (2021) palaeogeography also allow the trajectories of large terranes to be projected through geologic history, producing a more robust model foundation during a time of hypothesised metamorphism following the breakup of the supercontinent Rodinia (Merdith et al., 2017). For the remainder of the paper, the reconstructions based on Li et al. (2013) will be referred to as L+13, and those based on Merdith et al. (2021) will be referred to as M+21.

The land–ocean masks based on the reconstructions are re-gridded to the ocean-model grid with a horizontal resolution of $\mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$. Some grid cells have to be manually adjusted, for example to avoid isolated ocean cells within continents or narrow ocean channels. For the L+13 reconstruction, continental shelves and slopes with a typical width of one ocean grid cell each are assigned by hand around all land areas (Feulner and Kienert, 2014). For the M+21 reconstruction and the idealised configurations, we use an automatic procedure with shelves and slopes with a width of 500 km each. Ocean depth levels are randomly assigned and fall within the range 100–265 m for continental shelves, 349–993 m for continental slopes, and 1243–4987 m for the abyss. Land cells are set to either a constant elevation of 200 m or randomly distributed elevations between 200 and 400 m. Land-surface type is assumed to be bare soil unless there is snow cover.

The idealised continental configurations are constructed by shifting continental plates of the 635 Ma reconstruction by Merdith et al. (2021) using the GPlates program (Müller et al., 2018) and by subsequently performing the manual adjustments described above. We construct a “dispersed” configuration of small continents distributed evenly over the surface and a “polar supercontinent” configuration, where the single continent is the only landmass and is centred at the South Pole of the model grid. Depth levels of continental shelves, slopes, and deep-ocean cells are assigned in the same manner as for the M+21 reconstruction described above.

The use of different sources and the manual adjustments come with variations in the exposed continental-surface area. This is lower for the reconstructions, with the 720 Ma land fraction (L+13: 16.8 %, M+21: 17.1 %) being smaller than that at 635 Ma (L+13: 17.4 %, M+21: 20.0 %). The idealised dispersed configuration has a land fraction of 20.8 %, and that for the supercontinent is 24.1 %.

Table 1Overview of the simulations regarding the sensitivity to continental distribution, differentiating between those aimed at the temperature sensitivity at fixed CO₂ concentration and those aimed at the CO₂ threshold for snowball inception. All simulations are initialised in a warmer climate than the shown equilibrium states and are run with an orbital geometry of 23.5^∘ obliquity and zero eccentricity. S₀ is the solar constant. “L+13” denotes reconstructions based on Li et al. (2013); “M+21” denotes those based on Merdith et al. (2021).

Download Print Version | Download XLSX

Table 1 provides an overview of the simulations with differing continental configurations. We use the solar constant S₀=1283 W m⁻² for the Sturtian and 1290 W m⁻² for the Marinoan onset. Both values are derived from a standard solar model (Bahcall et al., 2001) using a present-day value of 1361 W m⁻². For the orbital geometry, we select an obliquity of 23.5^∘ and zero eccentricity. All unperturbed runs are initialised in a warm and ice- and snow-free state with a global-mean surface-air temperature of around 12 ^∘C.

In the simulations aiming at a comparison of the energy-balance effects from differing continental distributions, we run the model with an atmospheric CO₂ concentration of 250 parts per million by volume (ppm) for the Sturtian and 200 ppm for the Marinoan conditions. Non-CO₂ greenhouse gases are not considered in our simulations. Scanning for the CO₂ threshold is done by letting the model approach equilibrium for various CO₂ concentrations in steps of 10 ppm.

2.2.2 Sampling orbital geometries

In CLIMBER-3α, the incoming top-of-atmosphere solar radiation is determined by three orbital parameters: obliquity ε, eccentricity e, and argument of perihelion ω. $\mathit{\omega }=\mathrm{0}{}^{\circ }$ corresponds to perihelion at the September (southward) equinox and $\mathit{\omega }=\mathrm{90}{}^{\circ }$ to that at the December (southern) solstice. We sample the ε–e–ω parameter space by selecting the coordinates ε= (22, 23.5, 24.5^∘), e=(0, 0.03, 0.069), and ω= (0, 45, 90, 135, 180, 225, 270, 315^∘) and running the model for selected physically distinguishable combinations thereof. (Note that for e=0, the different values of ω are indistinguishable in terms of the annual and seasonal distributions of insolation.) The ranges of ε and e are confined by computations of Earth's orbital variations during the past 250 million years (Laskar et al., 2004). The climate in these simulations is initialised in the same way as for the set of simulations described above.

An overview of the simulations with varying orbital geometries is given in Table 2. For all 51 distinguishable configurations, we run the model at a CO₂ concentration of 150 ppm in order to compare the effects of orbital geometry on the surface-air temperature. Since these effects are significantly smaller than those of the continental configurations, we perform finer scanning for the CO₂ thresholds at nine different orbital geometries: after letting the model approach a steady state, which usually takes around 1000–3000 years, we reduce CO₂ by 1 ppm – or more if it is well above the threshold – and repeat this procedure until we have a snowball inception in each simulation. The continents and solar constant always correspond to the Sturtian onset (720 Ma) using the respective L+13 continental reconstruction.

Table 2Overview of the simulations regarding the sensitivity to orbital geometry, differentiating between those aimed at the temperature sensitivity at fixed CO₂ concentration and those aimed at the CO₂ threshold for snowball inception. All simulations are initialised in a warmer climate than the shown equilibrium states, run at a solar constant of 1283 W m⁻², and use the Sturtian L+13 palaeogeography. ε is the obliquity, e is the eccentricity, and ω is the argument of perihelion.

Download Print Version | Download XLSX

2.2.3 Sampling single global volcanic perturbations

We implement globally averaged volcanic perturbations as annually resolved transient variations in the solar constant. The radiative perturbation resulting from a single explosive eruption is idealised as follows: in the calendar year after the eruption, the perturbation reaches a peak forcing, which is followed by an exponential decrease ($\sim {\mathrm{e}}^{-\mathit{\lambda }{t}_i}$, where the decay rate λ has a value of 1 yr⁻¹) in successive calendar years t_i. In the calendar year of the eruption, the forcing is half the peak forcing, representing a linear buildup due to sulfate aerosols being injected into the stratosphere. This form of the perturbation, including the value of λ, is commonly used in the literature (e.g. Ammann et al., 2003; Gao et al., 2008; Metzner et al., 2014).

With the aim of sampling the peak forcing realistically, we construct our idealised scenarios by following ice-core-based reconstructions of the net solar radiative forcing ΔS_net(t) due to volcanism during the years 850–2000 common era (CE; Schmidt et al., 2012, based on Gao et al., 2008, and Crowley et al., 2008) as well as model-based forcings of the ∼ 74 000 a Toba (Timmreck et al., 2010) and 2.1 Ma Yellowstone eruptions (Segschneider et al., 2013). The annual- and global-mean peak forcings from these reconstructions cover a range of up to −20 W m⁻². We define our single-eruption scenarios using a nominal reduction of the global-mean net solar radiation, here called “nominal volcanic forcing” (NVF), with peaks of −2, −5, −10, −20, −30, and −40 W m⁻², thus also allowing for what we consider extreme and very rare perturbations beyond the reconstructions. Note that some previously applied (Jones et al., 2005; Robock et al., 2009; Gupta et al., 2019) estimates of the ∼ 74 000 a Toba eruption exceed a peak forcing of −100 W m⁻² but have been criticised on physical grounds as being much too strong (Timmreck et al., 2010, 2012).

The translation of a given NVF to a change in the solar constant results from a basic consideration of the top-of-atmosphere short-wave radiation balance: the daily- and zonal-mean net solar radiation is $S_{\mathrm{net}}\left(\mathit{\varphi }\right)=A\left(\mathit{\varphi }\right)g(\mathit{\varphi },{\mathit{\varphi }}_{\mathrm{z}})S_{\mathrm{0}}$, where ϕ is the latitude, A is the planetary coalbedo (1−albedo), g is a factor comprising geometric effects of Earth and its orbit, and ϕ_z is the latitude at which the Sun reaches the zenith at noon (e.g. Peixoto and Oort, 1992, their Sect. 6.3.2). For the annual and global mean, denoted by an overline, we approximate $\stackrel{\mathrm{‾}}{A\left(\mathit{\varphi }\right)g(\mathit{\varphi },{\mathit{\varphi }}_{\mathrm{z}})}$ by $A_{\mathrm{eff}}\stackrel{\mathrm{‾}}{g}$, where A_eff is a geometry- and radiation-weighted effective global coalbedo and $\stackrel{\mathrm{‾}}{g}=(\mathrm{4}\sqrt{\mathrm{1}-e^{\mathrm{2}}}{)}^{-\mathrm{1}}$ (Berger and Loutre, 1994). We use the present-day global coalbedo of 0.7 (Goode et al., 2001) as an estimate for A_eff. Thus, we translate the NVF into the solar-constant forcing ΔS₀ via

where t_i represents successive calendar years.

Three comments on the construction of the volcanic forcing are appropriate at this point. First, the volcanic forcing is called “nominal” here because, in this construction, each eruption corresponds to a reduction in the solar constant of around 5.7 NVF for the used eccentricity range, independent of the actual climate state. This means that, while the NVF represents the volcanic effects in a Quaternary or modern climate, the actual reduction in the net solar flux is somewhat smaller in a colder climate with a lower coalbedo. For the sake of simplicity and comparability with other studies, however, we continue to specify the eruption scenarios using the NVF.

Second, the use of annual reductions in the solar constant implies that volcanic eruptions always start at the beginning of a calendar year, i.e. in northern-hemispheric winter.

Third, Macdonald and Wordsworth (2017) argued on the basis of a vertical atmospheric and volcanic-plume model that volcanic eruptions can have a particularly strong impact on snowball inceptions because the tropopause, the lower boundary of the stratosphere, is generally lower in cooler climates. We do not explicitly include such an effect here.

All perturbed runs, including those described in the sections below, have Sturtian boundary conditions, making use of the respective L+13 palaeogeography. Only single global perturbations are tested for different orbital geometries, which are {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0}, {$\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$, e=0}, and {$\mathit{\epsilon }=\mathrm{24.5}{}^{\circ }$, e=0.069, $\mathit{\omega }=\mathrm{270}{}^{\circ }$}. See Table 3 for an overview of these simulations.

Table 3Overview of the simulations regarding the sensitivity of the CO₂ threshold for snowball inception to single global volcanic perturbations. All simulations are run at a solar constant of 1283 W m⁻² and use the L+13 palaeogeography for 720 Ma. ε is the obliquity, e is the eccentricity, and ω is the argument of perihelion. The “initial partially ice-covered attractor” refers to the conditions at which a particular simulation was initialised.

Download Print Version | Download XLSX

2.2.4 Sampling single zonal volcanic perturbations

To emulate the latitudinally resolved effects of large volcanic eruptions, we prescribe a transient zonal distribution of aerosol optical depth (AOD) yielding a proportional reduction in the solar constant at each latitude band. The aerosol time series are taken from Ammann et al. (2003), who estimated the monthly and zonally resolved AODs for historical eruptions between 1890 and 1999 CE based on a simple model of sulfate-aerosol transport in the stratosphere. We exemplarily select periods of 6–9 years including and following three eruptions at different latitudes: the 1902 Gagxanul/Santa María eruption (at 15^∘ N), the 1932 Quizapú/Cerro Azul eruption (36^∘ S), and the 1912 Novarupta eruption (58^∘ N). We interpolate these data to our atmospheric grid and prolong them to cover 50 years, assuming exponential decay at a rate of 1 yr⁻¹ after the last available values (Appendix Fig. A1). Additionally, we flip the data at the Equator to eventually get AOD time series for eruptions at 15, 36, and 58^∘ in both hemispheres.

In order to produce comparable results, we scale the AOD data for each peak forcing in such a way that the global-mean reduction in the solar constant integrated over 50 years is equal to the 50 year integral of ΔS₀(t) for the respective global volcanic perturbation described above. Regarding their effects on the model response, the only differences from the global perturbations are the finer spatial (latitudinal) and temporal (monthly) resolution. Note in particular that, as in the global case, eruptions always start in January.

All zonal perturbations are initialised on both the warmest and the coldest partially ice-covered attractors for the orbital geometry {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0}. See Table 4 for an overview.

Table 4Overview of the simulations regarding the sensitivity of the CO₂ threshold for snowball inception to single zonal volcanic perturbations. All simulations are run at a solar constant of 1283 W m⁻², use the Sturtian L+13 palaeogeography, and have an orbital geometry of 22^∘ obliquity and zero eccentricity.

Download Print Version | Download XLSX

2.2.5 Sampling sequences of global volcanic perturbations

We generate sequences of global volcanic perturbations along the lines of Ammann and Naveau (2010), who construct such annual-mean time series based on two processes, one purely stochastic and one purely deterministic.

First, the occurrence of explosive eruptions ejecting sulfate aerosols into the stratosphere is modelled as a stochastic annual time series b_p(t_i) consisting of binary values (“no eruption” or “eruption”, i.e. a Bernoulli process), where “eruption” has the constant probability p and “no eruption” has the probability 1−p. Each such eruption is then assigned a nominal volcanic peak forcing Π(t_i) randomly drawn from a generalised Pareto distribution with probability density

where β>0, ξ are real numbers, u and x are given in W m⁻², u > x, and the term in brackets is positive. This function asymptotically describes excesses over large thresholds of many probability distributions (Embrechts et al., 1997), including the peak radiative forcing of large volcanic eruptions (Naveau and Ammann, 2005). Note that in our case, the distribution describes negative events x beyond a negative threshold u. The nominal volcanic forcing, i.e. reduction of the net solar radiation, is translated into a reduction in the solar constant in the same way as described in Sect. 2.2.3. The complete first process, describing subsequent peaks Σ in the perturbation of the solar constant, therefore reads

where “d” over the equality sign implies equality in distribution.

Second, each peak is followed by an exponential decay as in our idealised single perturbations, which is described by the deterministic process $\mathrm{\Delta }{S}_{\mathrm{0}}\left(t_i\right)=\mathrm{\Delta }{S}_{\mathrm{0}}\left(t_{i-\mathrm{1}}\right)/\mathrm{e}$.

Both processes are combined, yielding the max-autoregressive model

Taking the maximum each year instead of the sum ensures that the effective distribution of peak forcing follows the distribution $g_{\mathit{\beta },\mathit{\xi },u}$, even for overlapping events.

The values of the parameters p, β, ξ, and u (characterising the frequency of eruptions and the distribution of peak forcings) define our scenarios. Roughly oriented towards the assumptions and parameter-fitting results by Ammann and Naveau (2010), we define the reference scenario as {p=0.1, β=1, ξ=0.3, u=1.5 W m⁻²}, which reproduces the volcanic-forcing statistics of the past ∼ 1000 years by and large. Based on this scenario, we subsequently change one parameter at a time in order to define five more scenarios. They comprise raising the eruption probability p to 0.3, varying the shape parameter β between 0.5 and 1.5 (a slightly wider range than in Ammann and Naveau, 2010), and shifting the threshold u to 0 and −3 W m⁻² (as in Ammann and Naveau, 2010). Tripling the probability of eruptions with a relevant radiative effect is an attempt to represent periods of more frequent volcanic activity, e.g. during the formation of large igneous provinces (LIPs) such as the Franklin LIP related to the Sturtian snowball inception by Macdonald and Wordsworth (2017). While u=0 W m⁻² should be only seen as a theoretical case because the generalised Pareto distribution is usually applied for nonzero thresholds (Naveau and Ammann, 2005), $u=-\mathrm{3}$ W m⁻² is chosen as the other limiting case since some effects can already be seen from a peak forcing of −2 W m⁻², as reported in Sect. 5.1. Changing ξ within or slightly beyond the range reported in Ammann and Naveau (2010) results in distributions very similar to some scenarios obtained from varying β and u. We therefore chose to leave it fixed at ξ=0.3 for all scenarios, yielding a Fréchet-type distribution with heavy tails (Embrechts et al., 1997).

Figure 3(a) Different forms of the generalised Pareto distribution determining the nominal peak forcing in the max-autoregressive model, which are displayed as unnormalised probability-density functions, i.e. they have a maximum value of 1. Dots denote the upper limits of each function's domain and range. For example, the function with $u=-\mathrm{3}$ W m⁻² allows for peak forcings at or below −3 W m⁻² only, while all functions are unbounded from below. (b) Nominal peak volcanic forcing in the simulations with each of the distributions from panel (a) compared to the reference scenario (diagonal line) at each time step of the generated volcanic time series. (c) Nominal volcanic forcing in the reference scenario. (d) Same as panel (c) but for the scenario with increased eruption probability.

Download

Figure 3a shows the five different parametrisations of $g_{\mathit{\beta },\mathit{\xi },u}$. In order to make the interpretation of results easier, we compute the pseudorandom sequences with a length of 5000 years for all different scenarios based on the same seed, i.e. initial algorithm state. Upon plotting the peak forcings of the reference scenario against the peak forcings from other scenarios for each time step (Fig. 3b), we can see that there are perfect correlations across the time series. Moreover, we can clearly characterise the effects of changing β and u: through the scaling effect of β, the scenarios yield dramatically different forcings for the rarest events, while u shifts the forcings homogeneously and thereby has the strongest effects on the most frequent eruptions, i.e. the “background noise”. Due to the construction, the scenario with an increased eruption probability p contains almost all of the individual events, including their magnitudes, from the reference case but adds many additional events (Fig. 3c, d). The background noise is strongly intensified, leading to a nearly permanent effect of volcanoes on the radiative balance. Additionally, this scenario contains a peak forcing (at around the year 2000 in Fig. 3d) that goes beyond the largest forcing in the reference scenario.

Table 5Overview of the simulations regarding the sensitivity of the CO₂ threshold for snowball inception to sequential volcanic perturbations. All simulations are run at a solar constant of 1283 W m⁻², use the Sturtian L+13 palaeogeography, have an orbital geometry of 22^∘ obliquity and zero eccentricity, and are initialised on the warmest attractor supported at the respective CO₂ concentrations.

Download Print Version | Download XLSX

All simulations are initialised on the warmest partially ice-covered attractors with an orbital geometry of {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0}. The reason for choosing only the warmest attractors is that in many simulations, the sequential perturbations force the model from the warmest to the coldest attractor. In these cases, the thresholds would become the same as for simulations initialised on the coldest attractors. In other words, the volcanic “noise” seems to connect the multiple stable states and push the model towards the coldest one. After applying the volcanic forcing for 5000 years, we let the model run without perturbations for 1000 years in order to ensure that we do not miss an eruption-induced transition to a snowball state. Table 5 summarises the simulations with sequential volcanic perturbations.

3.1.1 Dispersed configuration versus L+13 configuration

Let us now look at the individual energy-balance components (net solar radiation, meridional heat flux, and long-wave emission) in order to understand the large difference in surface temperature between the dispersed configuration and the Sturtian L+13 reconstruction (Fig. 4a(i)–(v)). The former case is significantly colder due to a great influence of the changed net solar radiation, i.e. albedo, and a minor effect of the increased global emissivity; see Fig. A3 for the global contributions. Changes in the meridional heat flux mostly act to compensate for zonal changes in the surface temperature caused by the other energy-balance components. When considering global averages, the meridional heat flux has virtually no effect on the surface temperature (Fig. A3), as anticipated in Sect. 2.3. We now proceed to discuss the components in more detail.

Since neither the solar constant nor the distribution of insolation is changed between the compared setups, differences in net solar radiation only represent albedo effects. These effects can already be detected in the spin-up of the simulations (not shown): after initialisation with nearly identical climatic conditions, the higher albedo and lower thermal inertia of high-latitude continents in the dispersed configuration lead to a strong cooling effect amplified by snowfall over land. As the simulations progress, the climate of the dispersed configuration consistently remains cooler than that of the L+13 reconstruction. Finally, upon ice margins encroaching into the mid-latitudes, the slightly larger continental area in the tropics further supports the cooler conditions in the dispersed case, favouring snowfall even there. At equilibrium, the most drastic differences in surface temperature occur over the sea-ice margins (Fig. 4a(i)). The cooling peaks at around 40^∘ S/N in the short-wave contribution are correlated with the differences in cryosphere fractions (Fig. 4a(iii)) and are likely the results of ice–albedo feedbacks at the ice margins. At the poles, the short-wave contribution yields higher temperatures in the dispersed configuration, where the presence of continents prevents the formation of sea ice and there is insufficient snowfall to compensate for the change in albedo.

Figure 5Long-term-mean sea-ice cover (50 % contours shown as thick black lines), surface ocean velocity (arrows with lengths proportional to the velocity), and heat flux at the ocean surface (colours; positive values indicate that the ocean gains heat) for Marinoan (650 Ma) equilibrium simulations using the dispersed (a), M+21 (b), and polar-supercontinent (c) configurations. The CO₂ concentration is 200 ppm (a), 170 ppm (b), and 20 ppm (c). Blue rectangles in panel (b) indicate ocean regions featuring local sea-ice protrusions.

Download

Changing the continental configuration can cause changes in the ocean circulation responsible for the meridional heat transport from lower to higher latitudes in the oceans (e.g. Jenkins, 1999), as also seen for modern-day analogies such as the role of the North Atlantic Current in the ocean temperatures of western Europe (Palter, 2015). Such changes can result in regional temperature effects which are difficult to detect in the zonal averages which are shown in Fig. 4. For the dispersed configuration, the ocean-surface currents are both diverted from the Equator and slowed down as compared to the open ocean in the supercontinent configuration (shown in Fig. 5a and c for the Marinoan conditions before snowball inception). The situation is similar for the heat transport between the Equator and the ice margins in the L+13 configuration (not shown).

Differences in emissivity, i.e. the greenhouse effect, within the modelled climate appear to have a localised influence that is dependent on the surface temperature and primarily due to the water-vapour feedback: the lower temperature in the dispersed configuration causes increased rates of evaporation and atmospheric concentrations of water vapour, again yielding lower temperatures. This demonstrates that emissivity is affected by the continental distribution indirectly and is a result of a feedback rather than an initial cause of cooling or warming. The influence of emissivity on the global climate can be investigated further by distinguishing the localised influences of cloud fraction and water vapour. A comparison of Fig. 4a(i), (iv), and (v) demonstrates that the water-vapour concentration strongly affects the long-wave emission, which is further modulated by cloud cover. Over the poles, the reduced cloud cover contributes to a cooling from reduced emissivity, which exceeds the warming from reduced ice cover.

A comparison of the dispersed configuration with the Marinoan L+13 reconstruction shows a very similar overall picture (Figs. A4 and A5). The general tendency for cooling in a dispersed configuration is corroborated by Hyde et al. (2000), who adopted boundary conditions including a SH supercontinent originally derived by Dalziel (1997). By increasing the degree of dispersion between the continents, Hyde et al. (2000) note that the threshold atmospheric CO₂ concentration for a snowball inception is increased.

3.1.2 M+21 configuration versus L+13 configuration

The Sturtian L+13 and M+21 reconstructions differ by 2 ^∘C regarding their global-mean surface temperatures, with the M+21 configuration being the colder one. This difference is mainly caused by SH albedo differences and amplified by emissivity changes (Figs. 4b(i)–(v) and A3).

The largest difference in zonally averaged land fraction between the two configurations is found in middle and higher latitudes of the Southern Hemisphere, where we also observe a pronounced increase in snow- and sea-ice cover. The cooler Southern Hemisphere is therefore likely a result of the increased continental area.

For the Marinoan simulations, the overall picture is very similar, but surface temperatures in the L+13 and M+21 reconstructions differ more strongly, mainly due to albedo differences (Fig. A5). Considering zonal averages (Fig. A4b(i)–(v)), the Marinoan configurations show larger local absolute differences in continental area and cryosphere fractions than the Sturtian equivalents. Additionally, the slightly greater low- and middle-latitude landmass area present in the Marinoan M+21 relative to the L+13 reconstruction reduces the absorption of solar radiation between the sea-ice margins.

For both the Sturtian and the Marinoan, the emissivity and meridional heat flux reflect the zonal cooling or warming caused by albedo differences. Nevertheless, the Marinoan M+21 configuration provides an example of local protrusions in the sea-ice front: due to the positioning of continental clusters, several middle-latitude ocean parts between 0 and 120^∘ E are secluded from ocean currents and, without an influx of warmer equatorial water, adopt a cooler regional climate (blue rectangles in Fig. 5).

3.1.3 Polar supercontinent configuration versus L+13 configuration

The SH polar supercontinent stands out as the warmest configuration (Figs. 4 and A3). Again, the largest part of the temperature difference in comparison to the Sturtian L+13 configuration results from net short-wave differences (Fig. A3c(i)–(v)). Almost the entire supercontinental domain has surface temperatures of several degrees, with the warming due only to albedo changes. The differences in cryosphere fractions (Fig. 4c(iii)) imply that it is the absence of sea ice which increases the albedo, despite the presence of snow on land.

Apart from the albedo, the difference in meridional heat transport seems to play an important part: the continental margin obstructs equatorial heat from bypassing a latitude of $\sim \mathrm{30}{}^{\circ }$ S via ocean currents. The resulting warming in this zone exceeds the cooling from snow cover. This enforces a decreased temperature reduction during the SH winter months, minimising sea-ice advancement between the Equator and 30^∘ S due to the stored thermal energy, as illustrated in Fig. 5c. Note that this configuration is not represented anywhere in geologic history, and so this mechanism likely does not elucidate any climate phenomenon currently under study.

Emissivity changes only partly reflect a simple water-vapour feedback. Over the supercontinent, cloud cover is homogeneously lower than in the L+13 configuration and gives rise to increased long-wave emission. This effect is compensated for by increased water-vapour concentrations only over the South Pole. Despite this cooling through long-wave radiation, most of the Southern Hemisphere is warmer in the supercontinent configuration, and the warmer conditions in the Northern Hemisphere (NH) can be partly explained by an ice–albedo feedback around the sea-ice margin. Perhaps another reason for the warmer NH climate is the absence of snow on land.

Very similar observations can be made for Marinoan conditions (Figs. A4 and A5).

4.1.1 Obliquity

As an example, we consider the zonal differences in the surface energy balance due to changing the orbital geometry from {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0} to {$\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$, e=0} (Fig. 8a(i)–(v)). For example, a negative temperature difference in this diagram implies a cooling contribution upon increasing the obliquity $\mathit{\epsilon }=\mathrm{22}{}^{\circ }\to \mathrm{23.5}{}^{\circ }$.

Figure 8Absolute difference in surface temperature ΔT (a(i)–c(i)), relative difference in top-of-atmosphere incoming solar flux $\mathit{\delta }{S}_{\mathrm{TOA}}^{\downarrow }=\mathrm{\Delta }\langle S_{\mathrm{TOA}}^{\downarrow }\rangle /\langle S_{\mathrm{TOA}}^{\downarrow }{\rangle }_{\mathrm{ref}}$ (a(ii)–c(ii)), absolute differences in the land-snow and sea-ice fraction Δf_cryo (a(iii)–c(iii)) and the cloud fraction Δf_cloud (a(iv)–c(iv)), and the relative difference in column water vapour (a(v)–c(v)) between long-term-averaged attractor states upon changing individual orbital parameters. Shown in (a(i))–(c(i)) are the total temperature difference and differences due to changes in net solar radiation (“short-wave”), emissivity (“long-wave”), and meridional heat-flux divergence (“meridional flux”). Panels (a(i)–(v)) show the effects of an increased obliquity ($\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$; the reference is 22^∘) with fixed e=0; panels (b(i)–(v)) show the effects of an increased eccentricity (e=0.069; the reference is 0.03) with fixed $\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$, $\mathit{\omega }=\mathrm{225}{}^{\circ }$; and panels (c(i)–(v)) show the effects of a shifted argument of perihelion ($\mathit{\omega }=\mathrm{90}{}^{\circ }$; the reference is 180^∘) with fixed $\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0.069. All differences are computed for zonal-mean magnitudes. The remaining boundary conditions are S₀=1283 W m⁻², a CO₂ concentration of 150 ppm, and the Sturtian L+13 palaeogeography. In the inset schematics, circles depict Earth and asterisks the Sun; new configurations (the minuends) are drawn with solid lines while reference configurations (the subtrahends) are drawn with dotted lines; schematics are not to scale.

Download

Since our obliquities are well below 55^∘, the increase in ε generally yields a more evenly distributed annual-mean incoming solar flux at the top of atmosphere (e.g. Linsenmeier et al., 2015, their Fig. 2), as also seen in our example (Fig. 8a(ii)). This can explain why the change in the net solar flux causes a surface cooling around the Equator and a surface warming around the poles. The remaining features of the short-wave contribution are correlated with the differences in cryosphere fractions (Fig. 8a(iii)) and seem to illustrate processes of local ice–albedo feedback at least. This feedback can also explain minor reductions seen in the sea-ice concentration and snow cover south of around 45^∘ S. Nevertheless, the global temperature differences are dominated by shifts of the cryosphere margins, which lie within the region of reduced insolation.

Changes in emissivity or heat-flux divergence mainly cause smaller differences in the surface temperature than the net solar flux. As an exception, the heat-flux divergence contributes to a cooling on the poleward sides of the cryosphere margins, exceeding the warming solar contribution. We interpret this as resulting from the decreased meridional temperature gradient in these regions: the local cooling around the cryosphere margins reduces that gradient polewards and thereby reduces the poleward heat transport there.

The local emissivity changes generally correlate well with the local balance between the short-wave and transport contributions (not shown separately, but this is evident from Fig. 8a(i)). We can attribute these changes to two positive feedbacks. First, a decrease in surface temperature lowers the atmospheric water-vapour concentration and thus raises the effective emissivity; this well-known water-vapour feedback occurs mainly in the tropics and mid-latitudes, but it also acts to enhance the high-latitude warming (Fig. 8a(v)). Second, the cooling around the cryosphere margins comes with a local decline in stratus clouds, which also increases the effective emissivity (Fig. 8a(iv)).

It is noteworthy that we observed a cooling with increased obliquity, while many studies report the opposite effect (e.g. Mantsis et al., 2011; Nowajewski et al., 2018; Brugger et al., 2019; Landwehrs et al., 2021). While those studies only investigated comparatively warm climate states with ice margins above 55^∘, the cooling effect shown here is consistent with the very cold late-Carboniferous simulations of Feulner (2017). In fact, for climate states with an ice line above 55^∘, our model simulates mostly rising temperatures with increased obliquity (Fig. A8).

All the discussed features are present in almost all other comparisons with differing obliquities but with other parameters fixed (not shown). There are two exceptions where increasing the obliquity from 23.5 to 24.5^∘ raises the temperature (e.g. see Fig. A6), but they are not discussed in detail here.

4.1.2 Eccentricity

We investigate the effects of increased eccentricity by comparing the orbital geometries {$\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$, e=0.03, $\mathit{\omega }=\mathrm{225}{}^{\circ }$} and {$\mathit{\epsilon }=\mathrm{23.5}{}^{\circ }$, e=0.069, $\mathit{\omega }=\mathrm{225}{}^{\circ }$} (Fig. 8b(i)–(v)) corresponding to perihelion in early May. In this example, we find globally increased surface temperatures as well as a hemispheric asymmetry regarding the total temperature differences.

Changes in the net solar flux contribute almost globally to the warming. A background contribution of well below +1 ^∘C can be explained by globally increased top-of-atmosphere insolation (Fig. 8b(ii)); this is due to the smaller semi-minor axis of the higher-eccentricity orbit, which results in a reduced mean Sun–Earth distance. Subannually, however, this distance oscillates, and the increased eccentricity yields a more pronounced perihelion and aphelion. For most values of ω, such as the one shown here, this increases the already existing net insolation difference between the hemispheres. This difference induces accelerated ice melting and strong warming during spring and summer in one hemisphere and delayed melting and weak relative cooling in the other. In principle, we see the same effect upon varying ω, which will be discussed below.

Outside the cryosphere margins, the warming contributions from the meridional heat-flux divergence are larger than those from the net solar flux. The local warming around the ice margins affects the meridional temperature gradients, the latter now being reduced around the Equator and enhanced over the polar ice caps. In the model, the correspondingly changing meridional heat transport causes the warming in both regions.

Changes in emissivity largely reflect the positive feedbacks already seen upon increasing the obliquity (Fig. 8b(iv)–(v)).

Depending on ω, the hemispherically asymmetric warming upon increasing the eccentricity from 0.03→0.069 can vary significantly, as discussed below. The increased global-mean insolation, however, is seen in all cases (not shown).

4.1.3 Argument of perihelion

For changes in the argument of perihelion, we look at the two orbital geometries {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0.069, $\mathit{\omega }=\mathrm{180}{}^{\circ }$} and {$\mathit{\epsilon }=\mathrm{22}{}^{\circ }$, e=0.069, $\mathit{\omega }=\mathrm{90}{}^{\circ }$} (Fig. 8c(i)–(v)), where the choices of ω represent the largest difference in surface-air temperature for a given ε and e (see Fig. A6).

Shifting ω changes neither the meridional distribution nor the annual average global amount of incoming solar flux (Fig. 8c(ii)). However, it affects the seasonal differences in insolation. At $\mathit{\omega }=\mathrm{180}{}^{\circ }$, perihelion and aphelion occur in March and September, respectively, i.e. at times with intermediate insolation in both hemispheres. In comparison, at $\mathit{\omega }=\mathrm{90}{}^{\circ }$ – perihelion in December, aphelion in June – the seasonality of the solar flux is amplified in the Southern and dampened in the Northern Hemisphere. Thus, during the ice-melting seasons (spring through summer), the configuration $\mathit{\omega }=\mathrm{90}{}^{\circ }$ lowers the incoming solar flux in the Northern and raises it in the Southern Hemisphere, which affects the ice-melting rates and contributes to the observed hemispheric warming–cooling asymmetry. Since the absolute increase in solar flux is larger in the Southern Hemisphere than the absolute decrease in solar flux in the Northern Hemisphere, the SH warming turns out stronger than the NH cooling. This is the main reason for the overall warming upon shifting $\mathit{\omega }=\mathrm{180}{}^{\circ }\to \mathrm{90}{}^{\circ }$.

Annual average short-wave reflection from snow in southern middle and high latitudes, particularly over the largest continent, is reduced at $\mathit{\omega }=\mathrm{90}{}^{\circ }$ (Fig. 8c(iii)). This contributes to a significant warming at these latitudes and is a plausible reason for the configurations with $\mathit{\omega }=\mathrm{90}{}^{\circ }$ being slightly warmer than those with $\mathit{\omega }=\mathrm{270}{}^{\circ }$. In order to test this hypothesis, we ran a number of simulations at varying ω with the continents mirrored at the Equator. A comparison with the original results (Fig. A9) shows that the difference between the temperature maxima at ω=90 and 270^∘ can be explained by the continental configuration and, in particular, the location of the largest continent.

Changes in the meridional heat flux counteract the short-wave contributions around the ice margins and around the North Pole but play a significant role in warming the tropics and the South Pole.

Emissivity differences yield temperature changes of a smaller magnitude and the same sign as the total temperature difference, again reflecting the positive water-vapour and stratus-cloud feedbacks.

Comparing the configurations with ω=180 and 270^∘, the results look very similar except for the strong warming that occurs in the Northern Hemisphere.