Snowball Bistability Vanishes at Moderate Orbital Eccentricity

TL;DR

This study investigates whether Snowball climate bistability persists on exoplanets as orbital eccentricity increases. Using ExoCAM simulations of an aquaplanet with a 50 m slab ocean and a complementary low-order ice-thermodynamic model, the authors show that bistability vanishes at moderate eccentricities ($e \sim 0.25$–$0.3$) because the Snowball-to-Waterbelt transition flux $S^*_{ ext{SB}\rightarrow\text{WB}}$ decreases with seasonality, while the Waterbelt-to-Snowball flux $S^*_{ ext{WB}\rightarrow\text{SB}}$ remains roughly constant. The mechanism is attributed to ice self-insulation: melting occurs at the ice surface during summer, whereas freezing occurs at the ocean-ice interface at the bottom, making Snowball deglaciation harder under heightened seasonality. The ice-thermodynamic model reproduces the GCM results when key parameters ($\alpha_i$, $Q_{adv}$) are tuned, and robustness tests show the result holds across changes in obliquity, albedo, clouds, mixed-layer depth, and CO$_2$ forcing, with implications for the Rare Earth hypothesis and exoplanet habitability.

Abstract

Snowball episodes are associated with increases in atmospheric oxygen and the complexity of life on Earth, and they may be essential for the development of complex life on exoplanets. Sustained, stable Snowball episodes require a Snowball bifurcation and climate bistability between the globally ice-covered Snowball state and a state with at least some open ocean. We find that climate bistability disappears for an aquaplanet with a slab ocean in the global climate model ExoCAM when the orbital eccentricity is increased to 0.2-0.3. This happens because the Snowball state ceases to exist as seasonal insolation variations intensify, while the warm state remains stable due to the ocean's large heat capacity. We use a low-order ice-thermodynamic model to show that the Snowball state ceases to exist as seasonality increases because winter freezing at the ice bottom is reduced relative to summer melting at the ice top due to ice self-insulation. Combined with previous research showing that Snowball climate bistability diminishes for planets orbiting low-mass stars and ones with longer rotation periods, and that it disappears entirely for tidally locked planets, our work suggests that the Snowball climate bistability may not be as robust to planetary parameters as previously thought, representing one aspect of habitability more consistent with the rare Earth hypothesis than the Copernican principle.

Figures (9)

• Figure 1: Eccentricity versus annual-mean stellar flux for rocky exoplanets. Dots represent measured eccentricities, with vertical lines showing observational uncertainties; triangles denote planets for which only upper-limit constraints on eccentricity are available. Red symbols highlight planets whose upper eccentricity limits exceed 0.25, while black symbols represent the rest. The annual-mean stellar flux is computed using the mean eccentricity via $S \propto (a^2 \sqrt{1 - e^2})^{-1}$. Only potential rocky exoplanets with constrained eccentricity measurements are included in the sample. 72 out of 216 rocky planets have eccentricity upper limits exceeding 0.25. The sample is drawn from the confirmed exoplanet catalog on the NASA Exoplanet Archive nasa_exoplanet_archive_confirmed_2019: for transiting planets, the criteria for selecting "rocky planets" follow ji_cosmic_2025 (filled symbols); for planets without radius measurements, those with planetary masses less than 6 Earth masses are shown (open symbols).
• Figure 2: Snowball bifurcation diagrams for different eccentricities with the GCM ExoCAM. Blue circles represent simulations initialized with global ice coverage and red circles represent simulations initialized with no ice coverage. Seasonal variations in monthly-mean values are shown by vertical bars but are smaller than the circles in many cases. The gray shaded region indicates the range of stellar flux where Snowball bistability exists. The panels on the left ((a)-(f)) show the global-mean ice fraction ($F_{\text{ice}}$) as a function of annual-mean stellar flux. The panels on the right ((g)-(l)) show the global-mean surface temperature as a function of annual-mean stellar flux.
• Figure 3: Comparison between equatorial ice thermodynamics model (top: panels a–c) and GCM output at the equator (bottom: panels d–f) near the limiting case where the Snowball state ceases to exist ($S/S_0$=0.87, e=0.2), demonstrating qualitatively similar behavior. The panels show $F(t)$ (Eq. (\ref{['eq:EBM']})) (a,d), the ice thickness (b,e), and the surface temperature (c,f). $F(t)>0$ between times $t_1$ and $t_2$, leading to melting. We define the integral of $F(t)$ over this period as $E_1$, an important quantity discussed in the text. Similarly, freezing occurs when $F(t)<0$, and we define the integral of $F(t)$ over this period as $E_2$. The red and blue curves in panel (f) show the daily maximum and minimum temperature, while the thick purple line is the diurnal mean.
• Figure 4: Maximum stellar flux allowing a Snowball state as a function of eccentricity. Red dots represent default GCM results, corresponding to the right boundary of the bistability region in Fig. \ref{['fig:bifur-diagram']}. The error bars indicate the range in annual-mean stellar flux spanned by the two simulations closest to this boundary. We also plot GCM sensitivity tests with the ice albedo increased by a factor of 1.04 (blue dots) and the obliquity increased to 23 degrees (purple dots). The red line shows results from the ice-thermodynamic model with $\alpha_i$ and $Q_{adv}$ tuned to default GCM results at low eccentricity. The blue and purple lines show results of the ice-thermodynamic model with parameters adjusted for the GCM sensitivity tests. The horizontal gray line shows the approximate critical stellar flux for the transition from Waterbelt to Snowball in the GCM, corresponding to the left boundary of the bistability region.
• Figure 5: GCM bifurcation diagram for Earth-like obliquity cases. The format follows Fig. \ref{['fig:bifur-diagram']}, where red markers represent Warm-Start simulations and blue markers represent Cold-Start simulations. Solid lines indicate simulations with Earth-like obliquity ($23^\circ$), while dotted lines show corresponding zero-obliquity cases for comparison. From top to bottom, eccentricity increases.