A coupled climate model with idealized representations of atmosphere, ocean, sea ice, and land is used to investigate transitions between global climate equilibria. The model supports the presence of climates with limited ice cover (Warm), a continuum of climates in which sea ice extends down into the midlatitudes and the tropics (Cold), together with a completely ice-covered earth (Snowball). Transitions between these states are triggered through volcanic eruptions, where the radiative effect of stratospheric sulfur emissions is idealized as an impulse reduction in incoming solar radiation. Snowball transitions starting from the Cold state are more favorable than from the Warm state, because less energy must be extracted from the system. However, even when starting from a Cold climate, Toba-like volcanic events (cooling of order −100 W m−2) must be sustained continuously for several decades to glaciate the entire planet. When the deep ocean is involved, the volcanic response is characterized by relaxation time scales spanning hundreds to thousands of years. If the interval between successive eruptions is significantly shorter (years to decades) than the ocean’s characteristic time scales, the cumulative cooling can build over time and initiate a state transition. The model exhibits a single hysteresis loop that connects all three climate equilibria. When starting from a Snowball, the model cannot access the Cold branch without first transitioning to an ice-free climate and completing the hysteresis loop. By contrast, a Cold state, when warmed, transitions to the Warm equilibrium without any hysteresis.
Paleoclimate reconstructions have shown that Earth has experienced a wide array of climates over the last billion years. At its most extreme, the climate has ranged between hothouse worlds with little or no ice at the poles and deep ocean temperatures upward of 20°C (e.g., Tripati and Elderfield 2005), to Snowball periods in which ice covered most of Earth’s surface (e.g., Hoffman et al. 1998). The ability of the climate to exist in a number of dynamical states was discussed by Stommel (1961) with regards to the bistability of the oceanic meridional overturning circulation (MOC). This work suggested that the MOC could exhibit two stable states, namely, a vigorous circulation mode driven by temperature and a weak mode controlled by salinity. Another source of instability in the climate is the ice–albedo feedback, which was first investigated by Budyko (1969) and Sellers (1969) in simple energy balance models (EBMs). These studies found that the Earth could exist in a stable state with either a small amount of ice at the poles or with complete ice cover. Rose and Marshall (2009) modified the EBM to demonstrate that an additional state with a large ice cap can be stabilized by the meridional convergence of heat transported by the ocean subtropical cells. Indeed several general circulation models (GCMs) have shown that climates with ice extending down to the midlatitudes or the tropics are realizable when ocean heat transport (OHT) can arrest the advance of sea ice at those latitudes (e.g., Poulsen and Jacob 2004; Langen and Alexeev 2004; Ferreira et al. 2011; Rose et al. 2013) or when the meridional gradient in ice albedo can limit the strength of the positive ice–albedo feedback (Abbot et al. 2011).
Transitions between these multiple equilibria require global-scale forcings that push the climate away from its original equilibrium to the next stable state. Many of the long-lasting climatic shifts that occurred in the last billion years were governed by the balance between weathering of silicate rocks and the input of carbon dioxide by volcanic eruptions in the atmosphere, which typically spans 105–106 yr time scales (Walker et al. 1981). However, large perturbations in the atmosphere’s radiation budget acting over much shorter (~decadal) time scales could have also caused severe episodes of climate change. Perhaps one of the largest forcings ever experienced by Earth’s climate was the Chicxulub asteroid impact 66 million years ago, during which the amount of dust and sulfate aerosols in the atmosphere could have blocked more than half of the solar radiation reaching Earth’s surface for several years, subsequently leading to extreme cooling and the beginning of a mass extinction (Brugger et al. 2017; Kaiho and Oshima 2017).
The volcanic emission of sulfur particles into the stratosphere is also known to cool the surface of the Earth for several years after the eruption (Robock 2000). Over the last few thousand years, these volcanic forcings were too weak and short-lived to have caused a drastic shift in climatic state. However, there have been recorded instances of volcanic activity much larger than that of the Holocene, which may have triggered a significant glaciation (e.g., Rampino and Self 1992; Prueher and Rea 2001). The Toba supereruption around 73.5 ka before present had a sulfur loading that was two orders of magnitude larger than the 1991 Mount Pinatubo eruption, and could have reduced incoming solar radiation by a third for several years after the event (Timmreck et al. 2010). While GCM simulations with these types of forcings did not yield widespread glaciation (Jones et al. 2005), some authors have argued that a sequence of supereruptions separated by a short time interval could eject enough sulfur particles into the stratosphere to trigger a significant shift in Earth’s history, such as a mass extinction (Baresel et al. 2017) or a Snowball Earth (Stern et al. 2008; Macdonald and Wordsworth 2017). In particular, Macdonald and Wordsworth (2017) argue that starting from the cool Neoproterozoic background climate, a succession of supereruptions occurring over several decades could have initiated the Snowball event of the Sturtian (717–616 Ma). While this estimate was obtained using a single-column climate model with a heat capacity corresponding to the mixed layer depth of the ocean (~50 m), Voigt and Marotzke (2010) and Voigt et al. (2011) argue that the entire ocean must cool to the freezing point to initiate a Snowball climate. In the present study, we explore the extent to which the ocean constrains the forcing magnitude and time scale required for global climate transitions using a coupled atmosphere–ocean–ice GCM (which we refer to here as the “climate model”), and explore whether volcanic eruptions could successfully initiate a drastic climatic state change, such as a Snowball transition.
The presence of multiple equilibria in the global climate can also affect the transition behavior of the system. While Voigt and Marotzke (2010) did not find intermediate stable states between their reference “warm” climate and the Snowball, the following work of Voigt et al. (2011), Voigt and Abbot (2012), Yang et al. (2012), and others did find such states. The climate model used in this study allows the existence of a continuum of climates with sea ice extending from the midlatitudes to the tropics (45° to 25° latitude), owing to the stabilizing nature of meridional heat transport by subtropical oceanic cells (e.g., Rose and Marshall 2009; Ferreira et al. 2011). We initiate transitions starting from these different states and evaluate how the proximity to a transition threshold affects the time scale required to tip the climate to a new state. Moreover, a number of studies have suggested that a Waterbelt (or Slushball) climate with ice extending down to the tropics could better explain the survival of oceanic life during the Sturtian and Marinoan (650–635 Ma) glaciations events, as opposed to a complete Snowball (e.g., Hoffman et al. 2017). The climate model used in this study shows that a Waterbelt-like climate may occur for a significantly smaller threshold forcing than that required for complete glaciation, but unlike the Snowball, this state exhibits no hysteresis when the forcing is removed.
2. The coupled climate model and modeling framework
Numerical experiments are carried out with the coupled atmosphere, ocean, sea ice, and land model based on the MITgcm (Marshall et al. 1997a,b). The particular setup of the model used here comprises two 45° wide flat continents separated by a 90° angle such that the ocean is split into a narrow Atlantic-like basin and a wide Pacific-like basin connected to an unblocked “Southern Ocean” in the South (see Ferreira et al. 2018). The ocean is flat bottomed with a uniform depth of 3000 m. The atmospheric scheme is based on the simplified parameterizations of primitive equation dynamics (SPEEDY) from Molteni (2003) at five-level vertical resolution. It includes a four-band radiation scheme, a parameterization of moist convection, diagnostic clouds, and a boundary layer scheme. In the reference states, the solar constant (divided by 4) is 341.5 W m−2. Despite the idealized nature of the model, it captures many of the essential large-scale features of the climate such as an overturning circulation predominantly in the small basin, atmospheric storm tracks and hydrological cycle, gyres and a circumpolar current in the oceans, and a seasonal sea ice cycle. Additional modeling details are provided in the online supplemental material.
Perturbation experiments are conducted by reducing the top of the atmosphere (TOA) insolation. The forcing F is diagnosed as the net radiative imbalance at the TOA, averaged over the first year of the simulation. For step forcing experiments, the insolation is changed abruptly at t = 0 and subsequently held constant in time. For volcanic forcing experiments, we crudely emulate the effects of stratospheric aerosol injection by uniformly reducing the insolation for 1 yr and subsequently turning the forcing off (as in Gupta and Marshall 2018). We also conduct simulations comprising successive eruptions separated by a specified time interval.
3. Stable climate states of a coupled climate model
a. Equilibrium solutions
Figure 1 shows the three equilibrium branches found in our model: Warm, Cold, and Snowball. Note that we reserve the term “state” for individual points on a given equilibrium branch. The Warm reference state has some sea ice in the Southern Hemisphere (SH) but none in the Northern Hemisphere (NH). It has a globally averaged ocean potential temperature of θ = 8.3°C and a global mean sea surface temperature (SST) of 19.6°C. The Cold equilibrium consists of a continuum of states with sea ice edge latitude ranging from 45° to 25° and θ between 0.5° and −1.7°C. The mildest states on this Cold branch have an ice cap extending down to the midlatitudes, which is reminiscent of the glacial climate of the Pleistocene (Ferreira et al. 2018). The most extreme states on the Cold branch have sea ice reaching down to the tropics, which is analogous to Waterbelt or Slushball states found by other studies (e.g., Rose 2015; Abbot et al. 2011). Since these Waterbelt-like climates are plausible analogs for the state of the Earth during the Neoproterozoic glaciations, it is instructive to treat them somewhat separately than the rest of the Cold branch, particularly in the context of state transitions. In what follows, we thus choose to refer to Cold states with θ < −1.5°C as Waterbelts, and others simply as Cold. Our discussion does not depend strongly on this definition, as long as it is chosen close enough to the freezing point. Finally, the Snowball is fully ice covered with θ = −1.8°C (the mean freezing point of the ocean). Simulations that reach this state are stopped when the whole ocean reaches the freezing temperature, to avoid runaway ice growth in the absence of geothermal heat flux.
Figure 1 also presents the various equilibrium states in terms of their energy levels E, approximated by the sum of ocean heat content and latent energy stored in the sea ice, all measured relative to the Warm state.
where ρ is the seawater density, c is the thermal heat capacity of water, Htot is the globally area-averaged ocean depth, θ is the ocean mean potential temperature, A is Earth’s surface area, Lf is the latent heat capacity of ice, M is the total mass of ice, and the subscript w indicates the Warm reference state. The energy stored in the ice only represents at most 10% of E. Therefore, θ and E can be used almost equivalently to describe the energy content of the system. E is an energy per unit area, which we express in (W m−2 yr) (=3.1 × 107 J m−2) to aid the comparison with the energy extracted by volcanic eruptions. For reference, the 1991 Mount Pinatubo eruption imposed a globally averaged forcing of approximately −4 W m−2 at the TOA for approximately 1 yr, thus equivalent to an energy extraction of −4 W m−2 yr.
Figure 1 shows that the Warm and Cold reference states have a large energetic gap (~2500 W m−2 yr), which is a significant barrier to overcome for any transition between them. The energy gap between the Cold reference state and Waterbelt climates is smaller (~700 W m−2 yr) but significant enough to distinguish them energetically, owing to the existence of a tropical thermocline in the former, but not in the latter. On the other hand, the Waterbelt and Snowball states are not energetically distinguishable, because most of the ocean is at freezing temperature in both states.
Figure 2 shows cross sections of ocean potential temperature and overturning streamfunction, along with the sea ice edge latitude and OHT for various stable states of our model. The shallow subtropical cells (STCs) are driven by the trade winds, which act upon the temperature gradient of the thermocline to transport heat to the midlatitudes. This heat flux allows the stabilization of the sea ice edge and facilitates the realization of multiple equilibria in the model (see Rose and Marshall 2009; Ferreira et al. 2011; Rose 2015). In the Warm state, the OHT reaches high latitudes in the NH, but only midlatitudes in the SH. The temperature of the deep ocean (below 1000 m) is everywhere above 6°C, including at the poles where fresh surface waters allow a temperature inversion with a surface colder than the abyss. In the Cold state, all the waters below 1200 m and poleward of 50° latitude in both hemispheres (~50% of the total ocean volume) are at the freezing point. This results from temperature-induced convection having to bring the entire water column to the freezing point before ice can grow at the surface (see Rose et al. 2013). The thermocline, however, remains well stratified, supporting poleward heat transport by the subtropical cells, which arrests the sea ice in the midlatitudes. The maximum strength of poleward heat transport in both the Warm and Cold states is of order 2 PW, peaking around 20° latitude. In the Waterbelt state, the thermocline stratification is eroded and the sea ice extends farther equatorward, until only a shallow (~300 m) surface layer of the tropical ocean remains above freezing and the ice edge reaches around 25° latitude in both hemispheres. The strength of poleward heat transport falls to 1.2 PW and the peak moves equatorward to 10° latitude. Finally, in the Snowball, OHT is essentially reduced to zero and the entire ocean is at the freezing point. In the following section, we explore the time scales associated with changes between these states.
b. Triggering transitions between equilibrium states
In this section, we conduct step (constant) forcing simulations in our climate model with the aim of initiating transitions between the various climate states presented in section 3a. Figure 3 shows these results on a bifurcation diagram, where each scatter represents the NH sea ice latitude in the equilibrated state of a given simulation with forcing F. This forcing magnitude is expressed with respect to the reference solar constant S0 = 341.5 W m−2. Simulations start from either the Warm (red circles), Cold (blue squares), Waterbelt (green triangle), or Snowball (black diamond) states. The solid lines describe the range over which these states exist, with the same color convention. The individual states labeled W (Warm), C (Cold), and S (Snowball) are chosen as reference climates because they exist for F = 0.
1) Warm start
Starting from W, a small step forcing of −1.5 W m−2 induces a slow transition to the Cold equilibrium over several thousand years. This implies that the Warm reference state is at the edge of a critical threshold and is thus particularly sensitive to changes in its radiative budget. Larger forcings of −7 and −11 W m−2 both lead to a Cold equilibrium and a forcing of −36 W m−2 produces a Snowball. We also force the Warm reference state with a positive forcing of +3 W m−2, which leads to a Very Warm (VW) climate that has no ice anywhere on the globe, a global mean ocean potential temperature of 14°C and a global mean SST of 24°C. We find no stable states between 90° and 45° NH sea ice extent. The absence of a stable climate with a small ice cap, known as the small ice instability (SICI), has been discussed in the context of simple EBMs (e.g., Held and Suarez 1975; North 1984) and more sophisticated models (e.g., Huang and Bowman 1992; Matteucci 1993; Lee and North 1995; Morales Maqueda et al. 1998; Langen and Alexeev 2004). While these studies concluded that the SICI tends to disappear with better representation of continental geometry, unforced variability, ice sheets, and ice dynamics, further work would be required to investigate the lack of a small ice cap climate in the model used here.
2) Cold start
Starting from C, we find that step forcings of −2, −3, −5, and −9 W m−2 lead to progressively colder climates without abrupt state transitions. For a forcing Fc/s = −20 W m−2, the climate initially transitions to a Cold state that remains stable for about 150 yr, but then collapses into a Snowball. Between the forcings Fwb = −5 W m−2 and Fc/s, we find Cold climates with θ < −1.5°C and tropical sea ice, which we refer to as Waterbelts. Starting from such a Waterbelt climate (WB) and relaxing the forcing to 0, we find that the system returns to C. Finally, we conduct warming experiments starting from C with F = +1.5 and +8 W m−2, which lead back to a Warm climate. We conclude that there is no identifiable hysteresis loop between the Cold and Warm branches.
3) Snowball start
Starting from a Snowball, we relax F to 0 and find that the climate remains in a Snowball state (S). Thus, the Cold branch is not directly accessible from a Snowball state. Starting from S, we attempt to exit the Snowball state by imposing a forcing of +7 W m−2, but find that the surface remains frozen everywhere. We expect that a larger positive forcing could overcome the strong ice–albedo feedback and bring the system back to the Warm branch, but the magnitude of this forcing is not quantified in this study. Note that if this threshold is high enough, the climate may not stabilize into a Warm state, but instead transition directly to a moist or runaway greenhouse (Yang et al. 2017).
The plot on the right of Fig. 3 summarizes these results in a simplified schematic (not to scale). The arrows show two hypothetical paths along the bifurcation diagram. Starting from Warm (red arrows) and slowly reducing the solar constant, the system would reach the Warm–Cold transition threshold Fw/c ~ −1.5 W m−2 and move down to the Cold branch. The system would then smoothly transition to a Waterbelt-like climate around Fwb ~ −5 W m−2, before reaching the Cold–Snowball threshold Fc/s ~ −20 W m−2. From a Snowball state, a large increase in the solar constant would be required to switch back to the Warm state and complete the loop. When starting from the Cold branch (blue arrows) and increasing the solar constant, the system can move to the Warm branch without any hysteresis. In the appendix, we further interpret this bifurcation diagram in terms of a 1D EBM developed by Rose (2015), which shows characteristics that are reminiscent of our climate model.
4. Transition time scales between equilibrium states
a. Transitions in response to step forcings
Figure 4 shows state transition times in the climate model for the step simulations starting from the Warm (red) and the Cold (blue) reference states as a function of the forcing magnitude. These are the same experiments as used in Fig. 3. State transitions to both a Waterbelt and a Snowball are considered. As expected, the transition time shortens as the forcing magnitude increases, and there is a minimum threshold forcing below which no transition occurs. Transitions from the Warm state take longer than from Cold (for the same forcing) because more energy must be extracted in the process. Voigt and Marotzke (2010) (hereafter, VM2010) argued that as the planet approaches a Snowball, deep convection and the transport of cold waters into the abyss eventually lead to a well-mixed ocean. This prompted them to adopt the following 1-box model (here expressed in terms of anomalies) to estimate state transition times:
where F is the TOA forcing, T is the temperature anomaly relative to the starting state, λ is the climate feedback parameter, and H is the ocean depth relevant to the transition. The model does not take account of the change in albedo as the transition occurs since all parameters are constant in time. The solution for the step response (constant F) of Eq. (2) is
where the time scale τ is given by
and the equilibrium climate sensitivity is
where here F is negative, corresponding to extraction of energy.
As we now explore, state transition times in our climate model are well captured by the 1-box model, when considering the difference in potential temperature between the initial and final states Δθ. The time taken by the 1-box model to achieve a temperature difference Δθ through a step forcing F is obtained from Eq. (3) as follows:
From Eq. (5), the forcing required to obtain a temperature response Δθ at equilibrium is
It follows that forcings weaker than Fmin cannot achieve this temperature change, even when sustained indefinitely. Therefore, F must be larger than or equal to Fmin to achieve a temperature change Δθ.
One can also estimate the efficiency of energy transfer required to achieve Δθ, as a function of forcing magnitude. The total energy Etot supplied by the forcing during the transition is
From Eq. (1), the net energy ΔEnet extracted from the system for a temperature change Δθ is
where we ignore the energy stored in the atmosphere, land and sea ice. The efficiency of energy extraction ε is then
As the forcing amplitude increases (F → −∞), the efficiency ε tends to 1, but drops to 0 as F tends toward Fmin. For large forcings, climatic feedbacks have limited time to act, and hence the energy extraction occurs efficiently. Conversely, as the forcing magnitude decreases, the climatic feedbacks have enough time to counter the effect of the forcing, and as F drops below Fmin, the feedbacks become too strong to allow a temperature change Δθ altogether. In the following paragraphs, we investigate to what extent these simple ideas apply to the climate model transitions.
For simulations starting from Warm, the temperature anomalies spread into every part of the ocean, including the abyss. Thus, the relevant depth H of the 1-box model is Htot = 2416 m (the total area-averaged depth of the ocean). The climate sensitivity parameter in the Warm reference state is found using the method outlined in Gregory et al. (2004) (details in supplemental material), which gives λw = 0.72 W m−2 K−1. This value is associated with the SST response of the model, but the assumption of a well-mixed ocean allows us to extend it to the ocean-mean temperature. For a transition to a Snowball, the initial temperature is θi = 8.3°C and the final temperature is θf = −1.8°C, the freezing point of water. Inserting these values into Eq. (6) gives the red curve in Fig. 4, which is a reasonable approximation to the coupled model results.
Figure 4 also shows (in yellow squares) the Snowball transition times obtained by VM2010, when starting from a modern-day climate in a comprehensive GCM. To facilitate comparison with our simulations starting from Warm, we scale the VM2010 time scales by a factor fvm, to account for the different heat capacities of the two models:
where Hvm = 2603 m is a realistic globally averaged depth for the modern ocean (including land surfaces) and θvm = 4.4°C, as reported in VM2010. This gives a scaling parameter fvm = 1.5. Figure 4 shows that the transition time scales obtained by VM2010 are larger than the ones corresponding to the Warm start simulations and this difference is accentuated by the factor fvm. Using the 1-box model with the parameters provided by VM2010 (yellow line), one finds the effective climate feedback parameter in their model to be λvm = 3.3 W m−2 K−1 (see supplemental material). This is significantly larger than λw and results in longer transition times. Note that λvm is larger than the climate feedback parameters typically reported by state-of-the-art GCMs (~1 W m−2 K−1) such as the one used in VM2010. However, a direct comparison may not be possible because λvm is obtained from a 1-box model fit to the Snowball transition times associated with θ, rather than from the more conventional “Gregory” method associated with surface temperatures. Therefore, λvm may encapsulate different climate processes than typical feedback parameters. In particular, the representation of nonlinear atmospheric feedbacks (e.g., clouds, water vapor) and the restructuring of the ocean may differ.
In the simulations starting from Cold, the temperature anomaly spreads into the part of the ocean that is not already at the freezing temperature, namely, the tropical thermocline. The relevant ocean depth is the top 1000 m (approximately), which has an average potential temperature of 3.1°C. The climate feedback parameter corresponding to the Cold reference state is again obtained by the Gregory method and found to be somewhat higher than for Warm: λc = 0.95 W m−2 K−1 (see details in supplemental material). This corresponds to a time scale of 159 yr (see section 1). Using these values along with θi = 3.1°C and θf = −1.8°C in Eq. (6) gives the solid green line in Fig. 4. Overall, the box model can capture the climate model’s transition times starting from Cold, but tends to underestimate them as F approaches the threshold forcing.
b. Transitions in response to volcanic impulse forcings
1) Single eruptions
We now investigate the behavior of the climate model following large idealized volcanic eruptions. The largest eruption simulated from the Warm state has a magnitude of F = −105 W m−2, which is on the same order as estimates of the Toba mega-eruption (e.g., Jones et al. 2005). In the year following the simulated eruption, the global-mean SST drops by 5.4°C, the global sea ice fraction increases by 2% and both quantities recover their original values within 75 yr. We also conduct a simulation of comparable magnitude (F = −89 W m−2) starting from Cold, which leads to a global-mean SST drop of 3.1°C, a shift in the ice edge latitude from 44.4° to 42.2° (a 4% increase in sea ice fraction) and a recovery within 400 yr. Despite the large magnitude of the forcings considered here, the net energy extracted from the system (~100 W m−2 yr) is not large enough to cause a shift in climate equilibrium (see Fig. 1). The longer recovery time for SST and sea ice fraction in the Cold start eruption is likely due to a stronger ice albedo feedback, since the Cold state has a larger sea ice fraction (33%) than the Warm state (11%). However, we now show that unlike the recovery times of surface quantities, the time scales associated with the global mean ocean potential temperature θ are larger in the Warm state than the Cold one.
Figure 5 shows time series of θ for all the single eruption simulations conducted in this study. The experiments are run until θ relaxes back to within the 2σ levels of unforced variability. One can estimate the peak cooling θpeak at the end of the first year by equating the energy extracted by the forcing to the change in heat content of the ocean, and neglecting the effects of the feedbacks as follows:
where Δt is the forcing duration, which is assumed small compared to the time scale associated with climatic feedbacks. Equation (11) is a good approximation for all the simulations shown in Fig. 5, irrespective of how deep the anomaly penetrates into the ocean.
For Cold starts, we consider two eruptions with F = −89 and −29 W m−2 (blue and purple curves, respectively). Because the deep ocean is already at the freezing temperature at the starting state, the temperature anomaly remains within the tropical thermocline. When normalized with respect to θpeak (Fig. 5b), both responses follow the same decay, which can be fit by a decaying exponential with an e-folding time of τc = 159 yr (blue dotted line). According to the 1-box model, the ocean depth associated with this time scale is [see Eq. (4)]
with λc = 0.95 W m−2 K−1 corresponding to the climate feedback for the Cold state. This gives Hc = 1149 m (roughly 47% of the total area-averaged ocean depth).
For Warm start simulations, Fig. 5b shows that forcing pulses of magnitude F = −35 and −105 W m−2 lead to a relaxation back to equilibrium over several thousand years (orange and red curves, respectively). The normalized plot shows that the two curves follow a similar path to recovery, but the one corresponding to F = −35 W m−2 reaches background noise levels within 2000 yr, versus at least 5000 yr for F = −105 W m−2. We therefore choose to fit to the F = −105 W m−2 experiment and find that the following two-exponential model gives a good approximation of the impulse response (red dotted line):
with T1 = −0.10°C, τw1 = 500 yr, T2 = −0.13°C, and τw2 = 5911 yr.
The millennial-scale decay τw2 cannot be captured by the 1-box model, since the longest time scale that can arise from Eq. (4) is 443 yr, when using the full ocean depth H = 2416 m and the inferred climate sensitivity parameter λw = 0.72 W m−2 K−1. To better understand the origin of τw2, we conduct an additional three numerical experiments with (i) a smaller magnitude F, (ii) a positive F, and (iii) a large negative F starting from the Very Warm state described in section 3:
When F is reduced to −11 W m−2 with a Warm start (yellow curve in Fig. 5a), the magnitude of the perturbation is not large enough for the θ anomaly to remain outside the range of natural variability for more than a couple of decades. The ocean temperature anomalies remain in the top 1000 m and are rapidly dissipated by climatic processes near the surface. This suggests that there is a minimum amount of forcing required (somewhere between −11 and −35 W m−2) to induce the millennial time scale of decay in θ.
A positive forcing experiment was conducted with F = +31 W m−2 (black curve), which likely does not correspond to a particular geophysical mechanism, but is useful for assessing the asymmetry in our model’s response between warming and cooling. Unlike in the corresponding negative forcing case (−35 W m−2), the decay of θ does not have a millennial time scale as it reaches background noise levels within less than 1000 yr.
To assess the sensitivity to the starting state, a mega-eruption simulation with F = −108 W m−2 was conducted from the Very Warm state (pink curve). In this experiment, θ responded with an e-folding time of 437 yr, which corresponds to the decay time scale of the 1-box model when considering the full depth of the ocean (pink dotted line).
In the following paragraph, we briefly explore aspects of ocean dynamics that can lead to the millennial decay time scale τw2. Figure 6 displays vertical profiles of ocean potential temperature for a subset of the single eruptions simulations conducted for this study. Figure 6a shows the F = −105 W m−2 simulation starting from Warm. After 100 yr, the temperature profile exhibits two peaks in cooling: the first around 500-m depth and the second at the bottom of the ocean. The shallower peak is likely a result of anomalously cold water parcels being driven down into the thermocline by Ekman pumping in the midlatitudes. This mechanism can shield the cold temperature anomaly away from atmospheric damping processes acting at the surface for decades to centuries (e.g., Stouffer 2004; Stenchikov et al. 2009; Gupta and Marshall 2018). The peak cooling in the abyss is set by enhanced deep convection around the poles, particularly at the sea ice margin in the SH (see Fig. S2). After 500 yr, the peak cooling in the thermocline is damped away by a combination of atmospheric feedbacks and mixing into the deeper ocean. At that time, the temperature anomaly profile is almost linear with depth and the mixed layer temperature anomaly reaches unforced variability levels (−0.1°C). The subsequent evolution of the temperature profile occurs over 1000-yr time scales, as the deep ocean temperature anomaly slowly diffuses back up toward the surface. This upward diffusion process is illustrated using a simple 1D ocean diffusive model with constant diffusivity κ = 3 × 10−5 m2 s−1 (the uniform background value used in our climate model) and boundary conditions ensuring no flux at the bottom of the ocean and a temperature anomaly fixed to 0 at the top. The latter implicitly assumes that climatic feedbacks at the surface act on much shorter time scales than the various ocean processes bringing the temperature anomaly back up. We initialize the 1D model with a linear temperature profile, roughly corresponding to the temperature anomaly profile found in our climate model at t = 500 yr. The 1D model temperature evolution (in the dotted lines) shows a roughly similar relaxation behavior to the climate model and has a time scale that scales with H2/κ, which is O(1000) yr.
Figure 6b shows θ anomaly profiles for the F = −3.5 W m−2 impulse simulation starting from Warm. At year 100, there is again a peak cooling around 500-m depth and a small peak in the abyss. Figure S2 shows signatures of Ekman pumping transporting cold waters into the thermocline and enhanced deep convection in the SH high latitudes, although not as large as in Fig. 6a. Subsequently, diffusive processes flatten the temperature profile and the deep temperature anomalies reach back up to the surface over several thousands of years.
Figure 6c shows temperature anomaly profiles for F = +31 W m−2 starting from Warm. After 100 yr, midlatitude Ekman pumping has driven warm temperature anomalies down into the thermocline. However, a smaller fraction of the anomaly than in Fig. 6b penetrates into the deep ocean because of the higher buoyancy of these anomalously warm waters. This results in an overall ocean relaxation time scale of several hundred years only.
Figure 6d shows temperature anomaly profiles for F = −108 W m−2 starting from Very Warm. After 100 yr, there is a large peak cooling at the surface and a small one in the abyss. Some of the temperature anomaly has penetrated into the thermocline, but the surface cooling has not been damped as significantly as in Fig. 6a. To explain this behavior, one should note that the Very Warm state does not contain ice anywhere on the globe. The impulse cooling promotes sea ice formation in the SH, but the resulting deep convection at the ice margin is not strong enough to carry a large volume of cold waters into the abyss (see Fig. S2). In the midlatitudes and the tropics, a fraction of the temperature anomaly does penetrate into the deeper ocean, but it is a smaller share of the total cooling than in Figs. 6a and 6b.
In summary, the millennial decay time scale τw2 is a result of cold temperature anomalies reaching the deepest layers of the ocean and slowly diffusing upward over several thousand years. These anomalies can reach the seafloor by enhanced deep convection in the SH high latitudes over the first few decades after the eruption. Notably, the millennial time scale does not arise for a positive impulse or when starting from the Very Warm state.
2) Succession of eruptions
We now investigate how the long oceanic time scales can favor a large buildup of the cooling response and potentially initiate state transitions when forced by repeated volcanic eruptions. We first establish whether the transition behavior for volcanic eruptions of magnitude F and duration Δt (assumed 1 yr) repeated every τi years is equivalent to that of a step simulation with average forcing Fav:
Figure 7 shows the globally averaged θ and sea ice fraction evolution for two simulations in which the model is forced by an idealized eruption every τi = 10 yr, with F = −105 W m−2 starting from Warm and F = −89 W m−2 starting from Cold. For comparison, we also consider the corresponding step simulations with Fav = −11 W m−2 starting from Warm and Fav = −9 W m−2 starting from Cold [using Eq. (14)]. The results show that the θ and sea ice fraction time series obtained from the repeated eruption simulations closely follow the corresponding step simulations, for both Warm and Cold starts. All simulations eventually transition to a Waterbelt state with approximately 60% ice fraction and the sea ice edge stabilizes around 27° latitude in both hemispheres. The step response is a good approximation for repeated eruptions because the interval between them is much smaller than the ocean relaxation time scales, and because each incremental cooling from a given eruption is much smaller than the overall Δθ.
Figure 7 also shows that the 1-box model of Eq. (2) faithfully represents the evolution of θ in all the cases considered. There are nonlinearities associated with the transition, but since those only significantly affect θ near the freezing point (θf = −1.8°C), the 1-box model provides a good estimate of the state transition times. For the Cold start simulations, both single eruption and step simulations were well approximated by the 1-box model with time scale τc = 159 yr and a climate feedback parameter λc = 0.95 W m−2 K−1. For Warm start simulations, however, there is an apparent discrepancy between the two-time-scale model used in Fig. 5 and the single time-scale model used in Fig. 7. This is resolved by noting that here, a state transition occurs after 550 yr, which is too short to activate the millennial time scale. Therefore, the 1-box model is a good approximation for transitions occurring within several centuries, but otherwise the two-time-scale model is more appropriate.
In Fig. 8, we explore the state transition times for repeated uniform eruptions over a range of impulse forcings F (lasting 1 yr) and eruption intervals τi. One can estimate the transition threshold by assuming a step response with constant forcing Fav because the system time scales are much longer than the interval between eruptions. This implies that scenarios with the same Fav have the same transition time, which results in the straight contours in Fig. 8. The dotted lines show the following transition threshold boundaries: Warm–Cold in black (Fw/c = −1.5 W m−2), Cold–Waterbelt in red (Fwb = −5 W m−2), and Cold–Snowball in white (Fc/s = −20 W m−2).
For the Cold start case, we assume a 1-box model with time scale τc = 159 yr and can thus estimate transition times to either a Waterbelt or Snowball climate using Eq. (6) and the parameters calculated in section 4a, with F = Fav. As noted in that section, the 1-box model tends to underestimate the GCM transition time scales for Cold start simulations when F is of low magnitude, so the transition times of Fig. 8a should be regarded as a lower bound. When the repeated eruptions start from Warm, the transition times can be estimated similarly, but in this case there are two time scales (τw1 and τw2) involved in the response. When the forcing is large enough for a transition to occur within several hundred years, Fig. 7 showed that the 1-box model was a suitable approximation. However, a more general estimate of the step response can be obtained by integrating Eq. (13) and scaling it appropriately to obtain the two-time-scale step response from Warm:
where Fref = −105 W m−2. The time taken to transition from an initial temperature θi = 8.3°C to the final temperature θf is calculated numerically using Eq. (15) and shown in Fig. 8b. For transitions to either a Waterbelt or Snowball climate, the final temperature is θf = −1.8°C (the freezing point), whereas for transitions to the Cold state it is θf = 0.5°C. Note that we only consider transition times smaller than 104 yr, because over longer time scales, greenhouse gas emissions from volcanic eruptions may start to overcome the cooling (e.g., Walker et al. 1981).
Figure 8 shows that transition to a Snowball requires a sequence of almost continuous Toba-like events (~−100 W m−2) for at least several decades, even when starting from a Cold climate. This differs from the results of Macdonald and Wordsworth (2017), who suggested that such a transition could occur within 3 yr for volcanic eruptions with a peak forcing as weak as −10 W m−2. This difference arises because the heat capacity of their model consists of the mixed layer depth of the ocean (50 m), whereas in our climate model, the entire thermocline (~1000 m) plays a role when starting from Cold. Transitions to a Waterbelt climate may occur for a weaker Fav, but over several hundred years for a Cold start and around 1000 or 2000 years for a Warm start. In this climate model, such Waterbelt-like states do not exist for F = 0 (see Fig. 3), so the system would relax back to Cold when the eruptions stopped. Nevertheless, the hysteresis associated with Waterbelt states may differ between models (e.g., Rose 2015), and therefore the possibility of a long-term transition to such a state is still plausible. Finally, volcano-induced transitions from Warm to Cold may occur for small-magnitude forcings (~−1.5 W m−2) because of the proximity of the Warm reference state to a transition threshold. However, such transitions take several thousands of years, during which the buildup of greenhouse gases from intense volcanic activity could start mitigating the cooling effect.
Our study has explored transitions between the global equilibria of a complex coupled climate model with representations of atmosphere, ocean, sea ice, and land. We study whether cooling associated with volcanic eruptions—acting alone (single pulse) or in concert (repeating pulse or step)—could trigger transitions between equilibrium states. Although the models employed are rather idealized, we offer the following general conclusions that hopefully do not depend on those idealizations:
A series of Toba-like events (F ~ −100 W m−2) must occur almost continuously for a minimum of several decades to initiate a transition to a Snowball climate by volcanic cooling acting alone. This is a significantly longer forcing period than reported by Macdonald and Wordsworth (2017), who ignore the thermal capacity of the ocean beneath the mixed layer. In our model, the entire water column must reach the freezing temperature of seawater before ice can form at the surface. This result is consistent with the work of VM2010 and Voigt et al. (2011). The relevant system time scales are associated with the full ocean depth when starting from a Warm climate and the tropical thermocline when starting from a Cold one. Assuming that volcanic forcings are no greater in magnitude than Toba-like events, the interval between eruptions must be on the order of years to decades (much smaller than the oceanic time scales) to successfully initiate any state transition between them. Consequently, the transition behavior can be predicted from the time-average forcing Fav without knowledge of the detailed history of volcanic forcing.
The presence of multiple equilibria in the climate system could in theory enable a multistage transition to a Snowball. Indeed, if the climate were already in a stable Cold state with sea ice extending down into the midlatitudes, then volcanic cooling could more readily initiate a Snowball transition, since only the tropical thermocline would remain to be brought to the freezing temperature. Nevertheless, decades of intense volcanic cooling would still be required, even in this most favorable case.
The proximity of the equilibrium state to a transition threshold controls the amount of forcing required to initiate a transition. In our model, the Snowball transition threshold is Fc/s = −20 W m−2 (a 5.8% reduction in incoming solar radiation with respect to the modern constant). This is comparable to the values found in other studies (e.g., VM2010; Voigt et al. 2011; Yang et al. 2012). The Warm reference state, which resides close to a critical threshold, may evolve into a Cold state for a forcing as small as −1 W m−2, although this forcing must be sustained for several thousands of years.
The most extreme climates on the Cold branch of our bifurcation diagram are analogous to Waterbelt (or Slushball) climates, in which sea ice extends down to the tropics. The energy required to access such states is similar to that of a Snowball, since both have a mean ocean temperature close to freezing. The threshold forcing for the Waterbelt state (Fwb = −5 W m−2) is significantly lower than for a Snowball (Fc/s = −20 W m−2), but if the forcing is relaxed after a Waterbelt transition, the system smoothly returns to a climate with no tropical sea ice. This differs from the results of Rose (2015), who find two distinct branches for the Cold and Waterbelt states, and a significant amount of hysteresis associated with the latter.
A simple 1-box model can capture transition times for both Warm and Cold start simulations relatively well. However, close to the Warm–Cold transition threshold, the climate model exhibits a millennial decay time scale that cannot be explained by the assumption of a well-mixed ocean. This is caused by enhanced deep convection in the SH, which produces cold temperature anomalies in the abyssal ocean that diffuse back up to the surface over thousands of years. In this case, a 2-box model that includes this millennial decay time scale provides better estimates of transition times.
The authors acknowledge support from the MIT-GISS collaborative agreement (NASA) and are grateful for informative conversations with Brian E. J. Rose. We also thank Dorian S. Abbott, Jun Jang, and an anonymous reviewer for their insightful comments on the paper.
Interpretation in Terms of a Budyko–Sellers Energy Balance Model
A simple Budyko–Sellers EBM modified to include the effect of OHT (as in Rose and Marshall 2009) can be used to illuminate the multiple equilibria found in the climate model. We employ the following model, developed by Rose (2015), which describes the zonal-mean energy balance of the planet in terms of a characteristic surface temperature T:
where x = sin (ϕ) and ϕ is latitude, T(x) is the zonal-mean equilibrium surface temperature, Ka is a coefficient of large-scale heat diffusion for the atmosphere, B is the longwave radiative damping, Q0 is the solar constant (divided by 4), Fo(x) is OHT convergence, and a(x, T) is a nonlinear coalbedo function that depends strongly on T. The term Q0[1+s2P2(x)] gives a reasonable spatial distribution of the incoming radiation, where P2(x) is the second Legendre polynomial and s2 = −0.48. The OHT is held constant in time and parameterized as follows:
where N is a positive integer and ψ is a constant (in units of PW) setting the amplitude of OHT. The parameter N controls the shape of the OHT profile, with higher values shifting the peak of OHT and its convergence closer to the equator and thus helping the stabilization of the ice edge at lower latitudes. The ocean heat convergence is then
The coalbedo switches abruptly from ice-free value (a0) to the ice-covered (ai), depending on the local temperature as follows:
A key feature is that the threshold Txi is not tied to a fixed temperature, but depends on the local energy balance at the ice-water interface as follows:
Equation (A4) is derived by equating the vertical heat flux through the ice to the ocean heat convergence at the ice edge. The parameter hmin is the minimum sea ice thickness to have a substantial effect on the albedo and ki is the thermal conductivity of sea ice. Therefore, Txi represents the maximum surface temperature possible to maintain an ice thickness larger than hmin for a given ocean heat convergence Fo, and is used as the crude threshold between ice-free and ice-covered conditions in the EBM. Incorporating this piece of ice physics helps stabilize climates with intermediate ice covers.
Figure A1 shows a range of bifurcation diagrams obtained by solving the EBM with typical parameter values, as in Rose (2015). The analytical solution is obtained using the method detailed in North (1975) and the numerical solution is obtained with a Crank–Nicholson scheme. The numerical estimates only capture the stable (positive slope) states whereas the analytical solution gives all equilibria. Solutions are plotted for increasing values of δ = Ka/B, which measures the efficiency of atmospheric meridional heat transport relative to radiative damping. Figure A1 shows that as δ increases, the slopes become steeper and thus more prone to instability, consistent with the analytical results of Held and Suarez (1975). However, as argued by Rose and Marshall (2009), the presence of OHT can help stabilize the ice edge as δ increases. Figure A1 also shows the OHT profile for Ψ = 4 PW and N = 4, which has a peak value of 0.91 PW that occurs at 22° latitude, not untypical of our own climate. We find that δ values between 0.4 and 0.5 give a bifurcation diagram with somewhat similar characteristics to the one obtained from the climate model, including steep slopes and a continuous set of stable equilibria at intermediate latitudes. The δ = 0.5 curve also gives a bifurcation point close to Q0 = 341.5 W m−2, which is a feature seen in our climate model. However, the stable states of this curve are not actually accessible from the Warm state, as cooling from the Warm branch directly leads to a Snowball.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0883.1.s1.