Resilience of Snowball Earth to Stochastic Events

Earth went through at least two periods of global glaciation (i.e., ``Snowball Earth'' states) during the Neoproterozoic, the shortest of which (the Marinoan) may not have lasted sufficiently long for its termination to be explained by the gradual volcanic build-up of greenhouse gases in the atmosphere. Large asteroid impacts and supervolcanic eruptions have been suggested as stochastic geological events that could cause a sudden end to global glaciation via a runaway melting process. Here, we employ an energy balance climate model to simulate the evolution of Snowball Earth's surface temperature after such events. We find that even a large impactor (diameters of $d \sim 100\,\mathrm{km}$) and the supervolcanic Toba eruption ($74\,\mathrm{kyr}$ ago), are insufficient to terminate a Snowball state unless background CO$_2$ has already been driven to high levels by long-term outgassing. We suggest, according to our modeling framework, that Earth's Snowball states would have been resilient to termination by stochastic events.


Introduction
The Snowball Earth hypothesis (Kirschvink, 1992;Hoffman et al., 1998Hoffman et al., , 2017) ) proposes global surface coverage of Earth in a thick ice layer.In contrast, the most recent glacial episode, the Last Glacial Period, witnessed only partial ice coverage extending down from the poles but not reaching equatorial regions (e.g., Batchelor et al., 2019).Snowball events have occurred on at least two occasions during the Neoproterozoic era (1000 − 538.8 Ma).In particular, there is good evidence for a snowball event during the Cryogenian period (720 − 635 Ma), with other potential Snowball states also being suggested (e.g., Chumakov, 2009).Evidence for Snowball Earth's existence comes in many forms, including paleomagnetic data in Australia (Harland, 1964;Hambrey & Harland, 1981), Australian cap carbonates (Kennedy, 1996), and extreme carbon isotope excursions in Norway and Greenland (Knoll et al., 1986).Geologic measurements have been complemented by modeling efforts (e.g., Caldeira & Kasting, 1992) to understand both the origins and demise of the Snowball state.Under Earth's current insolation, a Snowball state is one of three possible stable climate regimes for Earth, alongside a temperate state and a hot steam atmosphere state (Shields et al., 2014;Turbet et al., 2021).However, as shown in Brunetti et al. (2019), the number and nature of these stable climate regimes can vary depending on feedbacks considered within the model.The Snowball state is stable against moderate variations in insolation through the high albedo of ice, which acts to prevent increases in global temperature (Shields et al., 2014;Hoffman et al., 2017).This creates a bi-stability between the temperate Earth and Snowball Earth through the difference of albedo between ice and liquid water (Ghil, 1994).It is, therefore, impossible to deglaciate Snowball Earth under such conditions with-out other climate forcings.One such forcing that is commonly invoked is the build-up of atmospheric CO 2 (Walker et al., 1981;Menou, 2015).Magmatic volatile emissions from volcanism, coupled with the absence of CO 2 sinks of Earth's temperate-state carbon cycle (e.g., silicate weathering being prevented by the covering of ice) leads to a strong greenhouse effect that warms the planet, triggering the Snowball termination (Kirschvink, 1992;Caldeira & Kasting, 1992).However, the time required to accumulate sufficient CO 2 to trigger melting (4 − 30 Myr, Hoffman et al., 1998) is too long to account for the short Marinoan glaciation (Rooney et al., 2015).Suggested alternatives include impact events (Kring, 2003) and massive volcanic eruptions (Lan et al., 2022), which we hereafter refer to as stochastic events due to the nature of their occurrences.However, quantitative estimates of climate after stochastic events remain relatively unexplored.
Impacts of large asteroids can generate global climate effects (e.g., Brugger et al., 2017;Pierazzo et al., 1998;Artemieva et al., 2017;Turbet et al., 2020).Simulations of such impacts have estimated the masses of water vapor injected into the atmosphere by these events, but the climatic implications have not been investigated (Koeberl & Ivanov, 2019) or have been inconclusive due to uncertainties on the radiative effect of water vapor (Erickson et al., 2020).Nonetheless, the temporal coincidence between the youngest Paleoproterozoic glacial deposits and the Yarrabubba event has sparked interest in the role of impacts in climate evolution (Erickson et al., 2020).Further, given their large kinetic energy, impacts can melt the ice table covering the planet surface and vaporize rocks from the underlying crust.Therefore, among other effects, impacts can (1) increase the planet's surface temperature by hundreds of degrees for large distances from the impact site, and (2) deposit ejecta material onto the icy surface far from the impact site, affecting the planet's surface albedo.Supervolcanic eruptions are also suggested candidates to escape the Snowball Earth state.Such events are defined by either their magnitude (M ≥ 8, proportional to the erupted mass) or their ejecta volume (1, 000 km 3 , de Silva & Self, 2022), and are thousands of times larger than typical volcanic eruptions (e.g., Hansen et al., 1978).The largest known super-eruption is Toba (M = 9.1), which is recorded in the 74,000-year-old Youngest Toba Tuff (YTT, Rose & Chesner, 1990).The ash ejected by the event covered ∼ 40 million km 2 of Earth's surface at a depth greater than 5 mm (Costa et al., 2014).Toba-like super-eruptions are expected to happen 1 − 2 times per million years (Mason et al., 2004;Cisneros de León et al., 2022).Supervolcanic eruptions have significant global climate effects (e.g., Hansen et al., 1978).Large volumes of volcanic gases are injected into the atmosphere, such as (1) CO 2 , which acts as a classical greenhouse gas, (2) SO 2 , which acts as a solar-reflecting coolant, and (3) H 2 O, which can act as a greenhouse gas, but can also form radiative-active clouds.Further, eruptions deposit ash and dust on the icy surface, reducing the albedo in comparison to clean ice and hence increasing the amount of sunlight absorbed by the planet's surface.Although often cited as a possible cause of the onset of global glaciation (Macdonald & Wordsworth, 2017), there is no detailed exploration of the viability of this scenario.
In this study, we quantify the heating effects of stochastic events on a globally glaciated Earth, testing the potential of such events to generate a runaway melting process that would lead to the termination of the Snowball state.Through impact simulations and ash-/ejectadispersal modeling (Sections 2.2 & 2.3), we quantify the changes in albedo and surface temperature produced by such events.Using a 1-D Energy Balance Model (EBM; Section 2.1), we then analyze the radiative perturbation caused by such temperature and albedo changes, following the evolution of the post-event surface temperature profile and ice coverage.By varying the scale of our stochastic events (e.g., impactor size, magnitude of supervolcanic eruption), we assess the capabilities of such events to trigger global deglaciation (Section 3).We discuss the likelihood of occurrence of our event scales during the Neoproterozoic, and thus the likelihood of deglaciation by stochastic events, in Section 4.

