Ice Caps and Ice Belts: The Effects of Obliquity on Ice−Albedo Feedback

Obliquity has a profound effect on the meridional and seasonal distribution of insolation. Figure 2 shows the complete patterns of daily mean insolation for four representative values of obliquity: $\beta =0^\circ$ (the perpetual equinox), $\beta =23\buildrel{\circ}\over{.} 45$ (Earth's present-day value), $\beta =55^\circ$ (the critical value at which the annual mean equator-to-pole insolation gradient reverses), and $\beta =90^\circ$ (the subsolar point is at the North Pole at its summer solstice). In addition to the reversal of the annual mean gradient for $\beta \gt 55^\circ$, the seasonality of insolation is much stronger at high obliquity.

Zoom In Zoom Out Reset image size

In the following section, we describe an approximate series expansion for the daily average insolation and show the dependence of the expansion coefficients on the obliquity angle β. We do not treat the diurnal cycle since we use a zonally averaged model.

In pursuit of a non-dimensional formalism for the planetary albedo feedback problem, we begin by expanding the daily average insolation $S(x,t)$ in a Fourier–Legendre series (e.g., North & Coakley 1979):

$$\begin{eqnarray}&&S(x,t)={Qs}(x,t),\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&s(x,t)=\displaystyle \sum _{l=0,\,k=0}({a}_{{lk}}\cos k\omega t+{b}_{{lk}}\sin k\omega t){P}_{l}(x),\end{eqnarray} \tag{ 1b }$$

where Q is the global annual average insolation in W m⁻² (4Q is known as the solar constant for Earth), $s(x,t)$ is the normalized daily mean insolation (unit global, annual mean), P_l(x) is the lth-order Legendre polynomial, $\omega =2\pi /{t}_{\mathrm{year}}$, where t_year is the length of the year, and we use the independent variable $x=\sin \phi$, where ϕ is the latitude. As shown by North & Coakley (1979), all odd coefficients in this expansion aside from l = 1 vanish. For simplicity we will limit our analysis to circular orbits (zero eccentricity) for which the first harmonic is sufficient when l = 1, and the second harmonic is sufficient for l = 2. Therefore, we truncate the series to

$$\begin{eqnarray}&&s(x,t)=1+{s}_{11}\cos (\omega t){P}_{1}(x)+({s}_{20}+{s}_{22}\cos (2\omega t)){P}_{2}(x),\end{eqnarray} \tag{ 2 }$$

where ${P}_{1}(x)=x,{P}_{2}(x)=\tfrac{1}{2}(3{x}^{2}-1)$ are the first and second Legendre polynomials. Here we are choosing to set t = 0 at the northern hemisphere winter solstice.

For a planet in a circular orbit, s₂₀ (annual mean equator−pole insolation gradient), s₁₁ (amplitude of the annual cycle), and s₂₂ (amplitude of the semiannual cycle) are all simple functions⁴ of the obliquity angle β:

$$\begin{eqnarray}&&{s}_{11}=-2\sin \beta ,\end{eqnarray} \tag{ 3a }$$

$$\begin{eqnarray}&&{s}_{20}=-\displaystyle \frac{5}{16}(2-3{\sin }^{2}\beta ),\end{eqnarray} \tag{ 3b }$$

$$\begin{eqnarray}&&{s}_{22}=\displaystyle \frac{15}{16}{\sin }^{2}\beta .\end{eqnarray} \tag{ 3c }$$

These coefficients are plotted in Figure 3 for β between 0° and 90°. The distribution of the errors associated with expansion (2) is shown with dashed lines in Figure 2. The rms error of the global annual mean is less than 7% for any obliquity (lower panel of Figure 3). Most of this error lies in the spatial structure of the annual mean at obliquities near zero and in the annual cycle at higher obliquity, so improving this fit would require higher-order Legendre polynomials at low obliquity and higher harmonics of the annual cycle at high obliquity.

Zoom In Zoom Out Reset image size

Since the bulk of this paper focuses on the annual mean energy balance, we note here that the annual mean insolation is

$$\begin{eqnarray}&&\overline{s(x,t)}=1+{s}_{20}{P}_{2}(x),\end{eqnarray} \tag{ 4 }$$

where the overbar denotes an annual average. The coefficient s₂₀ is negative for low, Earth-like obliquity (insolation decreases poleward), reaches zero at the critical obliquity ${\beta }_{c}=\arcsin (\sqrt{2/3})\approx 55^\circ$, and is positive for high obliquity ($\beta \gt {\beta }_{c}$) for which the poles receive more sunlight annually than the equator. Since β can range between 0° and 90°, s₂₀ ranges between −5/8 and +5/16. We emphasize again that Equation 3(b) is a good approximation for the annual mean insolation for all but very low obliquities (for which s goes to zero at the poles).

We now adopt the diffusive EBM (North 1975b; North & Coakley 1979), which expresses the zonal mean energy budget for the climate system in terms of the zonal mean surface temperature $T(x,t)$. In its seasonally varying form, the model can be written as

$$\begin{eqnarray}&&C\displaystyle \frac{\partial T}{\partial t}={aQs}(x,t)-[A+{BT}]+\displaystyle \frac{K}{{R}^{2}}{{\rm{\nabla }}}^{2}T,\end{eqnarray} \tag{ 5 }$$

where for convenience we define a dimensionless meridional Laplacian operator on the sphere:

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\equiv \displaystyle \frac{\partial }{\partial x}\left[(1-{x}^{2})\displaystyle \frac{\partial }{\partial x}\right].\end{eqnarray} \tag{ 6 }$$

The lhs of Equation (5) is the seasonal heat storage, with C a column heat capacity in J m⁻² °C⁻¹. The first term on the rhs is the absorbed solar radiation, with a the co-albedo (the absorbed fraction of the incident solar radiation). The insolation function $s(x,t)$ depends on the obliquity angle β and length of year ${t}_{\mathrm{year}}=2\pi /\omega$ as described in Section 2. A + BT is a linear parameterization of the Outgoing Longwave Radiation (OLR), with dimensional parameters A and B governing its efficiency (dependent on atmospheric properties such as greenhouse gas concentration, cloudiness, lapse rate, etc.). The third term on the rhs represents the convergence of heat transport due to atmospheric and oceanic motions. The transport is parameterized as a diffusive flux down the large-scale temperature gradient, with the diffusivity parameter K setting its efficiency.⁵ Finally, R is the planetary radius.

Albedo feedback is introduced into the model by invoking a temperature dependence of the co-albedo, $a=a[T(x,t)]$. Following many classic studies (e.g., Budyko 1969; Held & Suarez 1974; North 1975b), we adopt a simple step function in which the ice and snow line is tied to a particular isotherm T₀:

$$\begin{eqnarray}a[T(x,t)]={a}_{\lrcorner }=\left\{\begin{array}{cc}{a}_{0}, & T(x,t)\gt {T}_{0}\\ {a}_{1}, & T(x,t)\lt {T}_{0}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where ${a}_{0}\gt {a}_{1}$. The temperature dependence of ${a}_{\lrcorner }$ introduces a nonlinearity into the EBM and raises the possibility of multiple equilibria and unstable ice growth. For the seasonal model, the natural choice for threshold temperature T₀ would be the relevant freezing point (0°C or about $-2^\circ {\rm{C}}$ for Earth's land and ocean surfaces, respectively). Note that the actual freezing point for seawater is salinity dependent and would be lower for a hypersaline ocean (as low as −21.1°C for a NaCl−water brine).

Many studies using an EBM to analyze features of Earth's climate have focused on annual mean conditions driven by the annual mean insolation $\overline{s(x,t)}$. Averaging the seasonal EBM, Equation (5), over a steady seasonal cycle yields

$$\begin{eqnarray}&&0=Q\overline{{as}}-[A+B\overline{T}]+\displaystyle \frac{K}{{R}^{2}}{{\rm{\nabla }}}^{2}\overline{T},\end{eqnarray} \tag{ 8 }$$

which is the starting point for many classic studies of Earth's energy balance (e.g., North 1975a, 1975b).

It is often supposed that seasonal covariance between insolation and albedo can be ignored or parameterized. We can define an "effective co-albedo" $\tilde{a}$ for the annual mean energy budget as the ratio

$$\begin{eqnarray}&&\tilde{a}=\displaystyle \frac{\overline{{as}}}{\overline{s}}.\end{eqnarray} \tag{ 9 }$$

The insolation term in Equation (8) is then replaced with $Q\tilde{a}\overline{s}$, with $\overline{s}$ given as a function of latitude and obliquity by Equation (4). The annual mean EBM is then closed with an assumption about the dependence of $\tilde{a}$ on the annual mean temperature $\overline{T}$.

Historically, the step-function form of Equation (7) has been applied to the annual mean model, with a threshold temperature T_f typically less than the instantaneous freezing temperature T₀. Specifically,

$$\begin{eqnarray}\tilde{a}={a}_{\lrcorner }=\left\{\begin{array}{cc}{a}_{0}, & \overline{T}(x)\gt {T}_{f}\\ {a}_{1}, & \overline{T}(x)\lt {T}_{f}\end{array}\right.,\end{eqnarray} \tag{ 10 }$$

where the canonical choice is ${T}_{f}=-10^\circ {\rm{C}}$ (e.g., North et al. 1981). Why is the relevant threshold −10°C rather than 0°C? This question has received relatively little attention in the EBM literature. The choice of threshold appears to go back to the seminal work of Budyko (1969), who noted that the transition to permanent year-round ice and snow cover in the Arctic occurs at a mean latitude of 72°N, and −10°C is roughly the observed annual mean temperature at this latitude.

Presumably, the transition to permanent ice and snow on a planet with very weak seasonal temperature variations would occur at a location with annual mean temperature closer to 0°C. We thus hypothesize that the magnitude of the difference ${T}_{0}\mbox{--}{T}_{f}$ is linked to the amplitude of the seasonal cycle. Although preliminary theoretical investigation of this relationship appeared promising, numerical solutions revealed that a step-function parameterization of the co-albedo is unlikely to work when the amplitude of the seasonal cycle is large. We thus defer deeper analysis of Equation (10) to future work. For now, we will simply note that our formal analysis (beginning below) relies on T_f or T₀ as a reference point for the non-dimensionalization.

For the seasonal model (5), a total of 11 dimensional parameters have been introduced: Q, a₀, a₁, A, B, K, R, β, T₀, t_year, and C. Similarly, the annual model (8) depends on nine parameters—the first eight in the seasonal list, plus T_f (which we hypothesize to be an implicit function of the seasonal parameters). For present-day Earth, many of these parameters are known or can be estimated from observations⁶ , and ${s}_{20}=-0.48$ gives a good approximation of the present-day annual mean insolation. Thus, the annual model is often presented as having only one free parameter, the heat diffusion constant K, which is tuned to reproduce the modern-day climate (e.g., North 1975a). However, these parameters can vary widely among exoplanets. It is therefore important to identify all of the key non-dimensional parameters in order to investigate the full range of possible climates embodied in the EBM.

We introduce the following dimensionless constants:

$$\begin{eqnarray}&&\tau =\omega t=2\pi t/{t}_{\mathrm{year}},\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{C\omega }{B},\end{eqnarray} \tag{ 11b }$$

$$\begin{eqnarray}&&\delta =\displaystyle \frac{K}{{R}^{2}B},\end{eqnarray} \tag{ 11c }$$

$$\begin{eqnarray}&&q=\displaystyle \frac{{a}_{0}Q}{A+{{BT}}_{\mathrm{ref}}},\end{eqnarray} \tag{ 11d }$$

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{a}_{0}-{a}_{1}}{{a}_{0}},\end{eqnarray} \tag{ 11e }$$

and non-dimensionalize the surface temperature and OLR with

$$\begin{eqnarray}&&{T}^{* }(x)=\displaystyle \frac{A+{BT}(x)}{A+{{BT}}_{\mathrm{ref}}}.\end{eqnarray} \tag{ 12 }$$

The reference temperature T_ref in Equations 11(d) and (12) is taken to be the temperature threshold—either T₀ for the seasonal model or T_f for the annual model, so that in either case ${T}^{* }=1$ at the ice edge.

The seasonal EBM can then be written as

$$\begin{eqnarray}\gamma \displaystyle \frac{\partial {T}^{* }}{\partial \tau }-\delta {{\rm{\nabla }}}^{2}{T}^{* }+{T}^{* }={qs}(x,\tau )\left\{\begin{array}{cc}1, & {T}^{* }\gt 1\\ 1-\alpha , & {T}^{* }\lt 1\end{array}.\right.\end{eqnarray} \tag{ 13 }$$

Substituting in the insolation expansion (2) and using Equation (3), we conclude that the seasonal model has five apparently independent parameters: $\gamma ,\delta ,q,\beta ,\alpha$. These five are usually taken to be uniform in latitude.

The annual EBM can be written as

$$\begin{eqnarray}\delta {{\rm{\nabla }}}^{2}{T}^{* }-{T}^{* }=-q[1+{s}_{20}{P}_{2}(x)]\left\{\begin{array}{cc}1, & {T}^{* }\gt 1\\ (1-\alpha ), & {T}^{* }\lt 1\end{array}.\right.\end{eqnarray} \tag{ 14 }$$

The annual model has just four independent parameters: $\delta ,\alpha ,q,{s}_{20}$ (with s₂₀ uniquely determined by obliquity β). However, the definitions of T* and q are not identical in Equations (13) and (14) because the reference temperatures may differ.

Here we describe the parameters defined in Equation (11) in physical terms. We provide typical Earth-like values for each parameter and note some factors governing their plausible ranges for habitable exoplanets. The Earth values are derived from dimensional parameters in North (1975b). More accurate or up-to-date values are certainly possible, but our main interest here is in situating the non-dimensional model within the existing EBM literature.

• T*: Proportional to both temperature and OLR. ${T}^{* }=1$ by definition at the ice edge.
• $\gamma :$ Seasonal heat capacity of the system relative to the radiative decay of temperature over one year. γ decreases with length of year and radiative damping, and increases with fractional ocean coverage and efficiency of ocean mixing. With the dimensional parameters from North (1975b), we have γ ≈ (0.55 m⁻¹)H, where H is the depth of water. $\gamma \approx 1$ for a dry Earth whose heat capacity is dominated by land and atmosphere. A realistic value for Earth is in between 5 and 20, with larger values appropriate for the ocean-dominated southern hemisphere. To apply a uniform gamma requires a compromise across this range.
• $\delta :$ Efficiency of dynamical heat transport. δ measures the relative importance of transport versus local radiative damping in response to a localized heat source (Stone 1978). Transport will smooth out meridional temperature variations over length scales smaller than $\sqrt{\delta }R$ (Lindzen & Farrell 1977). δ depends on atmospheric properties such as mass, greenhouse gas levels and cloud cover, and dynamical factors such as rotation rate and planetary radius (Williams & Kasting 1997; Vallis & Farneti 2009). Other factors might include the topographic forcing of atmospheric stationary waves (e.g., Cook & Held 1988) and the temperature dependence of latent heat transport (e.g., Caballero & Langen 2005). The role of oceans on the effective value of δ is complex due to the multiple spatial scales of ocean heat transport and their tight coupling to the ice extent (e.g., Rose & Marshall 2009; Rose et al. 2013; Ferreira et al. 2014; Rose 2015). Ocean heat transport has been found to increase with rotation rate in uncoupled simulations (Cullum et al. 2014) but decrease in coupled atmosphere−ocean simulations (Vallis & Farneti 2009). For present-day Earth, we take $\delta \approx 0.31$ following North (1975b). δ could vary widely across different planets.
• ${\text{}}q:$ Radiative forcing. The numerator in Equation 11(d) is the global mean absorbed shortwave for an ice-free planet. The denominator is OLR at the ice edge. q is sensitive to solar irradiance (Q), greenhouse gas amount (A), and planetary albedo, as well as choice of reference temperature. q could of course range widely across planets with different orbital distances and stellar output. The region of interest is the habitable zone near $q\approx 1$, since planets at large or small q would be locked into very warm or Snowball climates, respectively. For Earth, $q\approx 1.2$ in the annual mean model and 1.1 in the seasonal model.
• $\alpha :$ A measure of the potential ice−albedo feedback. α is bounded between 0 (no feedback; no change in albedo across the ice edge) and 1 (extreme feedback). α will be reduced for a cloudier planet and vice versa. α is also systematically reduced for cooler host stars with longer emission wavelengths (Shields et al. 2013). For Earth, North (1975b) used values corresponding to $\alpha =0.44$.

Before presenting our analysis of the stability of the annual mean model (14) at high and low obliquity, we take a brief detour into the effects of the seasonal cycle in the time-dependent seasonal EBM. An explicit goal of this analysis is to determine the non-dimensional parameter regime in which the seasonal effects are small, which then justifies the use of the annual mean model.

We begin by removing the nonlinearity from Equation (13) by taking $\alpha =0$ (uniform, constant albedo). This simplification makes Equation (13) a linear PDE for T*, and allows us to obtain periodic seasonal solutions of the form

$$\begin{eqnarray}{T}^{* } & = & \overline{{T}^{* }(x)}+{T}_{11}^{* }{P}_{1}(x)\cos (\tau -{{\rm{\Phi }}}_{11})\\ & & +{T}_{22}^{* }{P}_{2}(x)\cos (2\tau -{{\rm{\Phi }}}_{22}),\end{eqnarray} \tag{ 15 }$$

where $\overline{{T}^{* }(x)}$ is the annual mean temperature, ${T}_{11}^{* }$ and ${{\rm{\Phi }}}_{11}$ are the amplitude and phase lag (relative to insolation) of the annual cycle, and ${T}_{22}^{* }$ and ${{\rm{\Phi }}}_{22}$ are the same quantities for the semiannual cycle. We non-dimensionalize ${T}_{11}^{* }={{BT}}_{11}/(A+{{BT}}_{0})$ and ${T}_{22}^{* }={{BT}}_{22}/(A+{{BT}}_{0})$ such that they have zero annual mean.

Solutions for $\overline{{T}^{* }(x)}$ are discussed later in Section 4.1. In Appendix A, we derive solutions for the seasonal departures. For the annual cycle, we find

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{11}=\arctan \left(\displaystyle \frac{\gamma }{1+2\delta }\right),\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&{T}_{11}^{* }=\displaystyle \frac{{{qs}}_{11}}{\sqrt{{(1+2\delta )}^{2}+{\gamma }^{2}}},\end{eqnarray} \tag{ 16b }$$

which implies that temperature is coherent across latitudes in phase and increases in amplitude linearly with $| x|$ toward the pole in each hemisphere. The solutions for the semiannual component are

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{22}=\arctan \left(\displaystyle \frac{2\gamma }{1+6\delta }\right),\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{T}_{22}^{* }=\displaystyle \frac{{{qs}}_{22}}{\sqrt{{(1+6\delta )}^{2}+{(2\gamma )}^{2}}}.\end{eqnarray} \tag{ 17b }$$

The above solutions show the relative roles for local heat storage and meridional heat transport in the seasonal cycle. Setting $\delta =0$ gives a local radiative equilibrium solution in the absence of transport and emphasizes the role of γ in both damping the temperature seasonal amplitude and shifting its phase relative to the insolation.

The deep-water limit is $\gamma \gg 1+2\delta$, appropriate for a planet with a deep mixed ocean layer or short solar year. In this limit, the phase shifts for both annual and semiannual components are nearly $\pi /2$ and the temperature seasonal amplitude is weak. Present-day Earth is much closer to this limit than to the shallow-water limit to be discussed below, as evidenced by the observed phase shift of our seasons. To leading order in ${\gamma }^{-1}$, the deep-water limits for the annual and semiannual components are

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{11}\approx \displaystyle \frac{\pi }{2}-\displaystyle \frac{1+2\delta }{\gamma }\qquad {T}_{11}^{* }\approx \displaystyle \frac{{{qs}}_{11}}{\gamma },\end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{22}\approx \displaystyle \frac{\pi }{2}-\displaystyle \frac{1+6\delta }{2\gamma }\qquad {T}_{22}^{* }\approx \displaystyle \frac{{{qs}}_{22}}{2\gamma }.\end{eqnarray} \tag{ 18b }$$

It is notable that the (weak) seasonal amplitude becomes independent of δ in this limit.

Conversely, the shallow-water limit $\gamma \ll 1+2\delta$ is appropriate for a dry planet with a very long solar year. In this case, the temperature has large seasonal amplitude nearly in phase with the Sun. Specifically we get, to leading order,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{11}\approx \displaystyle \frac{\gamma }{1+2\delta }\qquad {T}_{11}^{* }\approx \displaystyle \frac{{{qs}}_{11}}{1+2\delta },\end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{22}\approx \displaystyle \frac{2\gamma }{1+6\delta }\qquad {T}_{22}^{* }\approx \displaystyle \frac{{{qs}}_{22}}{1+6\delta }.\end{eqnarray} \tag{ 19b }$$

It is notable that even in this limit the seasonal amplitude is substantially smaller than the local radiative equilibrium value ($\delta =0$) due to transport from other latitudes. Physically, the amplitude of the annual cycle of insolation increases poleward, but this is offset by the poleward increase of the seasonally varying convergence of heat transport from lower latitudes. There is a tradeoff between the damping of the seasonality due to local heat storage versus that due to meridional transport. This damping is stronger for the semiannual than for the annual component, which contributes to the smallness of the semiannual component away from the equator.

In the next section, we will present our analysis of the nonlinear annual mean model with albedo feedback. The annual model is independent of the heat capacity parameter γ, and through Equation (10) assumes a sharp boundary between the ice-covered and ice-free regions. The relevance of these results in the presence of substantial seasonal migration of the snow and ice line is questionable. We therefore limit ourselves to consideration of the deep-water regime $\gamma \gg 1+2\delta$ in which those migrations are small, noting that this restriction is more stringent for small, slowly rotating planets (for which δ is larger). We will revisit this question with numerical integrations of a seasonal model in Section 6.

Before considering partially ice-covered planets, we first consider solutions with constant albedo across the whole domain. In this case, Equation (14) is linear and simple exact solutions can be written:

$$\begin{eqnarray}&&{T}^{* }=q\left(1+\tfrac{{s}_{20}}{1+6\delta }{P}_{2}(x)\right)\end{eqnarray} \tag{ 20a }$$

for ice-free conditions, and

$$\begin{eqnarray}&&{T}^{* }=(1-\alpha )q\left(1+\tfrac{{s}_{20}}{1+6\delta }{P}_{2}(x)\right)\end{eqnarray} \tag{ 20b }$$

for Snowball conditions.

An ice-free planet requires the coldest temperature anywhere to be at or above T_f, so ${T}^{* }\geqslant 1$ everywhere. The coldest temperature occurs at the pole (x = 1) for low obliquity but at the equator (x = 0) for high obliquity. We define q_free as the minimum radiative forcing for which an ice-free planet can exist. Noting that ${P}_{2}(1)=1$ and ${P}_{2}(0)=-1/2$, the condition is

$$\begin{eqnarray}{q}_{\mathrm{free}}^{-1}=\left\{\begin{array}{lc}1+\tfrac{{s}_{20}}{1+6\delta }, & {s}_{20}\lt 0\\ 1-\tfrac{{s}_{20}}{2(1+6\delta )}, & {s}_{20}\gt 0\end{array}\right.,\end{eqnarray} \tag{ 21 }$$

where we remind the reader that s₂₀ ranges between −5/8 (for $\beta =0^\circ$) and +5/16 (for $\beta =90^\circ$) (Figure 3).

Notice that this says that ice-free conditions can exist for weaker radiative forcing on high-obliquity planets compared to low-obliquity planets with the same heat transport efficiency, both because $| {s}_{20}|$ is small (even at extreme $\beta =90^\circ$ it is still smaller than for present-day $\beta =23\buildrel{\circ}\over{.} 5$), and because of the extra factor of 1/2 when ${s}_{20}\gt 0$ ($\beta \gt 55^\circ$). Equivalently, warm climates can exist with weaker heat transport efficiency on high-obliquity planets compared to low-obliquity planets at the same radiative forcing. Physically, this arises because the annual mean insolation is larger at the equator for high obliquity than it is at the poles for low obliquity.

As a concrete example, if $\delta =0.31$ for an Earth-like planet, then condition (21) says that an ice-free climate is possible for $q\geqslant 1.2$ for $\beta =23.5^\circ$. Thus, present-day Earth is near a radiative forcing that supports an ice-free climate. On the other hand, for $\beta =90^\circ$, the condition is $q\geqslant 1.06$, meaning the planet could support an ice-free climate even under a 10% reduction in insolation (or a similar reduction in global mean temperature from reduced greenhouse gases). Figure 4 shows q_free plotted as a function of β and δ. The smallest radiative forcing for which ice-free conditions are possible occur at values of obliquity close to ${\beta }_{c}=55^\circ$. Here, the annual mean insolation has zero meridional gradient, and ${q}_{\mathrm{free}}\to 1$. Under these conditions, the surface is essentially isothermal and just above the freezing point.

Zoom In Zoom Out Reset image size

More generally, Figure 4 shows that q_free is substantially smaller for high obliquity than low obliquity over a wide range of transport parameters δ. All non-dimensional parameters being equal, a high-obliquity planet is more likely to be perennially ice free than a low-obliquity planet.

We can similarly define q_snow as the maximum radiative forcing for which a completely ice-covered climate can exist:

$$\begin{eqnarray}{q}_{\mathrm{snow}}^{-1}=(1-\alpha )\left(1+\displaystyle \frac{| {s}_{20}| }{1+6\delta }\left\{\begin{array}{lc}1/2, & {s}_{20}\lt 0\\ 1, & {s}_{20}\gt 0\end{array}\right\}\right).\end{eqnarray} \tag{ 22 }$$

It is well-known from EBM studies of Earth's climate that both climatic extremes (ice-free and Snowball) are possible over a wide range of radiative forcings, i.e., ${q}_{\mathrm{snow}}\gg {q}_{\mathrm{free}}$. The equilibrium climate is multiple valued and capable of undergoing hysteresis, which is central to the Neoproterozoic Snowball Earth hypothesis (e.g., Hoffman et al. 2017). How general is this result? What combination of parameters permits the coexistence of the ice-free and Snowball solutions? From Equations (21) and (22), both solutions are possible for a range of q values as long as the albedo parameter α is sufficiently large—that is, if the planet has a strong enough albedo feedback. Specifically, we can define a minimum value of α that permits both solutions for some range of radiative forcings (whether we are actually in the multiple equilibrium regime would then depend on q). The condition on α is

$$\begin{eqnarray}\alpha \geqslant \displaystyle \frac{3}{2}\left\{\begin{array}{lc}{\left[\tfrac{1}{2}-\tfrac{1+6\delta }{{s}_{20}}\right]}^{-1}, & {s}_{20}\lt 0\\ {\left[1+\tfrac{1+6\delta }{{s}_{20}}\right]}^{-1}, & {s}_{20}\gt 0.\end{array}\right.\end{eqnarray} \tag{ 23 }$$

Or, since ${s}_{20}/(1+6\delta )$ is a small number, to a crude approximation we can write

$$\begin{eqnarray}&&\alpha (1+6\delta )\gtrsim \displaystyle \frac{3}{2}| {s}_{20}| \end{eqnarray} \tag{ 24 }$$

for either sign of s₂₀. For the Earth-like parameter values cited above, we have $\alpha (1+6\delta )\approx 1.3$. This condition is then satisfied for any obliquity. We conclude that multiple equilibria are a rather common property of planets with albedo feedback, perhaps especially in the high-obliquity regime since $| {s}_{20}|$ is smaller.

When temperature crosses the ice threshold ${T}^{* }=1$ somewhere in the domain, Equation (14) becomes nonlinear. However, if we assume that the ice edge x_s is known, then the problem is piecewise linear in the warm and icy regions. We can then solve for the parameter values that would give the assumed ice edge. This is the basis for the analytical solutions given by North (1975b) and North et al. (1981). Here we extend these solutions to the high-obliquity case. Details are given in Appendix B, but we outline the solution method here.

Because the icy region exists on opposite sides of x_s in the low- and high-obliquity regimes, the appropriate form of the rhs of Equation (14) is slightly different depending on the sign of s₂₀. The upper form applies to ${s}_{20}\lt 0$ and $x\lt {x}_{s}$, or for ${s}_{20}\gt 0$ and $x\gt {x}_{s}$, while the lower form with the factor $(1-\alpha )$ applies to ${s}_{20}\lt 0$ and $x\gt {x}_{s}$, or for ${s}_{20}\gt 0$ and $x\lt {x}_{s}$.

These are second-order ODEs (forms of Legendre's equation) for ${T}^{* }(x)$. Since the annual mean model is symmetric about the equator, we solve on the hemisphere $x\in [0,1]$. We split this domain into two regions $x\lt {x}_{s}$, $x\gt {x}_{s}$, and solve a two-point boundary value problem in both regions. The particular solution is proportional to ${P}_{2}(x)$ following the insolation, and the general solution to the homogeneous equation is written in terms of the Legendre functions on both sides. Four boundary conditions are determined by conservation of energy: zero heat transport at the pole (x = 1) and the equator (x = 0), and continuous transport at x_s (which requires that both T* and its derivative are matched at x_s). The fifth condition is the ice-edge condition ${T}^{* }({x}_{s})=1$, from which we derive a relationship between unknown model parameters valid for a given x_s (see Appendix B):

$$\begin{eqnarray}{q}^{-1}=\left\{\begin{array}{cc}1+\tfrac{{s}_{20}}{1+6\delta }{P}_{2}({x}_{s})-\alpha F({x}_{s}), & {s}_{20}\lt 0\\ (1-\alpha )\left(1+\tfrac{{s}_{20}}{1+6\delta }{P}_{2}({x}_{s})\right)+\alpha F({x}_{s}), & {s}_{20}\gt 0\end{array}\right.,\end{eqnarray} \tag{ 25 }$$

where $F({x}_{s})$ is a special function computed in terms of hypergeometric functions, and depends on model parameters δ and ${s}_{20}/(1+6\delta );$ the formula is given in Equation (43). Equation (25) gives the radiative forcing q required for a given ice edge x_s. Although it seems perverse to solve for $q({x}_{s})$ rather than x_s(q), the inverse of Equation (25) is multiple valued and does not have a closed form.

The expressions in Equation (25) define the solution space for the EBM. If the planetary properties $\delta ,\alpha$, and β (and hence s₂₀) are known, then from Equation (25) we can calculate the radiative forcing q necessary for any ice edge. We might also think of q as known and use Equation (25) to compute the albedo or transport efficiency necessary for a given ice edge. The point is that Equation (25) defines a relationship between the four model parameters for any ice edge x_s.

Note that s₂₀ only appears in the solutions (including in the definition of F) as the ratio ${s}_{20}/(1+6\delta )$. This factor arises because the particular solutions to any forcing term ${P}_{2}(x)$ in the diffusion Equations (14) are damped by a factor $1/(1+6\delta );$ diffusion acts to smear out the local forcing, producing a weaker temperature gradient.

In Figure 5, we plot the solutions $q({x}_{s})$ for a range of values of δ and α, and for both Earth-like and high obliquity. Because the location of the icy region is reversed at high obliquity (Figure 1), we have reversed the y-axis in the $\beta =90^\circ$ plots in the lower row so that up is always cold and toward an ice edge (if it exists). This simplifies the visual interpretation of the stability criterion discussed below. These solutions are multiple valued: there are anywhere from one to five different equilibrium ice edges for any given radiative forcing (which we can visualize by drawing a vertical line through the plot). At most three of these solutions are stable equilibria (and thus physically realizable).

Zoom In Zoom Out Reset image size

The physical basis of the stability condition is straightforward: the equilibrium size of the icy region should decrease with an increase in radiative forcing, so that small perturbations in the ice extent will excite a net negative radiative feedback and decay. For polar ice caps at low obliquity, the ice edge is stable if

$$\begin{eqnarray}&&\displaystyle \frac{{dq}}{{{dx}}_{s}}\gt 0,\end{eqnarray} \tag{ 26a }$$

a condition known as the "slope-stability theorem" (Cahalan & North 1979). Locations on the graphs where the slope ${dq}/{{dx}}_{s}$ changes sign therefore indicate bifurcation points of the ice−climate system. For the familiar low-obliquity case, there are both minimum and maximum sizes for stable polar ice caps (North 1984; Roe & Baker 2010). The associated instabilities have been called the Small Ice-Cap Instability (SICI) and Large Ice-Cap Instability (LICI) or simply Snowball instability.

The green curve in Figure 5(b) is reasonably close to Earth-like parameters (β = 23.5°, δ = 0.32, α = 0.44, q = 1.2) and serves to illustrate the classic Snowball Earth hysteresis. The analogue of the present-day climate sits on the stable ice-cap branch with x_s near 60° latitude in this case (though x_s is closer to 70° for Earth). A reduction in greenhouse gases or solar constant sufficient to reduce q by about 2% would cause the climate to cool and the ice cap to expand, following the green curve down to the bifurcation point around 38° latitude. A further reduction in q would then trigger unstable ice expansion down to the equator (the horizontal line at the bottom of the graph). Exiting the Snowball scenario would then require sufficient radiative forcing to raise the equatorial surface temperature past the melting point, i.e., $q\geqslant {q}_{\mathrm{snow}}$. From Equation (22), this threshold is about 1.65, or about a 38% increase over the present-day value of 1.2. The planet would then undergo an unstable transition to extremely warm and ice-free conditions (the horizontal line at the top of the graph). Qualitatively, the same statements can be made for any curve in the upper panels that features a stable ice cap over some range of q.

Analogous stability diagrams can be easily computed for any obliquity value. In particular, the perpetual equinox case $\beta =0$ is qualitatively similar to the upper row of Figure 5 but somewhat more stable (not shown). However, this result should be treated with caution because of the misfit of the series approximation to $\overline{s(x)}$ near the poles for $\beta =0$ (Figures 2 and 3).

For high obliquity, since the ice sits equatorward of x_s, the stability condition is reversed. A stable equatorial ice belt requires

$$\begin{eqnarray}&&\displaystyle \frac{{dq}}{{{dx}}_{s}}\lt 0.\end{eqnarray} \tag{ 26b }$$

Since we have reversed the axes in the high-obliquity plots in Figures 5(e)–(f), the stable solution branches are still readily identifiable in all cases as regions of positive slope. As in the low-obliquity case, there are minimum and maximum sizes for stable ice belts (latitudes at which the slope of the graphs in Figures 5(e)–(f) changes sign). Here we introduce new terminology to describe these bifurcations: the Small Ice-Belt Instability (SIBI) and Large Ice-Belt Instability (LIBI) near the top and bottom of the graphs, respectively.

Stable finite ice cover is far from universal throughout the parameter space at both high and low obliquity. There are, for example, no stable ice edges for very large δ (blue curves); more modest δ admits stable solutions only for small α (weak albedo feedback, e.g., magenta and green curves). Stable ice edges are possible when neither δ nor α is large. In fact, we have already established in Equations (23) and (24) that the product $\alpha (1+6\delta )$ governs the existence of multiple equilibria due to albedo feedback, and the existence of stable ice edges appears to require that this same product not exceed some threshold. This can also be illustrated through the following approximate stability condition, following Rose (2015).

Suppose that an equilibrated climate is subject to a small globally uniform warming or cooling ${\rm{\Delta }}{T}^{* }$, which displaces the ice area by ${\rm{\Delta }}{x}_{s}$. This perturbation is stable if and only if the negative temperature feedback outweighs the positive albedo feedback:

$$\begin{eqnarray}&&{\rm{\Delta }}{T}^{* }\gt \alpha {qs}({x}_{s})| {\rm{\Delta }}{x}_{s}| ,\end{eqnarray} \tag{ 27 }$$

(where the absolute value accounts for the reversed ice orientation at high obliquity). For small perturbations, the ice-edge displacement is determined by the local temperature gradient: ${\left.{\rm{\Delta }}{T}^{* }=-{\rm{\Delta }}{x}_{s}\tfrac{{{dT}}^{* }}{{dx}}\right|}_{{x}_{s}}$. Taking temperature gradients from the linear solutions (20) then gives

$$\begin{eqnarray}&&\alpha (1+6\delta )\lt \displaystyle \frac{3{x}_{s}| {s}_{20}| }{1+{s}_{20}{P}_{2}({x}_{s})}.\end{eqnarray} \tag{ 28 }$$

Evaluating this condition near x_s = 0.5 (half global ice cover) gives, to a good approximation,

$$\begin{eqnarray}&&\alpha (1+6\delta )\lt \displaystyle \frac{3}{2}| {s}_{20}| \end{eqnarray} \tag{ 29 }$$

(effectively the opposite of Equation (24)). This is clearly an overly stringent condition because the linear solutions (20) underestimate the temperature gradient near the ice edge, but it serves to illustrate the dependence of stable branches on α and δ. Furthermore, Equation (29) suggests that the condition is more stringent for high obliquity where insolation gradients are weaker, as found in Figure 5.

The exact stability condition can be derived from Equation (25) by setting ${dq}/{{dx}}_{s}=0$ and solving for α. Stable ice edges are possible wherever $\alpha \lt {\alpha }_{\mathrm{crit}}({x}_{s},\delta ,{s}_{20})$, where the critical value takes on different forms for low and high obliquity:

$$\begin{eqnarray}{\alpha }_{\mathrm{crit}}=\left\{\begin{array}{cc}\tfrac{{s}_{20}{P}_{2}^{{\prime} }({x}_{s})}{(1+6\delta ){F}^{{\prime} }({x}_{s})}, & {s}_{20}\lt 0\\ {\left(1-\tfrac{(1+6\delta ){F}^{{\prime} }({x}_{s})}{{s}_{20}{P}_{2}^{{\prime} }({x}_{s})}\right)}^{-1}, & {s}_{20}\gt 0\end{array}\right.,\end{eqnarray} \tag{ 30 }$$

where the primes refer to derivatives with respect to x evaluated at x_s. ${\alpha }_{\mathrm{crit}}$ is shown in Figures 5(d) and (h). Again, as either α or δ increases, the range of latitudes at which stable ice caps or ice belts can be found is reduced.

The model predicts that stable ice edges should be rarer under high obliquity than low obliquity. For example, for Earth-like $\alpha =0.44$ but 90° obliquity, stable ice edges occur only for weak heat transport ($\delta =0.16$, red curve, and smaller), and these are very sensitive (large ice-edge displacements for small changes in radiative forcing) relative to their corresponding solutions at Earth-like obliquity. For large parts of the parameter space, the only physically realizable solutions are the ice-free and Snowball states, which are plotted as horizontal lines at 0° and 90° in Figure 5.

Finally, Figure 5 also shows that stable ice edges are frequently found in regions of the parameter space in which both ice-free and Snowball solutions are also possible. This is particularly true for the stable ice belts at high obliquity. For example, the red curve in Figure 5(f) has a stable ice-belt branch, but the entire branch coexists with the ice-free and Snowball branches. There are thus no possible transitions from either ice-free or Snowball conditions into this stable ice-belt state involving a hysteresis in q (radiative forcing). In this case, the ice belt is possible but difficult to realize as there is no ready mechanism by which the climate can drift into this state from warmer or colder initial conditions. The abundance of inaccessible ice belts may help explain why ice-belt states have only rarely been found in GCM studies (see our introduction).

As we argued above, the annual mean EBM predicts that stable high-obliquity ice belts should be rarer than stable low-obliquity ice caps. Here we quantify this statement under some simple, speculative assumptions about the probability distributions of planetary parameters. We will evaluate the relative likelihood L_ice of finding a stable ice edge as a function of obliquity, where "relative" here means that L_ice will be normalized by the probability for the same planetary properties but Earth's current obliquity.

Suppose that we have a known probability distribution of planetary properties

$$\begin{eqnarray}&&{h}_{\mathrm{planet}}(q,\delta ,\alpha ).\end{eqnarray} \tag{ 31 }$$

In principle, given δ, α, and q we can calculate the ice edge x_s and its stability. Since we have a functional form for $q({x}_{s})$ rather than x_s(q), we compute probability of a stable ice edge by integrating over all possible ice edges:

$$\begin{eqnarray}&&{P}_{\mathrm{ice}}(\beta )\\ &&\ =\,\displaystyle \frac{{\displaystyle \int }_{0}^{1}{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{{\alpha }_{\mathrm{crit}}}{h}_{\mathrm{planet}}(q({x}_{s},\delta ,{s}_{20},\alpha ),\delta ,\alpha )\,d\alpha \,d\delta \,{{dx}}_{s}}{{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{1}{h}_{\mathrm{planet}}(q,\delta ,\alpha )\,d\alpha \,d\delta \,{dq}},\end{eqnarray} \tag{ 32 }$$

where the denominator is equal to unity. By taking the integral in the numerator only over the interval $0\lt \alpha \lt {\alpha }_{\mathrm{crit}}({x}_{s},\delta ,{s}_{20})$, we sample just the part of the solution space that contains stable ice edges. We then normalize by the value of P_ice at Earth's obliquity:

$$\begin{eqnarray}&&{L}_{\mathrm{ice}}(\beta )=\displaystyle \frac{{P}_{\mathrm{ice}}(\beta )}{{P}_{\mathrm{ice}}({\beta }_{\oplus })}.\end{eqnarray} \tag{ 33 }$$

There is one further complication involving the inaccessible stable solution branches discussed at the end of Section 5.1. We exclude from our probability calculation all states that are inaccessible from either the ice-free or Snowball branch through a hysteresis in q. The method is illustrated in Figure 6. For given values of δ and s₂₀, ${\alpha }_{\mathrm{warm}}$ (magenta curve) is defined implicitly by the solution of $q({x}_{s},\alpha ={\alpha }_{\mathrm{crit}})={q}_{\mathrm{free}}$. Similarly, ${\alpha }_{\mathrm{cold}}$ (cyan curve) is defined implicitly by the solution of $q({x}_{s},\alpha ={\alpha }_{\mathrm{crit}})={q}_{\mathrm{snow}}$ evaluated at ${\alpha }_{\mathrm{crit}}$. These conditions describe the parameter space boundaries for transitions to the stable ice-edge branch, respectively, from the ice-free and Snowball branches. We then reduce the limit of the integration in Equation (32) from ${\alpha }_{\mathrm{crit}}$ to a smaller value ${\alpha }_{\max }$ given by

$$\begin{eqnarray}&&{\alpha }_{\max }({x}_{s},\delta ,{s}_{20})=\min ({\alpha }_{\mathrm{crit}},\max ({\alpha }_{\mathrm{warm}},{\alpha }_{\mathrm{cold}})).\end{eqnarray} \tag{ 34 }$$

The maximum value in Equation (34) ensures that we include all stable ice edges that are accessible (through a radiative hysteresis) from at least one side.

Zoom In Zoom Out Reset image size

We are unable to solve analytically for ${\alpha }_{\max }$ (it involves a transcendental equation for x_s), but a numerical solution with a root-finding algorithm is straightforward. We thus compute two version of the likelihood L_ice: one in which all possible stable ice edges are accounted for using Equation (32), and another in which inaccessible stable states are excluded by replacing ${\alpha }_{\mathrm{crit}}$ with ${\alpha }_{\max }$ in the limits of integration.

Calculating L_ice requires a plausible form for the probability distribution h_planet. We assume here that the three parameters $q,\delta ,\alpha$ are independent of each other, so that their joint probability distribution is separable:

$$\begin{eqnarray}&&{h}_{\mathrm{planet}}(q,\delta ,\alpha )={h}_{q}(q)\,{h}_{\delta }(\delta )\,{h}_{\alpha }(\alpha ).\end{eqnarray} \tag{ 35 }$$

Alternatively, we might consider some parameter interdependence, e.g., α and δ might co-vary due to the effects of stronger albedo contrasts on atmospheric baroclinicity and eddy activity. We make the separable assumption here for simplicity. Recall that $q={a}_{0}Q{(A+{{BT}}_{\mathrm{ref}})}^{-1}$ and $\delta ={{KR}}^{-2}{B}^{-1}$ are products of several parameters—each of which may vary considerably—and both variables are also positive definite, with no upper bound. These features suggest that taking log-normal distributions for h_q and h_δ may be reasonable. We lack such a firm basis for the form of h_α, but as α is bounded on $[0,1]$, both uniform and beta distributions seem reasonable.

We adopt three different sets of assumptions, resulting in three different curves in Figure 7. For PDF0 (blue curve), we assume the following: h_α is uniform on [0, 1]; h_δ is log-normal on $[0,\infty ]$ with shape parameter 1.0, scale parameter 1.0, and location parameter 0 (mode at $\delta =0.37$, median at $\delta =1$); and h_q is log-normal on $[0,\infty ]$ with shape parameter 0.5, scale parameter 1.0, and location parameter 0 (mode at q = 0.78, median at q = 1). PDF1 (green curve) is the same as PDF0 except h_δ gives more weight to both small and large values of δ (log-normal with shape parameter 2.0 and scale parameter e; mode at $\delta ={e}^{-3}$, median at $\delta =e$). Finally, PDF2 (red curve) is the same as PDF1, except h_α is a parabolic beta distribution on $[0,1]$ with a mode at $\alpha =0.5$, thus giving a higher probability for moderating albedo feedback. In each case, L_ice is computed through numerical integration of Equation (32) for a given obliquity value and normalized by its value at Earth obliquity.

Zoom In Zoom Out Reset image size

Figure 7 shows that in all cases the likelihood of a stable ice cap decreases with increasing obliquity, approaching zero at the critical obliquity of 55°. The likelihood of a stable ice belt then increases with obliquity between 55° and 90°, but in all cases is lower than the likelihood of a stable ice cap at Earth-like obliquity. The details appear to be most sensitive to assumptions about ${h}_{\alpha }$ (i.e., the red curve is substantially different from the green and blue curves). However, in all cases the shape of L_ice is qualitatively similar to the annual mean equator-to-pole insolation gradient $| {s}_{20}|$. As this gradient becomes small near $\beta =55^\circ$, the transport will have an increasingly strong tendency to overwhelm the insolation gradient and produce a nearly isothermal climate either above or below the ice threshold.

At high obliquity, the exclusion of inaccessible stable states has a larger impact on L_ice than at low obliquity (as anticipated by the shapes of the $q({x}_{s})$ curves in Figure 5). The relative rareness of radiative hysteresis loops that access the stable ice-belt configuration means that such states are less likely to be observed (or modeled). The tentative conclusion (based on these highly speculative assumptions about the distributions of planetary parameters) is that ice belts will be less probable than ice caps, by a factor of 2 or so.

We can go one step further toward a general statement about the observability of ice belts versus ice caps by making an assumption about the distribution of obliquities. If, for example, obliquity is uniformly distributed between 0° and 90°, the area under the curves in Figure 7 on either side of the critical 55° define relative likelihoods of finding stable ice caps versus ice belts. Under this assumption, the proportion of all hypothetically observable stable ice edges comprised of ice belts is then 22%, 24%, or 17% using PDF0, PDF,1 and PDF2, respectively. So, about three-fourths to four-fifths of all planets with partial ice cover would be in the form of Earth-like low-obliquity polar ice caps. Note that these numbers just represent relative likelihoods of ice caps versus ice belts; they do not represent probabilities of finding any of these stable ice-edge states relative to ice-free or Snowball conditions.

Is uniformly distributed obliquity a reasonable prior? Miguel & Brunini (2010) simulated planetary accretion and found a primordial obliquity distribution peaking at 90°. Taking ${h}_{\beta }=\sin \beta$ to fit these results better gives 46%, 48%, or 39% ice belts using PDF0, PDF1, and PDF2, respectively. However, it is not clear that observable obliquities have the same distribution as primordial obliquities, as primordial obliquities can be altered by a variety of mechanisms (e.g., Laskar & Robutel 1993; Correia & Laskar 2001). In particular, high obliquities will be reduced by strong tidal interactions for planets close to their host star, yielding a distribution closer to uniform. Better constraints on obliquity distributions and all planetary parameters are clearly warranted.

A question of relevance to life on other planets is whether the Snowball catastrophe occurs at lower radiative forcing at low or high obliquity. There is some conflict in the literature on the role of the Snowball bifurcation in planetary habitability. On the one hand, defining habitability in terms of surface liquid water suggests that a Snowball climate is uninhabitable, which has led various authors to propose metrics of fractional or seasonal habitability for planets with partial ice cover (e.g., Williams & Pollard 2003; Spiegel et al. 2008). On the other hand, not only did photosynthesis persist through Snowball events in Earth' history, but the events may have crucially shaped the subsequent evolution of complex life (e.g., Hoffman & Schrag 2002; Hoffman et al. 2017; Laasko & Schrag 2017). The traditional habitable zone concept assumes a planet with a functioning silicate-weathering feedback and a positive CO₂ greenhouse effect (Kasting et al. 1993; Kopparapu et al. 2013). Global glaciation may be triggered on such a planet through a rapid drawdown of atmospheric CO₂ that reduces q below the thresholds at which the non-Snowball states disappear. In Earth history, this seems to have occurred through accidents of tectonics (Hoffman et al. 2017). These events are self-terminating through the suppression of silicate weathering and accumulation of a strong CO₂ greenhouse. Such transitions can in principle occur anywhere within the habitable zone. Our simple model is unsuited to the tasks of diagnosing the inner edge of habitability (where the relevant physics are the runaway water vapor greenhouse and hydrogen loss to space) or the outer edge (where the relevant physics are CO₂ condensation and Rayleigh scattering). However, we can compute the q value at which the bifurcation occurs as a function of obliquity and other model parameters, and whether the transition into the Snowball state occurs from a partially ice-covered state (cap or belt) or directly from the ice-free state.

We therefore define q_hab as the smaller of the minimum q necessary for ice-free conditions (given by Equation (21)) and the minimum q for which a stable ice edge exists (if any; a quantity we will call q_stab), i.e.,

$$\begin{eqnarray}&&{q}_{\mathrm{hab}}=\min ({q}_{\mathrm{stab}},{q}_{\mathrm{free}}).\end{eqnarray} \tag{ 36 }$$

As discussed above, for large values of δ and α, there are no stable ice edges, so q_stab is undefined. Elsewhere it is given implicitly by the solution of $\alpha ={\alpha }_{\mathrm{crit}}({x}_{s})$ for the critical ice-edge latitude x_crit (the bifurcation point at which ${dq}/{{dx}}_{s}=0$ and the stable branch ends in Figure 5). We evaluate this numerically for given values of ${s}_{20},\delta ,\alpha$.

Figure 8 shows contour plots of q_hab as a function of obliquity and heat transport efficiency δ, for the same three values of α used in Figure 5. These are plotted so that darker colors (red to black) indicate planets that remain Snowball-free under weaker radiative forcing.

Zoom In Zoom Out Reset image size

How does the ice−albedo feedback influence the location and nature of the Snowball bifurcation? Figure 8 illustrates a number of different issues. The Snowball catastrophe is moistly avoided on planets with very weak albedo feedback (small α). Such planets can maintain very narrow bands of open water even when the radiative forcing is weak, and this is nearly equally true for low- and high-obliquity planets. On the other hand, for more than a trivial albedo change across the ice edge, some clear differences emerge between low- and high-obliquity worlds. The differences are best illustrated by the black contours, which show the value of δ above which ${q}_{\mathrm{hab}}={q}_{\mathrm{free}}$. Below this line, the Snowball bifurcation is reached from partial ice cover (ice caps or ice belt). Above this line, the bifurcation is reached from ice-free conditions. Figure 8 shows clearly that this threshold occurs at substantially smaller δ values for high obliquity.

For a concrete example, consider the middle panel of Figure 8 with $\alpha =0.44$. At the Earth-like value of $\delta =0.31$, we have ${q}_{\mathrm{hab}}=1.18$ for $\beta =23.45^\circ$, but ${q}_{\mathrm{hab}}={q}_{\mathrm{free}}=1.05$ for $\beta =90^\circ$. The high-obliquity planet remains Snowball-free for substantially weaker radiative forcing, and remains entirely ice-free before transitioning to the Snowball for all but very small δ. With any significant albedo feedback, a given planet is less likely to fall into a Snowball state at mid-to-high obliquity than at low obliquity.

In summary, there are essentially two rather different scenarios that favor Snowball avoidance at low radiative forcing: either weak albedo feedback and inefficient heat transport, or efficient heat transport at mid- to high obliquity. The first scenario might be associated with large, rapidly rotating planets orbiting M-dwarfs, where a partially unfrozen region is maintained by relatively weak albedo feedback and inefficient heat transport due to strong zonal flow. The second scenario might be associated with smaller, slowly rotating planets at high obliquity, which would remain completely ice-free even with a potentially strong albedo feedback (large α) from shorter wavelength insolation. Albedo feedback is more likely to drive large climate changes under low obliquity than high obliquity.

We emphasize that these calculations are relevant for planets with oceans anywhere within the habitable zone. Near the outer edge of the habitable zone, a planet with a functioning silicate-weathering feedback should have multiple bars of atmospheric CO₂. The resultant strong atmospheric scattering would effectively mask the surface albedo, such that the relevant α value would be small (Wordsworth et al. 2011).