Method
2.1 Climate modeling Budyko (1969) and Sellers (1969)'s independent seminal works on the Earth's climate established the usefulness and robustness of 1-D EBMs.In a latitudinal EBM, the prognostic surface temperature is zonally averaged, thus leaving only the meridional dimension.This gives a remarkably good first-order estimate to solving the climate of a rocky planet (North & Coakley Jr, 1979;Spiegel et al., 2008;Dressing et al., 2010).The 1D-nature of the planet is further exemplified by parameterizing the meridional atmospheric heat transport as temperature-driven diffusion.Because bi-stability between the temperate Earth and Snowball Earth is a radiative balance process that exists through the difference of albedo between ice and liquid water (Ghil, 1994), we use an EBM focused on radiative balance (Williams et al., 1998;Spiegel et al., 2008;Dressing et al., 2010) to model the evolution of snowball deglaciation.The corresponding temperate and snowball temperature profiles are shown in Supplementary Figure S2.As 1-D EBMs are not longitudinally defined, we developed a method to account for spatially-located events, correcting input calculations from impact simulations and estimations of ash-dispersion from eruptions (see Supplementary Text S4).
For our simulations, we assume a planet following a circular orbit at 1 AU around a Sun-like star, with a modern Earth solar constant.We fix the planetary obliquity to 0°and also assume zero eccentricity.Land and oceans are uniformly distributed across the planet such that each latitude cell consists of a 70:30 modern Earth ocean-to-land ratio.Heat capacity and albedo parameterizations between ice and liquid water are smooth transition functions taken from Williams and Kasting (1997) and Spiegel et al. (2008) (see Supplementary Text S4).To constrain our study within a simple framework, we neglect the effects of aerosols and ice-melting latent heat: potential caveats from our assumptions are discussed in Section 4.However, the neglected processes (e.g., dust enrichment of the stratosphere due to impacts, SO 2 production of volcanism) tend to cool the atmosphere, making deglaciation more difficult to achieve.Even in these conditions, a significant warming of the planet is hardly reachable as shown in Section 3. Nonetheless, we account for variable relative humidity (RH), cloud forcing (Abbot, 2014) and the potential addition of CO 2 released by volcanic eruptions.We also evaluate the influence of the location of stochastic events by modeling eruptions and impacts at the equator and at 40°N.The effect of impacts is modeled by modifying the initial temperature distribution as per hydrocode simulations (see Section 2.2).Eruptions do not modify the temperature distribution, but affect the surface albedo due to ash-dispersion.

Impact simulations
We use the shock physics code iSALE (Collins et al., 2016) to model asteroidal impacts on a Snowball Earth state and produce 1-D thermal profiles for the surface of the planet in the aftermath of the impact.We additionally analyze the spread of the impact ejecta curtain in order to estimate the reduction in surface albedo associated with covering the ice with ejecta material.
We initialize simulations based on estimates of Snowball Earth's surface state (Hoffman et al., 1998).A 5 km thick global ice sheet thus sits at the surface of the target planet, with a rocky crust, mantle, and iron core sitting below.An average surface temperature of 227 K is prescribed for the pre-impact temperature profile, matching the steady-state temperature found by the EBM in the absence of a stochastic climate forcing event.Impactors consist entirely of dunite, matching closely the thermodynamic behaviors of the chondritic material that are expected from such objects (Benz et al., 1989).We run simulations for impactors with diameters between 40 − 100 km to assess the effects of impact scale.The impact velocity is constant for all impacts at 1.5× mutual escape velocity (v esc ), representing the most likely impact velocity (Chyba, 1991;Le Feuvre & Wieczorek, 2011), and we consider only head-on collisions, which is an approximation accounting for the small mass of the impactors relative to the Earth (more informations about the impact setup are provided in Supplementary Text S2).
Simulations are run until the motion in the mantle region excavated by the impact settles out (i.e., the mantle has rebounded and subsequent oscillatory motions have terminated).A 1-D surface temperature profile is then extracted, and used as the initial perturbation of the EBM.Because we use a 1-D latitudinal EBM, the initial temperature at the impact site results from the average of the high-temperature impact site longitudes and the low-temperature unperturbed longitudes (see Supplementary Text S4 for more details).As a result, the initial temperature in the EBM (zonal mean temperature, 200 − 350 K) is lower than the temperature at the impact site (local temperature, > 1000 K).We find that this approach slightly facilitates deglaciation by marginally overestimating the surface albedo (< 1 %).
An albedo profile is similarly determined, accounting for the exposure of crustal material after melting of the ice layer and the deposition of material from the impact ejecta curtain on top the ice table outside of the impact crater.We assume that ejecta reaches distances up to 5× the crater radius (e.g., Richardson et al., 2005), corresponding to distances of 2250 km, 3850 km, 3900 km and 5000 km from the impact site for the 40 km, 60 km, 80 km and 100 km impactors, respectively.The albedo for ejecta-covered areas is set to 0.2, similar to that of cryoconite (Hotaling et al., 2021).A step-like transition is then assumed between dark ejecta and ice albedo (0.7).Finally, we re-scale the temperature increase with the surface area that is affected by the impact and the size of latitudinal band, assuming the impact energy is conserved during the zonal transport.The temperature rise is submitted as the input to the EBM.

Supervolcanic eruptions
We implement the formalism of Pyle (1989) to model the variation of volcanic ash thickness (d) with distance from the eruption center.The exponential relationship between d and the isopach area A iso is given by: where d 0 is the ash maximum thickness (units of m) and k is the rate of ash thinning with isopach area, also in m (Pyle, 1989).We find k using (d, A iso ) values of volcanic ash fall from Toba, reconstructed through several tens of thickness measurements of the YTT tephra deposit (Costa et al., 2014) and digitized with QGIS (QGIS Development Team, 2024).
The ash deposition region is divided into fully-(d > 0.014 mm) and partially-covered (0.002 < d < 0.014 mm), the latter known as dusty ice (Le Hir et al., 2010).Assuming an elliptical geometry for ash-dispersal, the ash deposition area is converted in latitude and longitude.For Toba, the fully-covered area extends for ±43 • in latitude and ±86 • in longitude; the dusty ice area spans a wider range between ±50 • in latitude and ±100 • in longitude.At visible wavelengths, coarse dust has an albedo of 0.18 (Warren, 1982).Thus, assuming a clean ice albedo of 0.7, we assign the dusty ice albedo as changing linearly between clean ice and dust (Supplementary Figure S3).Such values are consistent with measurements of bare ice mean albedo (Warren et al., 2002).Fine dust correspond to higher albedo values, thus making deglaciation harder (Flanner et al., 2021).

Results
The EBM is able to provide the temporal evolution of surface temperature for up to several years after the moment of climate forcing.We find, however, that temperatures tend to converge within around two years (e.g., Figure 1).Geospatial variation in temperature, due and impacts of varying size (panel E-H).For the eruptions, each panel corresponds to a different set of initial conditions, highlighting the differences that such effects produce under the same climate forcing event.For the impacts, all panels show the same cloud-free atmosphere with 400 ppm CO2 of (A).The global-mean temperature (solid lines), minimum and maximum temperatures (dashed lines), and the melting temperature of the surface ice (dotted lines) are all shown.to the regional nature of the stochastic events, can be well represented by recording global temperature minima, maxima, and mean averages (Figure 1).
We first analyzed the effects of a Toba-like eruption (Figure 1A-D) for a cloud-free atmosphere similar to the present-day Earth atmosphere, with variable water vapor, 1 bar of nitrogen and 400 ppm of CO 2 , and with the eruption happening at the equator (Figure 1A).Additional simulations then included cloud radiative effects (14 W m −2 , Abbot, 2014) and increased atmospheric CO 2 at 800 ppm (Figure 1B) in order to test the robustness of the result to initial atmospheric conditions, and to account for greenhouse gases released by the eruption.Finally, we repeated these set simulations with the climate forcing occurring at 40 • N (Figure 1C).In none of these simulations does a Toba-like eruption change surface temperatures sufficiently to escape the Snowball state.In the best-case scenario for deglaciation (i.e., accounting for the cloud-warming effect and with 800 ppm of CO 2 ), the maximum equatorial temperature is only 255 K after relaxation from forcing.Other effects were further tested, but were found to contribute minorly to results within realistic bounds, including: the addition of humidity due to volcanic moisture, weaker atmospheric circulation, different ice albedo, and greater effective stellar insolation (see Supplementary Text S5 and Figure S1).
In the case of our impactors, the largest (d = 100 km) generates a transient water belt (note this may not be a true water belt but rather a function of the EBM's spatial structure), with maximum extents of 4.5 • N and 4.5 • S, and lasting for around two months (Figure 1H).During this time, the surface temperature of the planet decreases due to thermal emission.Once the freezing point is reached, the change of albedo induces a rapid runaway re-glaciation of the transient water belt.We thus find that even the 100 km impactor cannot melt a surface large enough to break the ice albedo feedback, which prevents the Snowball Earth deglaciation.The lower resultant temperatures and reduced geospatial domain of our smaller impactors (d = 40 − 80 km) results in no production of a water belt and indeed limited melting of the planet's surface (Figure 1E-G), with the same result of no deglaciation.Impactors larger than 100 km will melt more ice and will generate greater temperature perturbations.However, such impacts would also melt part of the planet crust and/or mantle, for which there is little evidence in the geologic record, and such massive impactors are highly unlikely by the time of the Neoproterozoic (Section 4).A simple analytical calculation tends to show that the energy required to deglaciate could only be delivered by an impactor that is larger than 100 km (Text S1).
The equilibrium spatial temperature profiles (i.e., spatially resolved profiles for the last times shown in Figure 1A-H) can provide further insight into the climate forcing effects of our stochastic events (Figure 2).Without any stochastic forcing, we observe the usual latitudinal distribution of surface temperatures, with a maximum of ∼ 205 K at the equator and a minimum of ∼ 190 K at the poles.
For a Toba-like eruption, we find that deposition of the volcanic ash leads to greatest temperature rises near the deposition site itself (Figure 2A).Ash deposition affects the surface albedo on a limited area, locally increasing the absorbed flux.For this reason, considering the eruption at 40 • N (Figure 2A, red line) centers the peak of the temperature profile to this location, and lower the maximum temperature of about 5 K due to a weaker incoming flux at this latitude.Additionally, for both cases, heat diffusion (atmospheric circulation) allows a global increase of temperature.By accounting for the radiative effect of the clouds and by assuming 800 ppm of CO 2 , the extra warming is around 10 K, independent of latitude and of the eruption location.We note that the results shown here do not account for the ash injected into the stratosphere, which reduces the bond albedo of the planet and thus induces global cooling (e.g., Abbot & Pierrehumbert, 2010).We discuss such processes, as well as others that we have chosen to neglect in the presented results, and their effect on our conclusions in Section 4.
For impact events, the converged temperature profiles are warmer than the non-perturbed simulation (Figure 2B), due to ejecta affecting the surface albedo.The rise in temperature is a function of the area covered by ejecta: the greater the area covered, the more the albedo of the surface is affected.However, the maximum temperatures for all impact scenarios are colder than those found in scenarios forced by the supervolcanic eruptions.The injection of heat that is provided by the impactors is a relatively minor effect in comparison to the change in surface albedo of the ice that is brought about by the ejecta curtain.Evidence for this can also be found in the temporal evolution of temperatures after the impacts (Figure 1E-H).In the four cases, the large radiative imbalance at the impact site dissipates heat rapidly relative to the timescale of global temperature changes.Therefore, the colder converged temperatures of our impacts are due to the diminished ejecta coverage that they produce in comparison to the coverage of the supervolcanic ash (see Supplementary Figure S3).Interestingly, for 80 km and 100 km impactors, we find a cold patch in the converged temperature profile at the impact site.We determine this to be due to the transient water belt that forms after the impact: when the planet cools, this liquid water freezes at the impact site as clean ice with greater albedo than the ejecta-covered ice around the crater.

Limitations of the climate model
H 2 O and CO 2 are efficient greenhouse gases and should facilitate the warming and deglaciation of the ice sheet.However, we find that variation of the gas partial pressures, within sensible bounds for the climate of a Neoproterozoic Snowball Earth, cannot push our stochastic events into regimes where they cause global deglaciation (see Supplementary Text S5).
There are a variety of effects that are not included in the EBM that could be considered as important.We neglect the formation of stratospheric aerosols, which would act to cool the planet (Macdonald & Wordsworth, 2017).In order to minimize the degrees of freedom in the model, we do not account for ice thickness or the latent heat released by ice melting, which could be important in accounting for the time-scales over which the water belt induced by the 100-km impactor freezes.In much the same way, oceanic circulation will affect the creation of ice in the water belt (e.g., Yang et al., 2020;Brunetti et al., 2019).However, a commonality between these effects is that their inclusion would serve only to increase the challenge in producing global deglaciation, which our models already suggest is not achievable given the scales of events considered.We thus suggest that our simplified EBM is able to point out possible directions for further works using more complex models.

Assumptions on impacts and large eruptions
Our simulations suggest that impacts larger than 100 km are necessary, albeit not sufficient, to initiate deglaciation on the planet.However, dynamical models indicate that, outside of the early Solar System, such impacts are expected to occur less than once every billion years (Koeberl & Ivanov, 2019).Moreover, geologic evidence for such event should be present, since a large projectile should have globally spread an ejecta layer analogous to the Chicxulub impact event (Alvarez et al., 1980).So far, no impact signature caused by a ∼ 100-km impactor has been found on Earth (Gyollai et al., 2014;Peucker-Ehrenbrink et al., 2016); as per Allen et al. (2022), the largest impact crater found is caused by a 20 − 25 km impactor.
The results we present are based on a single impact scenario, comprising of a fixed stratigraphy and an asteroidal projectile hitting the surface at 17 km s −1 .We chose an ice thickness of 5 km following Erickson et al. (2020) as a reasonable upper limit for the actual size of glaciers covering the planet.In reality, the ice thickness will depend on the location, as suggested by the geologic record (e.g., McMechan, 2000) and climate modeling (Hyde et al., 2000), ranging from meters to 5 km.Variation in ice thickness will produce a mechanically different response to impact.However, given the scale of impactor that our models find would be necessary to cause global deglaciation (i.e., d > 100 km), the ice thickness is relatively shallow in comparison, and we thus expect that the consequences of this effect would be minor.Further, we did not consider impactors of alternative compositions, such as icy bodies, in our models.While short-period comets (impact velocity of 30 km s −1 , Chyba, 1991) have kinetic energies similar to that of the asteroidal impacts we consider, long-period comets (50 km s −1 , Chyba, 1991) might change the impact mechanics, but are less likely to hit the Earth (Weissman, 2006).Such events would leave behind a distinct geochemical signal from a rocky impactor analogue that, while not yet discovered, should similarly perhaps not be excluded for now.Further modeling of such stochastic events could be entertained by future works, although will likely find challenges of scale similar to our work (i.e., requiring extremely unlikely comet sizes to cause deglaciation).Finally, we assumed zero porosity as a worst case scenario: for a given size, non-porous material ensures the maximum kinetic energy available to melt ice and vaporize target rocks.Similar to the climate choices made, therefore, relaxing this assumption only serves to make deglaciation more challenging.
In modeling volcanic eruptions, results are sensitive to the assumed dust distribution and the absence of lava warming and water vapor.Based on Icelandic dust and ash (Dragosics et al., 2016), a thickness of 15 mm is considered as the threshold for ice to melt due to the decrease in albedo.A thicker layer would thermally insulate the ice below, preventing melting.Increasing such a threshold would expand the dust-covered area subject to melting, favoring deglaciation.Furthermore, the local temperature increase due to erupting lava is disregarded, since its effect is negligible when averaged longitudinally.Similarly, we do not consider a long-term (i.e., multi-year) eruption, which would force a higher temperature around the volcano for a longer time.Additionally, as seen for the Hunga Tonga-Hunga Ha'apai eruption, a substantial amount (up to 10 11 kg) of water vapor can be injected into the stratosphere where it can have multi-year long radiative effects (net warming of the surface, Millán et al., 2022).Finally, aerosols ejected up to the stratosphere could increase snowfall, thus covering low albedo dusty ice.

Possible shortening of snowball events
In this work, we consider a low CO 2 pressure as the initial climate state.However, during a snowball event the carbon cycle is altered due to the ice shield preventing CO 2 (partial) dissolution into the ocean and restricted silicate weathering on land.On modern Earth, CO 2 outgassing is about 47 bar/Gyr (Catling & Kasting, 2017) through a combination of continental and submarine volcanism.Although there is no consensus so far on the exact volcanic outgassing flux for snowball Earth, multiple studies suggest that the value was much lower than present day.Fischer and Aiuppa (2020) estimate the arc production around 5.7 bar/Gyr while Dutkiewicz et al. (2024) argue that the mid-ocean ridge outgassing was even lower.We found that an extra radiative forcing of 50 W/m 2 (in addition to 14 W/m 2 of forcing induced by the clouds) is required, coupled with the effect of a 100 km impactor, to deglaciate the planet.That is corresponding approximately to 0.08 bar of CO 2 according to Wolf et al. (2018).In comparison, without any stochastic event the required forcing to deglaciate is equal to 90 W/m 2 (in addition to 14 W/m 2 of forcing induced by the clouds).Using modern day Earth outgassing from Fischer and Aiuppa (2020) as boundary values, it takes between 1.7 Myr and 14 Myr to reach a CO 2 pressure for which a 100-km impactor could induce a deglaciation, in our framework.An intermediate outgassing value of 25 bar/Gyr gives approximately 3.2 Myr to reach the required CO 2 pressure.These estimates should be considered as lower limits because we neglect the potential sinks of CO 2 existing on snowball Earth, as well as cooling processes of impacts (e.g., increase of the bond albedo due to aerosols).However, this timescale is shorter than the duration time of the Marinoan and Sturtian glaciations.Although a large impact during this period seems highly unlikely, the possibility of deglaciating an extra-solar planet in this way cannot be ruled out.

Conclusions
Earth's periods of global glaciation are traditionally conceived to have ended through a gradual build-up greenhouse gases in the atmosphere (e.g., CO 2 ).However, the shortest of these so-called Snowball Earth states were geologically short and stochastic events such as impacts and supervolcanic eruptions has been suggested as plausible deglaciation precesses.
Here, we present modeling on the response of Snowball Earth to such events under differing scenarios.We use impact simulations and ash spreading models to estimate the initial effect of each stochastic event, and follow these up by evolving Earth's surface environment through time using a 1-D EBM.
In our modeling framework, surface albedo of the planet is the dominant driving force of long-term temperature change in the aftermath of stochastic events.For impacts, this takes the form of the ejecta curtain stemming from the impact site.For eruptions, volcanic ash is transported over a wide range of latitude and longitudes.We thus observe that eruptions produce a greater change in Earth's surface albedo and hence a greater temperature response for typical scale events.However, even an impactor radius of 100 km is not sufficient to induce deglaciation of Snowball Earth.The warming at the impact site melts ice locally in the short-term, but due to a large radiative local imbalance, ice reforms in less than a year, returning the planet to the snowball state.Even if larger impactors could induce deglaciation, they are dynamically unlikely to have occurred on Earth at the time of its snowball episodes.Similarly, no recorded supervolcanic eruption could have single-handedly deglaciated Earth.The surface temperature increase due to dust deposition from the largest recorded event, Toba, does not reach the melting temperature of ice anywhere on the planet.
Our model thus places constraints on the magnitude of stochastic events that are incapable of triggering Snowball Earth termination.However, these constraints should be refined with more comprehensive modeling.For impacts, 3-D modeling can better resolve the ejecta distribution through accounting for oblique impacts.Impact simulations should also be able to account for the ice erosion effect due to the tsunami created upon impact, which would increase the low-albedo area.For the climate response, 3-D GCM modeling is required to include longitudinal and vertical resolutions, as well as to provide better representations of clouds, rain/snow, atmospheric dynamics, ocean dynamics, and distribution of aerosols ejected by the impactor (e.g., Abbot et al., 2011;Ashkenazy et al., 2013;Turbet et al., 2020) or released by the volcano (McGraw et al., 2024).For instance, even 2-D models could show more complex transitions between two steady states if specific eigenmodes are triggered by an impact happening at the right location (Mulder et al., 2021;Bastiaansen et al., 2022).
Lastly, in this work, we treated impacts and supervolcanoes as independent and instantaneous events all starting from the same initial state.Future analyses should investigate some of the prolonged effects associated with our stochastic events (e.g., multi-year lava eruptions).Additionally, the combined effect of impacts and eruptions should be considered, with impact-induced volcanism having been suggested as a potential supplement for the K-T extinction (e.g., Renne et al., 2015), although still debated (e.g., Bhandari et al., 1995).Multiple such events would likely stack their temperature responses, although unlikely to be in a linear manner, and could thus instigate a deglaciation of Snowball Earth.

Open Research
At present, the iSALE code (Collins et al., 2016) is not fully open-source.It is distributed via a private GitHub repository on a case-by-case basis to academic users in the impact community, strictly for non-commercial use.
The temperature profile after impact, used to perform climate simulations can be found from Chaverot et al. (2024).

Introduction
The following Supplementary Information provides first an analytical calculation of the order of magnitude of energy required to deglaciate Snowball Earth (Text S1).details about numerical setups of both the impacts simulations (Text S2, Table S2) and EBM simulations (Text S4, Table S3).Also, we provide additional discussion on the approximations made for impact simulations (Text S3), on the sensitivity texts done on the EBM setup (Text S5), as well as a discussion concerning implications of this work for exoplanets (Text S6).

Text S1: Analytical estimation of energy required for deglaciation
It is possible to estimate the order of magnitude of energy required to deglaciate the snowball Earth, using a simple analytical calculation.To do so, we first estimate equilibrium climates with a 0-D energy balance equation: where A is global mean albedo (reflectivity), F is global mean incoming solar flux, T s is global mean surface temperature, and a and b are constants.The left handside represents absorbed solar insolation; the right handside represents emitted infrared radiation.
July 12, 2024, 5:54pm Following Koll and Cronin (2018) and Wagner and Eisenman (2015), the parameters are , where T m = 273 K is the melting point of ice, a = −339.647W/m 2 , and b = 2.218 W/m 2 /K.Thus, 3 possible scenarios exist based on the choice of incoming solar flux F .When F is too small or too large, only one solution for climate can be found.When F is moderate, three solutions can be found.
The snowball bifurcation is represented by the warmest and coldest solutions (the middle solution T m is unstable).We find the range of F that allows the snowball bifurcation is Then we estimate the energy required for snowball deglaciation by picking up a typical solar forcing within the bifurcation range: The two stable solution can be solved analytically: To transition from the snowball state T snowball to T warm , we only need to heat the planet over the unstable point T m .Once T > T m , the planet will converge to T warm within finite time.The energy required for snowball deglaciation is: July 12, 2024, 5:54pm The iSALE shock physics code (Wünnemann et al., 2008) is an extension of the SALE hydrocode (Amsden et al., 1980).Modifications to the original SALE code include an elasto-plastic constitutive model, fragmentation models, various equations of state, and multiple materials (Melosh et al., 1992;Ivanov et al., 1997).More recent improvements include a modified strength model (Collins et al., 2004), a porosity-compaction model (Wünnemann et al., 2006;Collins et al., 2011), and a dilatancy model (Collins, 2014).
The thermodynamic behavior of each material is described by equations of state generated by the ANalytic Equation Of State program (ANEOS; Thompson & Lauson, 1974).The detailed properties of each material and the relevant equations of state are detailed in Table S1.
Our impactors are sufficiently large that the curvature of the Earth is important.However, for the information we wish to extract from the simulations, it is still appropriate for them to remain 2-D, which assists in the trade-off between resolution and computational resources.
Text S3: Discussion on impact modeling approximations Because of the strong effect of albedo on the climate evolution, our conclusions are sensitive to the assumptions made in estimating the dust distribution, namely assuming vertical impacts and no atmosphere present in the impact model.With a 2-D hydrocode, we are unable to reproduce oblique collisions, which affects (1) the amount of vaporized material, since shock waves would propagate in a larger volume of ice/crust in oblique impacts compared to vertical ones with the same kinetic energy, and (2) the shape of ejecta deposits, with oblique impacts preferentially distributing in the downrange direc-July 12, 2024, 5:54pm tion (Collins et al., 2020).Three-dimensional simulations are needed to better understand whether oblique impacts can cause a larger dust coverage.However, we again suggest that conclusions surrounding deglaciation would remain the same well within order of magnitude.In terms of not including an atmosphere in our impact simulations, the main effect is that ejecta follows ballistic trajectories that are unperturbed by atmospheric circulation.
Rock vapor is thus deposited on the surface without residing in the atmosphere.Consequently, the extension of dust we estimate is a lower limit, since we neglect spreading at further distances due to atmospheric transport (e.g., Kring, 2007).This would thus serve to decrease the planetary albedo, and should not affect our conclusions.

Text S4: EBM parameterization
We investigate the climate of the model planets using a 1-D time-dependent EBM introduced in Williams, Kasting, and Frakes (1998), Spiegel, Menou, and Scharf (2008) and further explored by Dressing, Spiegel, Scharf, Menou, and Raymond (2010).The equation governing the temporal evolution of the surface temperature T is given by: where x = sin λ gives the location of the grid point of latitude λ, t is time, C is the effective heat capacity of the atmosphere over this grid point, D = 0.5394 W m −2 K −1 is the diffusion coefficient which determines the efficiency of the latitudinal redistribution of heat, I(T, t) is the temperature-dependent infrared emission, S is the incoming stellar radiation By definition, the EBM is longitudinally-averaged which implies that modifying one grid cell affects all corresponding longitudes.As the considered stochastic events are limited to a particular region in space, we add a correction factor so that not all longitudes are impacted.More specifically, during the initialization of the simulation, we average the values of the albedo and of the temperature after the event with the non-impacted part of the planet at every latitude, accounting for the area of this event.By doing so, the initial maximum temperature of the EBM is much lower than the maximum temperature after the impact obtained from iSALE.

Text S5: Sensitivity tests
We perform several sensitivity tests in order to benchmark the original parameterization we introduced in the EBM for the thermal emission and to quantify the intrinsic effect of July 12, 2024, 5:54pm specific parameters.An exhaustive overview of these tests is presented in Table S3 and an overview of the most important ones is shown in Fig. S1.
First, the temperature profile we obtain using our parameterization of the thermal emission with 400 ppm of CO 2 is consistent with that of Dressing et al. (2010).Second, we tested the impact of several parameters of the EBM.The value of the diffusion factor (D in Eq. ( 8)) we used is taken from Dressing et al. (2010), but as mentioned in Pierrehumbert et al. (2011), cold climates have a reduced heat transport efficiency.Therefore, we tested D = 0.25 W m −2 K −1 as per Pierrehumbert et al. (2011) (see orange curve in Fig. S1).
The modification of the heat transport increases the temperature contrast between the equator and poles by a few Kelvin.We also added a cloud forcing equal to 14 W m −2 to account for their warming effect, according to Abbot (2014).This increases the global temperature by 10 K (difference between blue and red curves in Fig. S1).Due to our longitudinal correction, we overestimate the albedo by less than 1% depending on the albedo difference and on the temperature gradient along the longitude.

Text S6: Implications for Earth-like exoplanets
While the modeling in this study is predicated on Earth-like conditions, our conclusions may yet be useful for considerations of Earth-like exoplanets experiencing similar July 12, 2024, 5:54pm phenomena to Snowball Earth.For example, the question of planetary bi-stability within exoplanet demographics may be of particular relevance: if deglaciation through stochastic events is commonplace, the non-glacial planetary state will be heavily favored.
There are a variety of model components that may vary between worlds, such as stratigraphy, volcanic dust spreading through atmospheric circulation, atmospheric chemical composition and thermal structure, geologic record for impact and volcanic events, tectonic regime, and impactor size-frequency distributions.For instance, atmospheric dynamics and structure on planets with thicker or more turbulent atmospheres may produce ash spreading far beyond that found for Toba, enabling greater albedo change.Similarly, the thermo-mechanical response of basalt is different from that of granite (which we assumed to be the material underneath the ice layer), meaning that alternate crustal compositions and structures may produce wider dust spreading and/or surface temperature profiles.Finally, and perhaps most importantly, limiting factors on the scales of events for Earth do not necessarily apply to other planets.Supervolanic eruptions of greater magnitude than Toba may be more common, as may be projectiles of 100 km diameter or more.
Therefore, while we suggest that stochastic events are unlikely to have caused deglaciation of Snowball Earth, we do not rule out such events may be able to cause deglaciation on worlds where the parameters permit it.July 12, 2024, 5:54pm July 12, 2024, 5:54pm

Figure 1 .
Figure 1.Temporal evolution of surface temperature after a Toba-like eruption (panel A-D)

Figure 2 .
Figure 2. Latitudinal temperature profiles after (A) Toba-like eruptions, both at the equator (blue lines) and at 40 • N (red lines), and with both the cloud-free atmosphere with 400 ppm CO2 (solid lines) and an atmosphere including cloud radiative effects and 800 ppm CO2, and (B) impacts with impact diameters of 40 km (green), 60 km (blue), 80 km (orange), and 100 km (red), all striking at the equator and under cloud-free atmosphere with 400 ppm CO2 (solid lines).The red dashed line in panel B corresponds to a 100 km impactor under a cloudy atmosphere, assuming 400 ppm of CO2.The initial temperature profile (i.e., before climate forcing stochastic events) is shown for both types of events (dotted lines).
and A(T, t) is the temperature-dependent albedo of the grid point.The net radiative energy flux at a grid point is determined by the relationship between the irradiation July 12, 2024, 5:54pm Forget, Charnay, Wordsworth, and Pottier (2013), RH=0.45 is able to reproduce the infrared emission of a temperate Earth from a Global Climate Model.Higher RH values are used to represent the extra moisture accumulated in the atmosphere after a large eruption.This original parameterization of the thermal emission provides similar results to those obtained using the analytical parameterization from Dressing et al. (2010).The temperature profiles of the temperate and snowball regimes are shown in Fig. S2 (green and blue lines, respectively).Additionally, a stable water-belt configuration is shown for a diffusion coefficient D = 0.25 W m −2 K −1 (orange line).These different equilibrium states are obtained for a same setup by assuming different initial temperatures, highlighting climate multistability.

Figure
Figure S1 also shows that by accounting for 800 ppm of CO 2 instead of 400 ppm,

Figure S1 .
Figure S1.Latitudinal temperature without stochastic events for different parameterizations.

Figure S3 .
Figure S3.Latitudinal albedo parameterized for Toba-like supervolcanic eruptions with the eruption center at the equator (solid curve) and at 40 o N (dashed curve).The albedo values of ice

Table S1 .
Connection between impact energy Q and impact size D

Table S3 .
Overview of the EBM setups