Comparison of Two Analytic Energy Balance Models Shows Stable Partial Ice Cover Possible for Any Obliquity

The search for habitable exoplanets, perhaps hosting life, is one of the great endeavors of our time. To aid the search, it is important to understand what planetary factors contribute to a planet's habitability for life as we know it. In this direction, analytical and computational studies using climate models to investigate different planetary scenarios have been incredibly important for advancing our scientific understanding of exoplanet climates. Recent work has investigated how a planet's orbital parameters, such as obliquity (e.g., Armstrong et al. 2014; Rose et al. 2017), eccentricity (e.g., Kane et al. 2020), or spin–orbit resonance (e.g., tidal locking; Checlair et al. 2017, 2019a, 2019b), may affect the planet's climate. Other studies have looked at what role climatic elements such as sea ice drift (e.g., Yue & Yang 2020), land albedo (e.g., Rushby et al. 2020), or ocean heat transport (e.g., Checlair et al. 2019a) may play in the long-term habitability of a planet. This study adds to the rapidly growing body of literature on exoplanet climate by considering the role of albedo/temperature feedback in a zonally averaged energy balance model with explicit obliquity dependence. We compare the model output for two methods of modeling the energy transport across latitudes and find qualitatively similar results between them.

This study is an extension of the work in Rose et al. (2017). In that study, Rose et al. (2017) analyzed an energy balance model with explicit obliquity dependence and heat transport modeled with a diffusion term. They found an analytical solution for temperature as a function of the sine of the latitude and obliquity based on the methods given in North (1975a). Their work shows that partial stable ice cover is more likely for ice caps than for ice belts (the model scenario where ice advances from the equatorial region of the planet rather than the poles; e.g., as in Figure 1) and identifies regions of parameter space where the snowball catastrophe—where a small change in radiative forcing causes the planet to quickly become completely ice-covered—may be avoided.

Zoom In Zoom Out Reset image size

Here we compare the model used in Rose et al. (2017) in two ways: (1) to the same model but with a more accurate approximation to the insolation distribution and (2) to a similar energy balance model but with heat transport modeled with relaxation to the global mean average temperature. Comparisons between energy balance models with diffusive and relaxation terms for heat transport have been conducted in the past for Earth (e.g., North 1984; Roe & Baker 2010; Walsh 2017); however, an analysis agnostic to the specificity of particular planets has yet to appear in the literature.

In Rose et al. (2017), the incoming stellar radiation is approximated by a second-degree polynomial in the sine of the latitude and cosine of the obliquity. A second-degree polynomial is sufficient to capture the qualitative behavior of insolation distributions for planets with very low obliquities (between 0° and 45°) and obliquities close to 90° (between 65° and 90°) but does not capture the qualitative distribution for planets between these ranges. ³ To capture the behavior accurately, one needs at least a sixth-degree polynomial in the sine of the latitude and cosine of the obliquity (Nadeau & McGehee 2017; Dobrovolskis 2021). Here we show that the higher-order approximation is a necessary one; we would not find that stable partial ice cover is possible at any obliquity without it. In particular, we find that stable partial ice cover is possible for any obliquity for both diffusive and relaxation to the mean models.

Further, we find that the mode of heat transport does not change the qualitative distribution of the likelihood of planets with stable partial ice cover. Here we show that stable partial ice cover is less likely for high-obliquity planets than for low-obliquity planets for both diffusive and relaxation to the mean models, agreeing with the results already shown in Rose et al. (2017). We also show that low albedo contrast and low efficiency of heat transport favor stable partial ice cover. The qualitative similarities between the models are encouraging because although analytical solutions can be found in both cases, the analytical solutions for the relaxation to the mean model are more tractable, as they can be written using elementary functions.

Following Rose et al. (2017), we consider the stability of the equilibrium of the system as radiative forcing changes in the model. Physically, this could be caused by atmospheric effects, such as changing greenhouse gases. In particular, we consider how different parameters in the model affect the snowball catastrophe. We consider the severity of the snowball catastrophe bifurcation based on the distance between the bifurcation latitude and the snowball state. We identify the regions in parameter space where more severe bifurcations occur, since a severe snowball catastrophe may have stronger signals in a planet's climate record and may be more likely to be observed. We show that the snowball catastrophe is less severe in the relaxation to the mean model compared to the diffusion model.

Our paper is laid out as follows. In Section 2.1, we present the governing equations and a nondimensionalization of these equations to better quantify the effects of parameter changes on the behavior of the system. In Section 2.2, we derive the equations that relate the latitude of the saddle node bifurcation to the corresponding parameter values. In Section 3, we calculate the relative likelihood of stable partial ice cover. In Section 4, we classify the bifurcation into the snowball state based on its severity. A discussion of the results follows in Section 5, and we conclude in Section 6. In the Appendix, we briefly analytically compare the equilibrium solutions of the diffusion and relaxation to the mean models.

2.1. Annual Average Energy Balance

We consider a one-dimensional energy balance model of the form popularized by Budyko (1969), Sellers (1969), and North (1975b). These models describe the time evolution of the temperature on a planet depending on the incoming and outgoing radiation and heat transfer across latitudes. We consider the model

$$\begin{eqnarray}&&R\displaystyle \frac{\partial T}{\partial t}={Qs}(y,\beta )(1-\alpha (y,\eta ))-(A+{BT})+{ \mathcal F }(T),\end{eqnarray} \tag{ 1 }$$

where T = T(y, t) is the annual zonally averaged temperature as a function of the sine of latitude y and time t, and ${ \mathcal F }(T)$ will either be diffusion,

$$\begin{eqnarray*}&&{ \mathcal F }(T)=D{{\rm{\nabla }}}^{2}T=D\displaystyle \frac{\partial }{\partial y}(1-{y}^{2})\displaystyle \frac{\partial }{\partial y}T(y,t),\end{eqnarray*}$$

or relaxation to the global mean temperature,

$$\begin{eqnarray*}&&{ \mathcal F }(T)=-C(T-\overline{T})=-C\left(T(y,t)-{\int }_{0}^{1}T(\phi ,t)d\phi \right).\end{eqnarray*}$$

As noted in Roe & Baker (2010), both forms of ${ \mathcal F }(T)$ are parameterizations of the divergence of the poleward heat flux. Scientific discussions of the terms in the relaxation to the mean model can be found in the literature, e.g., Held & Suarez (1974) and Checlair et al. (2017). Readers interested in a mathematical discussion should see, for example, North (1975b), Tung (2007), Kaper & Engler (2013), Widiasih (2013), McGehee & Widiasih (2014), and Walsh (2017). In the above model, y is the sine of the latitude. Hemispheric symmetry is assumed here so that the latitude y ranges from zero (the equator) to 1 (the north pole). The ice line latitude (the boundary between the high- and low-albedo regions) is denoted by η. The mean annual amount of incoming solar radiation (insolation) is represented by Q. The insolation distribution function s(y, β) depends on the latitude and the planetary axial tilt β. The co-albedo function (1 − α(y, η)) determines the proportion of incoming solar radiation absorbed by the planetary surface at each latitude. Outgoing radiation A + BT is in a linearized form. The dimensional parameters and their values for the Earth can be found in many references, such as Tung (2007) or Kaper & Engler (2013). Nondimensionalized parameter values for Earth are given in Table 1.

Table 1. Nondimensional Parameters for Equations (4) and (5) ^a

Parameter	Definition	Brief Description	Value for Earth
γ	$\tfrac{R\omega }{B}$	Seasonal heat capacity	6.13
q	$\tfrac{(1-{\alpha }_{l})Q}{A+{{BT}}_{c}}$	Radiative forcing	1.27
$\overline{\alpha }$	$1-\tfrac{1-{\alpha }_{h}}{1-{\alpha }_{l}}$	Albedo contrast	0.44
ζ	$\cos (\beta )$	Cosine of obliquity	0.92
λ	$\tfrac{\rho (A+{{BT}}_{c})}{B\omega }$	Ice line response	Varies
N		2N is the degree of the insolation approximation	1
μ	$\tfrac{C}{B}$	Relaxation to mean efficiency of heat transport	1.6
δ	$\tfrac{D}{B}$	Diffusion efficiency of heat transport	0.31

Note.

^a The dimensional parameter values are the same as in Tung (2007), except D, which is chosen to be consistent with Rose et al. (2017) and tuned to reflect the conditions of the modern-day climate (e.g., North 1975b).

Download table as:  ASCIITypeset image

The albedo function α(y, η) is a piecewise constant function demarcating the regions of the planet with high and low albedo. Let α_p be the albedo poleward of the ice line and α_e be the albedo equatorward of the ice line; then,

$$\begin{eqnarray*}\alpha (y,\eta )=\left\{\begin{array}{ll}{\alpha }_{e} & 0\lt y\lt \eta \\ \displaystyle \frac{{\alpha }_{e}+{\alpha }_{p}}{2} & y=\eta \\ {\alpha }_{p} & \eta \lt y\lt 1.\end{array}\right.\end{eqnarray*}$$

In the case of ice caps, the polar region has a higher albedo than the equatorial region with α_p > α_e . For ice belts, the situation is reversed, and α_p < α_e (Figure 1). In the following, we will denote high albedo with α_h (ice-covered regions) and low albedo with α_l (non-ice-covered regions). In some energy balance models of Earth, the albedo values are set to α_l = 0.32 and α_h = 0.62 with α_e = α_l and α_p = α_h (e.g., in Widiasih 2013). The ice line–dependent albedo is used in McGehee & Widiasih (2014), Widiasih (2013), Walsh (2017), and Barry et al. (2017). This is in contrast with the temperature-dependent albedo used in Rose et al. (2017) and North (1975b). Here we use the ice-dependent version because, while the two forms are equivalent in the diffusion model (Cahalan & North 1979), the temperature-dependent version in the relaxation to the mean model can result in spurious temperature solutions (Walsh & Rackauckas 2015).

We model the annual average changes in the temperature profile by using zonally averaged mean annual insolation in Equation (1). The mean annual insolation distribution s(y, β) depends only on obliquity β and latitude y. The insolation distribution is given by the integral (e.g., Ward 1974; McGehee & Widiasih 2014)

$$\begin{eqnarray*}&&\begin{array}{l}s(y,\beta )\\ \,=\displaystyle \frac{2}{{\pi }^{2}}{\int }_{0}^{2\pi }\sqrt{1-{\left(\sqrt{1-{y}^{2}}\sin \beta \cos \gamma -y\cos \beta \right)}^{2}}d\gamma ,\end{array}\end{eqnarray*}$$

where γ is the planet's longitude and which is not amenable to analytical calculations in the models we consider. Instead, Nadeau & McGehee (2021) showed that we may write s(y, β) as a Legendre series,

$$\begin{eqnarray*}&&s(y,\beta )=\displaystyle \sum _{n=0}^{\infty }{a}_{2n}{p}_{2n}(\cos \beta ){p}_{2n}(y),\end{eqnarray*}$$

where p_2n is the Legendre polynomial of degree 2n and

$$\begin{eqnarray*}&&{a}_{2n}=\displaystyle \frac{{\left(-1\right)}^{n}(4n+1)}{{2}^{2n-1}}\sum _{k=0}^{n}\displaystyle \left(\genfrac{}{}{0em}{}{2n}{n-k}\right)\displaystyle \left(\genfrac{}{}{0em}{}{2n+2k}{2k}\right)\displaystyle \left(\genfrac{}{}{0em}{}{1/2}{k\,+\,1}\right)\end{eqnarray*}$$

using the standard notation

$$\begin{eqnarray*}&&\displaystyle \left(\genfrac{}{}{0em}{}{x}{k}\right)=\displaystyle \frac{x(x-1)\cdots (x-k+1)}{k!}.\end{eqnarray*}$$

Truncating the series for some fixed N gives a degree 2N polynomial approximation σ_2N to the true insolation distribution, where increasing N decreases the rms error of the approximation (Nadeau & McGehee 2021). Nadeau & McGehee (2017, 2021) showed that the sixth-degree approximation given by

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{6}(y,\beta ) & = & \displaystyle \sum _{n=0}^{3}{a}_{2n}{p}_{2n}(\cos \beta ){p}_{2n}(y)\\ & = & 1+{a}_{2}{p}_{2}(\cos \beta ){p}_{2}(y)+{a}_{4}{p}_{4}(\cos \beta ){p}_{4}(y)\\ & & +{a}_{6}{p}_{6}(\cos \beta ){p}_{6}(y),\end{array}\end{eqnarray} \tag{ 2 }$$

where a₂ = −5/8, a₄ = −9/64, and a₆ = −65/1024, and

$$\begin{eqnarray*}\begin{array}{rcl}{p}_{2}(y) & = & (3{y}^{2}-1)/2,\\ {p}_{4}(y) & = & (35{y}^{4}-30{y}^{2}+3)/8,\\ {p}_{6}(y) & = & (231{y}^{6}-315{y}^{4}+105{y}^{2}-5)/16\end{array}\end{eqnarray*}$$

is the smallest degree approximation that captures the qualitative shape of the distribution for all obliquities and approximates the true distribution to within 1.6% for all obliquities. In particular, it is the smallest degree approximation that captures the characteristic "W" shape of the insolation distributions for planets with obliquity between approximately 45° and 65° (Nadeau & McGehee 2017). Dobrovolskis (2021) showed that extending the approximation to degree 8 or 10 improves the approximation significantly at the poles for obliquities close to zero but only slightly at other latitudes or for higher obliquities. In Rose et al. (2017), the degree 2 approximation is used (a₄ = 0 and a₆ = 0). Notice that when $\cos \beta =\sqrt{3}/3$ (i.e., when β ≈ 54.74°), the degree 2 approximation is identically equal to 1 for all values of y. Here we use a generic truncation σ_2N in our analysis. In Section 3, we present results for the cases where N = 1 and 3.

The position of the ice lines, denoted by η, depends on the mean annual temperature of the ice line. Ice–albedo feedback is incorporated by the dynamic ice line equation that is coupled with Equation (1). We use the following dynamic ice line equation first formulated in Widiasih (2013) for ice caps on Earth as

$$\begin{eqnarray}&&\displaystyle \frac{d\eta }{{dt}}=\rho (T(\eta ,t)-{T}_{c}).\end{eqnarray} \tag{ 3 }$$

For ice belts, the right-hand side should be multiplied by –1. The physical boundaries at the pole and equator are built into the model; i.e., η cannot be greater than 1 or less than zero. A mathematical treatment of the resulting nonsmooth system using a projection rule and a Filippov framework can be found in Barry et al. (2017), where the invariance of the physically possible region is shown.

An intuitive interpretation of Equation (3) notes that the mean annual temperature at the ice line is denoted by T(η, t). The critical temperature T_c is the highest temperature at which multiyear ice can be present. If the ice line temperature is above T_c , then the ice cover shrinks. If the temperature at the ice line is below T_c , the ice cover grows. The response constant ρ controls the speed of the ice line response to a change in temperature. We are interested in the equilibrium position of the ice line, which is obtained when the ice line temperature is exactly T_c .

We nondimensionalize the system using transformations analogous to those in Rose et al. (2017), namely,

$$\begin{eqnarray*}&&\tau =\omega t=\displaystyle \frac{2\pi t}{{t}_{\mathrm{year}}},{T}^{* }=\displaystyle \frac{A+{BT}(y)}{A+{{BT}}_{c}},\end{eqnarray*}$$

where ω = 2π/t_year rescales time to be dimensionless, such that τ = 1 corresponds to 1 yr. The nondimensionalized temperature T* is proportional to temperature and outgoing longwave radiation. At the ice line at equilibrium given by Equation (6), T* is always equal to 1, i.e., T*(η) = 1.

The nondimensionalized parameters are summarized in Table 1. The parameter transformations (which are the same for both diffusion and relaxation to the mean) are

$$\begin{eqnarray*}\begin{array}{lll}\gamma =\displaystyle \frac{R\omega }{B}, & q=\displaystyle \frac{(1-{\alpha }_{l})Q}{A+{{BT}}_{c}}, & \\ \overline{\alpha }=1-\displaystyle \frac{1-{\alpha }_{h}}{1-{\alpha }_{l}}, & \zeta =\cos (\beta ), & \mathrm{and}\\ \lambda =\displaystyle \frac{\rho (A+{{BT}}_{c})}{B\omega }. & & \end{array}\end{eqnarray*}$$

These parameters have the following physical interpretations.

• 1.  
  γ: seasonal heat capacity of the system relative to the outgoing radiation over 1 yr.
• 2.  
  q: a measure of radiative forcing balance. It is directly proportional to the annual average incoming solar radiation and inversely proportional to the outgoing radiation at a critical temperature T_c .
• 3.  
  $\overline{\alpha }$: a measure of albedo contrast that changes the ice–albedo feedback. The minimum value of $\overline{\alpha }=0$ means that the high-albedo regions, α_h , and low-albedo regions, α_l , have the same albedo. The maximum value of $\overline{\alpha }=1$ is not necessarily the maximal albedo contrast but is instead the range of states where α_h = 1 (the regions with high albedo are infinitely reflective and absorb no energy) and α_l ≠ 1 (the regions with low albedo absorb some energy).
• 4.  
  λ: a measure of the speed of the ice line response to the changes in temperature.

The nondimensional parameter for heat transport in the diffusion model is

$$\begin{eqnarray*}&&\delta =\displaystyle \frac{D}{B},\end{eqnarray*}$$

and it is a measure of the efficiency of heat transport across latitudes (Stone 1978; Rose et al. 2017). The nondimensional parameter for heat transport in the relaxation to the mean model has the same form,

$$\begin{eqnarray*}&&\mu =\displaystyle \frac{C}{B},\end{eqnarray*}$$

and similar interpretation. Note that despite similar forms, μ does not directly correspond to δ. The discrepancy is caused by the fact that the diffusion coefficient in the diffusion model is not a linear scaling of the horizontal heat transfer coefficient C in relaxation to the mean model. North (1975a) commented that the divergence operator and the relaxation to the mean operator have the same effect on degree 2 polynomials when δ = μ/6. When N = 3, the relationship between these parameters is more complicated. We consider the relationship between equilibrium solutions to the two models in more detail in the Appendix. Throughout this work, we use δ = μ/6 when directly comparing the behavior of the two models for fixed values of the heat transport parameter.

Rewriting the albedo function with the nondimensionalization for ice caps yields

$$\begin{eqnarray*}{\alpha }^{* }(y,\eta )=\left\{\begin{array}{ll}1 & 0\lt y\lt \eta \\ \displaystyle \frac{2-\overline{\alpha }}{2} & y=\eta \\ 1-\overline{\alpha } & \eta \lt y\lt 1.\end{array}\right.\end{eqnarray*}$$

For belts, the positions of 1 and 1 − α are swapped with α*(y, η) = 1 for y > η and $\alpha * (y,\eta )=1-\overline{\alpha }$ for y < η.

The nondimensionalized version of the temperature model (Equation (1)) when heat transport is modeled as relaxation to the annual average temperature is given by

$$\begin{eqnarray}&&\gamma \displaystyle \frac{\partial {T}^{* }}{\partial \tau }=q{\sigma }_{2N}(y,\zeta ){\alpha }^{* }(y,\eta )-{T}^{* }-\mu ({T}^{* }-\overline{{T}^{* }}),\end{eqnarray} \tag{ 4 }$$

and when heat transport is modeled as diffusion, it is given by

$$\begin{eqnarray}&&\gamma \displaystyle \frac{\partial {T}^{* }}{\partial \tau }=q{\sigma }_{2N}(y,\zeta ){\alpha }^{* }(y,\eta )-{T}^{* }+\delta {{\rm{\nabla }}}^{2}{T}^{* }.\end{eqnarray} \tag{ 5 }$$

Note that T* is a function of y, η, and τ, but those dependencies are being suppressed in the above equations. In the relaxation to the mean model, $\overline{{T}^{* }}$ depends on η and τ (see Section 2.2 for more details). The nondimensionalized ice line equation for ice caps is

$$\begin{eqnarray}&&\displaystyle \frac{\partial \eta }{\partial \tau }=\lambda ({T}^{* }(\eta ,\tau )-1),\end{eqnarray} \tag{ 6 }$$

where T*(η, τ) is the temperature evaluated at the ice line y = η. Note that, as with the dimensional Equation (3), the right-hand side must be multiplied by −1 for ice belts. In the following, we do not explicitly consider the dynamics of the ice line, so λ plays no role in the following analysis.

2.2. Analysis of Temperature Equilibria

We expect to observe planets that have existed for a long time, and therefore we expect in the simplest case to find the corresponding planets at temperature–ice equilibrium. We focus on the equilibria of the ice line η as they undergo bifurcations in radiative forcing q and the associated hysteresis loops in our nonsmooth system. Previously, bifurcations in A and Q have been considered for Earth's range of parameter values (Widiasih 2013), and bifurcations in q have been studied in Rose et al. (2017). We follow the framework developed in Rose et al. (2017) and study bifurcations in q by considering an illustrative sample of combinations of the physical parameters. In Section 3, we will also conduct a parameter sweep in the obliquity, albedo contrast, and efficiency of heat transport to determine the likelihood of stable partial ice cover.

In order to compute the parameter values that correspond to the bifurcations, we first need to find an expression for the equilibrium temperature profile T*(y, τ). Below, we compute the equilibrium temperature profile for the relaxation to the mean model. This derivation can be found in many places (e.g., McGehee & Widiasih 2014), so we only summarize it here for convenience. Following these calculations, we summarize the method to compute the equilibrium for the diffusion model, which is shown in detail in North (1975a) and also summarized in Rose et al. (2017). We briefly compare these analytical results in the Appendix.

2.2.1. Relaxation to the Mean Equilibrium Temperature

To find ${T}_{\mathrm{eq}}^{* }(y,\eta )$, the equilibrium temperature at each latitude, we set $\tfrac{\partial {T}^{* }}{\partial t}=0$ and first find $\overline{{T}_{\mathrm{eq}}^{* }}(\eta )$, the global average equilibrium temperature, by averaging over the latitudinal range zero to 1. Since ${\int }_{0}^{1}{T}^{* }(y,\eta ,\tau ){dy}=\overline{{T}^{* }}(\eta ,\tau )$, the last term in Equation (4) cancels, and we can explicitly solve for $\overline{{T}_{\mathrm{eq}}^{* }}(\eta )$, namely,

$$\begin{eqnarray*}&&\overline{{T}_{\mathrm{eq}}^{* }}(\eta )={\int }_{0}^{1}q{\sigma }_{2N}(y,\zeta )({\alpha }^{* }(y,\eta )){dy}.\end{eqnarray*}$$

Note that T_eq(η) depends on η, since α^∗(y, η) depends on η, while the dependence on y is integrated out.

Since α^∗(y, η) differs for ice belts and ice caps, $\overline{{T}_{\mathrm{eq}}^{* }}(\eta )$ also differs. Both solutions are polynomials depending on obliquity and ice edge latitude. For ease of notation, we let

$$\begin{eqnarray*}&&{{\rm{\Sigma }}}_{2N}(\eta )={\int }_{0}^{\eta }{\sigma }_{2N}(y,\zeta ){dy}=\sum _{n\,=\,0}^{N}{a}_{2n}{p}_{2n}(\zeta ){P}_{2n}(\eta ),\end{eqnarray*}$$

where _Pi (y) = ∫_pi (y)dy is the antiderivative of the Legendre polynomial _pi (y). Then, the global average equilibrium temperature is

$$\begin{eqnarray}\overline{{T}_{\mathrm{eq}}^{\ast }}(\eta )=\left\{\begin{array}{ll}q\left((1-\overline{\alpha })+\overline{\alpha }{{\rm{\Sigma }}}_{2N}(\eta )\right), & \mathrm{ice}\ \mathrm{caps}\\ q\left(1-\overline{\alpha }{{\rm{\Sigma }}}_{2N}(\eta )\right), & \mathrm{ice}\ \mathrm{belts}.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

Note that $\overline{{T}_{\mathrm{eq}}^{\ast }}(\eta )$ is proportional to the nondimensionalized radiative forcing q. For fixed η, the more radiative forcing the planet receives, the warmer its mean equilibrium temperature.

Once $\overline{{T}_{\mathrm{eq}}^{\ast }}(\eta )$ has been found, we may solve for ${T}_{\mathrm{eq}}^{\ast }(y,\eta )$ in Equation (4) with $\tfrac{\partial {T}^{\ast }}{\partial t}=0$. Indeed, for y ≠ η,

$$\begin{eqnarray*}&&{T}_{\mathrm{eq}}^{\ast }(y,\eta )=\displaystyle \frac{1}{1+\mu }\left(q{\sigma }_{2N}(y,\zeta ){\alpha }^{\ast }(y,\eta )+\mu \overline{{T}_{\mathrm{eq}}^{\ast }}(\eta )\right).\end{eqnarray*}$$

Due to the discontinuity in α*(y, η), the temperature profile is discontinuous at the ice line. We define the value at the ice line to be the average of the left and right limits of ${T}_{\mathrm{eq}}^{* }$ as y approaches η, namely,

$$\begin{eqnarray*}&&{T}_{\mathrm{eq}}^{* }(\eta ,\eta )=\displaystyle \frac{\mathop{\mathrm{lim}}\limits_{y\to {\eta }^{+}}{T}_{\mathrm{eq}}^{* }(y,\eta )+\mathop{\mathrm{lim}}\limits_{y\to {\eta }^{-}}{T}_{\mathrm{eq}}^{* }(y,\eta )}{2}.\end{eqnarray*}$$

For ice caps, the left- and right-hand limits are

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{y\to {\eta }^{-}}{T}_{\mathrm{eq}}^{* }(y,\eta )=\displaystyle \frac{q{\sigma }_{2N}(\eta ,\zeta )+\mu \overline{{T}_{\mathrm{eq}}^{* }}(\eta )}{1+\mu },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{y\to {\eta }^{+}}{T}_{\mathrm{eq}}^{* }(y,\eta )=\displaystyle \frac{q{\sigma }_{2N}(\eta ,\zeta )(1-\overline{\alpha })+\mu \overline{{T}_{\mathrm{eq}}^{* }}(\eta )}{1+\mu },\end{eqnarray} \tag{ 9 }$$

and for ice belts, ${\mathrm{lim}}_{y\to {\eta }^{-}}$ and ${\mathrm{lim}}_{y\to {\eta }^{+}}$ are swapped. In both cases, the equilibrium temperature at the ice line is given by

$$\begin{eqnarray*}&&{T}_{\mathrm{eq}}^{* }(\eta ,\eta )=\displaystyle \frac{q{\sigma }_{2N}(\eta ,\zeta )(2-\overline{\alpha })+2\mu \overline{{T}_{\mathrm{eq}}^{* }}(\eta )}{2(1+\mu )};\end{eqnarray*}$$

however, the function for $\overline{{T}_{\mathrm{eq}}^{* }}(\eta )$ differs for ice caps and ice belts (see Equation (7)). Note that this definition of the temperature profile at the ice line coincides with the result derived using our definition of α^∗(η, η).

Ice line equilibria occur when T_eq*(η, η) = 1. We consider the response of equilibria to changes in radiative forcing q. The relaxation to the mean model allows us to solve exactly for the unique value of radiative forcing q as a function of the other parameters, namely,

$$\begin{eqnarray}&&{q}_{\eta }(\zeta ,\overline{\alpha },\mu )=\displaystyle \frac{2(1+\mu )}{{\sigma }_{2N}(\eta ,\zeta )(2-\overline{\alpha })+2\mu {T}_{x}(\eta ,\zeta )},\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray*}{T}_{x}(\eta ,\zeta )=\left\{\begin{array}{ll}(1-\overline{\alpha })+\overline{\alpha }{{\rm{\Sigma }}}_{2N}(\eta ,\zeta ) & \mathrm{ice}\ \mathrm{caps}\\ 1-\overline{\alpha }{{\rm{\Sigma }}}_{2N}(\eta ,\zeta ) & \mathrm{ice}\ \mathrm{belts}.\end{array}\right.\end{eqnarray*}$$

Note that ${T}_{x}(\eta ,\zeta )=\overline{{T}_{\mathrm{eq}}^{* }}(\eta )/q$.

For ice caps, a planet is ice-free when η = 1 and in a snowball state when η = 0. Following the conventions from Rose et al. (2017), for ice caps, we let q_free denote the lowest value of q for which q₁ is stable. Stability of the ice-free state is inferred from where q_η intersects the line η = 1. When q > q₁, the nondimensional temperature ${T}_{\mathrm{eq}}^{* }$ at the pole is greater than 1. Similarly, we let q_snow denote the location where q_η intersects the line η = 0. We can think of this as the highest value of q for which q₀ is stable. Note that for ice belts, q₁ is q_snow, and q₀ is q_free. See Figure 2 for a bifurcation diagram demonstrating q_free and q_snow. It should be noted that the definitions of q_free and q_snow used in this work take the average of the two sides of the discontinuity at the ice line. For ice caps, physically, the definition q_free = q₁ can be interpreted as vanishingly small ice caps, and q_snow = q₀ can be interpreted as a vanishingly small equatorial strip of water. Alternative definitions can be used. We elaborate on the implications of this choice in Section 5.

Zoom In Zoom Out Reset image size

2.2.2. Diffusion Equilibrium Temperature

For the remainder of the paper, we will use a tilde over a function or variable to denote that it is associated with the diffusion version of the model.

North (1975a) analytically solved the diffusive energy balance Equation (5) with ∂T/∂t = 0 for Earth with a degree 2 approximation to the insolation distribution. He found that solutions to the diffusive equation can be expressed in terms of hypergeometric functions. This analytical solution is generalized to the exoplanet case with degree 2 approximation to the insolation distribution in Rose et al. (2017).

For an arbitrary approximation to the insolation approximation, one may use the same methods described in North (1975a) and Rose et al. (2017). For ease of notation, let

$$\begin{eqnarray*}&&{H}_{2N}(x)=\sum _{n=0}^{N}\displaystyle \frac{{a}_{2n}{p}_{2n}(\zeta )}{1\,+\,2n(2n+1)\delta }{p}_{2n}(x).\end{eqnarray*}$$

As noted in North (1975a), qH_2N (x) is the particular solution to the diffusion equation on the interval [0, η), and $q(1-\overline{\alpha }){H}_{2N}(x)$ is the particular solution to the diffusion equation on the interval (η, 1]. To find the general solution, one must find the solution to the homogeneous equation. This derivation is given in North (1975a), so we do not give it here. Instead, we report the results of those computations that yield the equilibrium temperature at the ice line,

$$\begin{eqnarray*}{\tilde{T}}_{\mathrm{eq}}^{* }\left(\eta \right)=\left\{\begin{array}{ll}q\left({H}_{2N}\left(\eta \right)-\overline{\alpha }{F}_{2N}\left(\eta \right)\right), & \mathrm{ice}\ \mathrm{caps}\\ q\left(\left(1-\overline{\alpha }\right){H}_{2N}\left(\eta \right)+\overline{\alpha }{F}_{2N}\left(\eta \right)\right), & \mathrm{ice}\ \mathrm{belts}\end{array}\right.,\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{F}_{2N}(\eta ) & = & \displaystyle \frac{P{{\prime} }_{\nu }(\eta ){H}_{2N}(\eta )-H{{\prime} }_{2N}(\eta ){P}_{\nu }(\eta )}{P{{\prime} }_{\nu }(\eta ){f}_{\nu }(\eta )-f{{\prime} }_{\nu }(\eta ){P}_{\nu }(\eta )}{f}_{\nu }(\eta ),\\ {f}_{\nu }(x) & = & {\ }_{2}{F}_{1}(-\nu /2,(1+\nu )/2,1/2,{x}^{2}),\\ {P}_{\nu }(x) & = & {\ }_{2}{F}_{1}((1+\nu )/2,-\nu /2,1,1-{x}^{2}),\\ \nu & = & -\displaystyle \frac{1}{2}+\displaystyle \frac{{\left(1-4/\delta \right)}^{1/2}}{2},\end{array}\end{eqnarray*}$$

and ₂ F₁ is the hypergeometric function. Although P_ν (x) and f_ν (x) are complex-valued for δ < 4, taking the real parts of these functions yields linearly independent solutions, as noted in North (1975a). Solving for q gives

$$\begin{eqnarray*}{\tilde{q}}_{\eta }{\left(\zeta ,\overline{\alpha },\delta \right)}^{-1}=\left\{\begin{array}{ll}{H}_{2N}(\eta )-\overline{\alpha }{F}_{2N}(\eta ), & \mathrm{ice}\ \mathrm{caps}\\ (1-\overline{\alpha }){H}_{2N}(\eta )+\overline{\alpha }{F}_{2N}(\eta ), & \mathrm{ice}\ \mathrm{belts}.\end{array}\right.\end{eqnarray*}$$

When N = 1, the results from Rose et al. (2017) are recovered.

3.1. Defining the Region of Integration

Planets with stable partial ice cover are potential candidates for a snowball catastrophe bifurcation. We focus on quantifying the likelihood of stable partial ice cover depending on planetary obliquity. We follow Rose et al. (2017) and compute an estimate of the likelihood of stable partial ice cover based on the size of the region in parameter space that admits these types of solutions. Locations where q_η or ${\tilde{q}}_{\eta }$ have a zero derivative demarcate regions of stable partial ice cover because the same location in the transformed plot with η on the vertical axis is a saddle node bifurcation point (Figure 2).

Taking the derivative of q with respect to η yields

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {q}_{\eta }}{\partial \eta } & = & \displaystyle \frac{-2(1+\mu )}{{\left[(2-\overline{\alpha }){\sigma }_{2N}(\eta ,\zeta )+2\mu {T}_{x}(\eta )\right]}^{2}}\\ & & \times \,\left(\left(2-\overline{\alpha }\right)\displaystyle \frac{\partial }{\partial \eta }{\sigma }_{2N}\left(\eta ,\zeta \right)+2\mu \left(\pm \overline{\alpha }{\sigma }_{2N}\left(\eta ,\zeta \right)\right)\right)\end{array}\end{eqnarray} \tag{ 11 }$$

for caps and belts, where $\pm \overline{\alpha }{\sigma }_{2N}(\eta ,\zeta )$ is the derivative of T_x (η, ζ) for caps and belts, respectively. Since T_x (η) ≥ 0, this quantity is well defined. Parameter values that result in $\tfrac{\partial {q}_{\eta }}{\partial \eta }=0$ give the location of the saddle node. Setting Equation (11) to zero and solving for $\overline{\alpha }$ yields

$$\begin{eqnarray}{\alpha }_{\mathrm{crit}}(\zeta ,\mu ,\eta )=\left\{\begin{array}{ll}\displaystyle \frac{2\tfrac{\partial }{\partial \eta }{\sigma }_{2N}(\eta ,\zeta )}{\tfrac{\partial }{\partial \eta }{\sigma }_{2N}(\eta ,\zeta )-2\mu {\sigma }_{2N}(\eta ,\zeta )} & \ {\rm{for}}\,{\rm{caps,}}\\ \displaystyle \frac{2\tfrac{\partial }{\partial \eta }{\sigma }_{2N}(\eta ,\zeta )}{\tfrac{\partial }{\partial \eta }{\sigma }_{2N}(\eta ,\zeta )+2\mu {\sigma }_{2N}(\eta ,\zeta )} & \ {\rm{for}}\,{\rm{belts}}\end{array}\right.,\end{eqnarray} \tag{ 12 }$$

which is the critical value of the albedo contrast at the saddle node bifurcation latitude. Stable partial ice cover is possible whenever $\overline{\alpha }\lt {\alpha }_{\mathrm{crit}}(\zeta ,\mu ,\eta )$. At $\overline{\alpha }={\alpha }_{\mathrm{crit}}(\zeta ,\mu ,\eta )$, the saddle node bifurcation occurs. Note that since μ is theoretically unbounded, the function α_crit can become arbitrarily small. In the next section, we use α_crit to find the relative likelihood of stable partial ice cover.

Following the method in Rose et al. (2017) but with the general degree 2N insolation approximation, we find that in the diffusion model,

$$\begin{eqnarray}{\tilde{\alpha }}_{\mathrm{crit}}(\zeta ,\mu ,\eta )=\left\{\begin{array}{ll}\displaystyle \frac{H{{\prime} }_{2N}(\eta )}{F{{\prime} }_{2N}(\eta )} & {\rm{for}}\,{\rm{caps,}}\\ \displaystyle \frac{H{{\prime} }_{2N}(\eta )}{H{{\prime} }_{2N}(\eta )-F{{\prime} }_{2N}(\eta )} & {\rm{for}}\,{\rm{belts}}.\end{array}\right.\end{eqnarray} \tag{ 13 }$$

Recall that H_2N and F_2N both also depend on ζ and δ.

3.2. Relative Likelihood of Stable Partial Ice Cover

In this section, we integrate over the parameter region in order to estimate the relative likelihood of stable partial ice cover using the same method as was described in Rose et al. (2017). We summarize it here for convenience.

We assume that there exists a true distribution for each parameter and that the parameters are independent of each other. However, since the true probability distributions are not known, we consider a handful of candidate probability distributions. We integrate the composite probability density function (PDF) over the region of the domain where stable edges are present to obtain the likelihood of stable partial ice cover for a given value of obliquity ζ. We normalize our results by the likelihood value for the obliquity of the Earth, ζ = $\cos$ (23.5°). From the independence assumption, we can write the overall PDF h_planet as follows:

$$\begin{eqnarray*}{h}_{\mathrm{planet}}=\left\{\begin{array}{ll}{h}_{q}(q){h}_{\mu }(\mu ){h}_{\overline{\alpha }}(\overline{\alpha }), & \mathrm{relaxation}\ \mathrm{to}\ \mathrm{mean},\\ {h}_{q}(q){h}_{\delta }(\delta ){h}_{\overline{\alpha }}(\overline{\alpha }), & \mathrm{diffusion}.\end{array}\right.\end{eqnarray*}$$

We follow Rose et al. (2017) in our choice of candidate probability functions to test in order to facilitate the comparison with their results and include two additional probability functions to account for the difference in δ and μ. Since q, μ, and δ are nonnegative and unbounded, lognormal distributions are used to incorporate the possibility of a long tail. The general form of the PDF of lognormal distribution is

$$\begin{eqnarray*}&&h(x)=\displaystyle \frac{1}{\sigma \sqrt{2\pi }(x-\theta )}\exp \left(\displaystyle \frac{-{\left(\mathrm{ln}\left(\tfrac{x-\theta }{m}\right)\right)}^{2}}{2{\sigma }^{2}}\right),\end{eqnarray*}$$

where σ is referred to as the shape parameter, m is the scale parameter, and θ is the location parameter. In our two new PDFs, we use gamma distributions for μ. As with the lognormal distribution, the gamma distribution is a maximum-entropy probability distribution and so minimizes the amount of prior information included in the distribution. The use of the gamma distribution allows us to create PDFs that make Earth's value of μ = 1.6. While the gamma distribution has similar properties to those of the lognormal distribution, they are different enough to make them good candidates for sensitivity analysis (Wiens 1999; Iaci 2000). The general form of the PDF for the gamma distribution is

$$\begin{eqnarray*}&&h(x)=\displaystyle \frac{{b}^{a}{x}^{a-1}{e}^{-{bx}}}{{\rm{\Gamma }}(a)},\end{eqnarray*}$$

where a is the shape parameter, b is the rate parameter, and Γ is the gamma function. Since $\overline{\alpha }\in [0,1]$, uniform and beta distributions are used. The beta distribution has the PDF

$$\begin{eqnarray*}&&{h}_{\overline{\alpha }}(\overline{\alpha })=6\overline{\alpha }(1-\overline{\alpha })\end{eqnarray*}$$

and favors values of $\overline{\alpha }$ close to the value for Earth ($\overline{\alpha }=0.44$) compared to extremes of $\overline{\alpha }=0$ or 1. The PDF parameters are chosen so that Earth's parameter values are not unlikely.

For PDF0, h_α is uniform on [0, 1]; h_μ and h_δ are lognormal on [0, ∞ ] with shape parameter 1.0, scale parameter 1.0, and location parameter 0; and h_q is lognormal on [0, ∞ ] with shape parameter 0.5, scale parameter 1.0, and location parameter 0.

The PDF1 is the same as PDF0 except h_μ and h_δ are lognormal on [0, ∞ ] with shape parameter 2.0, scale parameter e, and location parameter 0. Compared to PDF0 and PDF2, PDF1 makes very small and large values of μ or δ more likely.

The PDF2 is the same as PDF0 except h_α is a parabolic beta distribution on [0, 1] with mode at 0.5. Compared to PDF0 and PDF1, PDF2 favors intermediate values of albedo contrast compared to extreme values.

The PDF3 is the same as PDF0 except h_μ and h_δ have a gamma distribution on [0, ∞ ] with shape parameter 4.0 and rate parameter 2.5. Compared to PDF0, PDF3 makes the μ value for Earth (1.6) more likely.

The PDF4 is the same as PDF1 except h_μ and h_δ have a gamma distribution on [0, ∞ ] with shape parameter 4.0 and rate parameter 1.85. Compared to PDF3, PDF4 makes larger values of μ more likely.

The (nonnormalized) likelihood is given by

$$\begin{eqnarray}&&{P}_{\mathrm{ice}}(\beta )=\displaystyle \frac{{\int }_{0}^{1}{\int }_{0}^{\infty }{\int }_{0}^{{\alpha }_{\mathrm{crit}}}{h}_{\mathrm{planet}}({q}_{\eta },\mu ,\alpha )d\alpha \ d\mu \ d\eta }{{\int }_{0}^{\infty }{\int }_{0}^{\infty }{\int }_{0}^{1}{h}_{\mathrm{planet}}(q,\mu ,\alpha )d\alpha \ d\mu \ {dq}}.\end{eqnarray} \tag{ 14 }$$

Note that the denominator in Equation (14) is equal to 1, since h_planet is a PDF. For higher values of μ and δ, the corresponding values of α_crit become vanishingly small for all latitudes η. Therefore, even though the region of integration is unbounded, numerical integration converges. The results of the integration are summarized in Figure 3.

Zoom In Zoom Out Reset image size

The qualitative results are similar across all models. Planets at low obliquity have a higher likelihood of stable partial ice cover than planets at high obliquity, and planets with moderate obliquity have the lowest likelihood of stable partial ice cover. The apparent discontinuities in the relative-likelihood plots for the relaxation to the mean model in Figure 3 are due to cutoffs for the obliquities where we use a model with caps or belts. As shown by Dobrovolskis (2021), the minima of the insolation distribution occur between the poles and the equator when the obliquity is between approximately 45° and 65.355°. This behavior can be thought of as a transitional regime between ice caps (where the insolation minima occur at the poles) and ice belts (where the insolation minimum occurs at the equator). While we do not specifically consider planets with varying obliquity in this work, we use both the cap and belt models for this range of obliquities.

In agreement with the degree 2 diffusion model in Rose et al. (2017), the likelihood of stable partial ice cover goes to zero for the degree 2 relaxation to the mean model. It is easy to see the reason for this by noting that the derivative of the insolation function appears in the numerator of α_crit and recalling that the insolation function σ₂ is constant for $\zeta =\sqrt{3}/3$, which is when the obliquity is 54.74°. For the models with higher-degree insolation approximations, the likelihood is never identically zero but nevertheless attains its minimum at mid-obliquities for all tested probability distribution functions.

For example, in Figure 4, we plot bifurcation diagrams for the diffusion model at β = 54.47°. For N = 3 shown in the figure, small values of $\overline{\alpha }$ and δ admit stable partial ice edges for both ice caps and ice belts, although it is only for the smallest values of these parameters that the edges are accessible in a hysteresis loop. Once α ≈ 0.3, there are no longer any stable partial ice edges in the diffusion model. In contrast, no value of albedo contrast or efficiency of heat transport will cause the diagram for N = 1 to have a saddle node bifurcation.

Zoom In Zoom Out Reset image size

In addition to assessing the likelihood of partial ice cover, we quantify the severity of the snowball state. Typically, one is interested in determining the inner and outer edges of the habitable zone based on the system parameters; however, as Rose et al. (2017) noted in their Section 5.3 (and references therein), the models are too simplistic to give good estimates for this range. It is possible to quantify whether a bifurcation from no or partial ice cover to the snowball state occurs for the system parameters.

The most severe snowball bifurcation occurs when q_free < q_snow (similarly, ${\tilde{q}}_{\mathrm{free}}\lt {\tilde{q}}_{\mathrm{snow}}$) and for all 0 < η < 1, q_free ≤ q_η ≤ q_snow. Although there may be stable ice line equilibria between zero and 1, they would inaccessible by varying q through a hysteresis loop in this situation. A less severe bifurcation occurs when the ice line continuously transitions from the ice-free state to small stable ice caps (or ice belts), and then the ice line drops to the snowball state. Behavior of solutions is consistent with the typical passage through a saddle node bifurcation (Strogatz 2018). The least severe scenario is where no saddle node bifurcation occurs and the stable ice line equilibrium depends continuously on the parameter q. This occurs when q_snow < q_free and for all 0 < η < 1, q_snow ≤ q_η ≤ q_free.

In Figure 5, we plot contours for the smallest value of the nondimensional incoming radiation q for which the snowball state is not globally attracting (i.e., the only stable solution). In Rose et al. (2017), this value was referred to as q_hab and computed by

$$\begin{eqnarray*}&&{q}_{\mathrm{hab}}=\min \{{q}_{\mathrm{free}},{q}_{\mathrm{stab}}\},\end{eqnarray*}$$

where q_stab is the minimum q value for which an ice line equilibrium exists between zero and 1. Frequently, this corresponds to the saddle node bifurcation, which causes the snowball catastrophe (if present; e.g., in Figure 4, only the lower left panel exhibits this bifurcation). If the model has stable partial ice edges between zero and 1 but no saddle node bifurcation causing the snowball catastrophe, q_stab = q_snow. When q_free > q_stab, the hysteresis loop generated by varying q causes a drop from the ice-free state to a partially ice-covered state instead of directly to the snowball state.

Zoom In Zoom Out Reset image size

In Figure 5, the regions above the dashed black curves are where the snowball catastrophe occurs directly from the ice-free state. In an extension of the work in Rose et al. (2017), here we consider the models with the degree 6 approximation to the insolation distribution instead of the degree 2 approximation for insolation in the diffusion model. For ease of comparison with the results in Rose et al. (2017), we restrict the ice cap models to obliquities between 0° and 55° and the ice belt models to obliquities between 55° and 90°.

It is straightforward to compare the top row of Figure 5 to Figure 8 in Rose et al. (2017). Analytically, the only difference is the degree of the insolation distribution used. On the whole, the figures are very similar, having nearly the same contours. They differ only slightly in the maximum and minimum values; the sixth-degree approximation decreases the extrema achieved for q_hab for the model with ice caps and increases the extrema for the model with an ice belt relative to the second-degree approximation. This is due to the increased accuracy of the insolation approximation at the poles (minima for ice caps, maxima for an ice belt). The largest difference is that for the smallest obliquities in the degree 6 diffusion model, the snowball catastrophe from the ice-free state occurs for larger values of the heat transport parameter. This means the snowball catastrophe is less severe for low-obliquity planets than in the degree 2 diffusion model. The transition from ice-free directly to the snowball state is minimally affected for obliquities above approximately 30°.

In the bottom row of Figure 5, we plot the contours for the degree 6 relaxation to the mean model. The relaxation to the mean model has stark contrasts with the diffusion model in the row above it. For low-to-moderate albedo contrast, the minima and maxima of both models are relatively similar. However, for large albedo contrast, the maximum values achieved for the relaxation to the mean model are significantly greater. Analytically, this is caused by the dependence of q_free on the albedo contrast. In the diffusion model, q_free is independent of $\overline{\alpha }$, but in the relaxation to the mean model, it is not. As $\overline{\alpha }$ increases, so does q_free for the relaxation to the mean model (this can be seen in Equation (10)).

In both models, increasing the albedo contrast increases the severity of the snowball catastrophe for all obliquities (the dashed black curves decrease), but the increase in severity in the relaxation to the mean model is minimal. However, for all obliquities in either model, lower heat transport efficiency always guards against the most severe transition (from ice-free directly to snowball). In the diffusion model for the range of heat transport plotted, this most severe transition is always present for the mid-obliquities and increasingly present for higher obliquities as the albedo contrast increases. To avoid the ice-free directly to snowball transition requires decreasing the heat transport efficiency at least another order of magnitude than what is plotted (i.e., to at least 0.001 for β = 90° and $\overline{\alpha }=0.7$). Conversely, in the relaxation to the mean model, there is only a small range of obliquities where this most severe transition is present for the range of heat transport considered.

In this paper, we have analyzed a one-dimensional energy balance model with heat transport modeled by relaxation to the global mean temperature and diffusion. The relaxation to the mean and diffusion versions of the Budyko–Sellers model provide two similar ways to model the multitude of processes involved in energy transfer between latitudes. Both methods ensure that energy is transported from latitudes that are "hot" to ones that are "cold." The diffusive heat transport is a local process that necessitates special treatment at the poles, while relaxation to the mean global temperature is a global process that does not require special boundary conditions (Widiasih 2013). This work may be interpreted for any rapidly rotating rocky planet with some physical mechanism of heat transport. This work applies to planets where the temperature affects the albedo; in particular, we assume that higher temperatures decrease the albedo as they do for ice/water on Earth.

We have also considered the effects of different approximations to the annual average insolation distribution on the results of the models. The second-degree approximation of the insolation distribution used in Rose et al. (2017) does not capture the qualitative distribution of mid-obliquity planets. Planets with obliquities between approximately 45° and 65° have a characteristic "W" shape that requires a degree 6 (or higher) polynomial approximation (Nadeau & McGehee 2017).

A main result from Rose et al. (2017) is that the likelihood of stable partial ice cover goes to zero at 55° obliquity. This result is due to the insolation approximation used in their study, which is constant for 55° obliquity. In the definition of ${\tilde{\alpha }}_{\mathrm{crit}}$ (Equation (13)), if H_2N is constant, then ${\tilde{\alpha }}_{\mathrm{crit}}=0$ for all values of the arguments. When N = 1, as in Rose et al. (2017), the function H_2N is constant exactly when the insolation is constant, i.e., when the obliquity is 54.47°. This means that the integral in the numerator of the likelihood calculation (Equation (14)) is zero. We see a similar problem with the relaxation to the mean model when N = 1. Here it is clear to see that when the insolation is constant, α_crit = 0 (see Equation (12)). Taking a higher degree approximation in either model, as we do here, avoids these problems.

In Figure 3, we see that low-obliquity planets are more likely to have stable partial ice cover than those with high obliquity, and planets at middle obliquities are least likely to have stable partial ice cover, which is qualitatively similar to the likelihood computations in Rose et al. (2017). As noted above, we find that stable partial ice cover is possible at all obliquities and that, in particular, the relative likelihood of finding a planet with partial stable ice cover ranges between 7% and 32%.

Comparing the relaxation to the mean model to the diffusion model, we note that the latter predicts a lower likelihood of stable partial ice cover at lower obliquities. This effect is due to the fact that the diffusion model can exhibit a second saddle node bifurcation at high values of η close to the poles (called the small ice cap instability) for Earth's obliquity. As the obliquity decreases to zero, the small ice cap instability shrinks in size and can disappear completely for small-to-moderate values of the albedo contrast $\overline{\alpha }$, increasing the range of η where stable ice edges are possible. In contrast, the relaxation to the mean model does not exhibit the small ice cap instability at Earth's obliquity for degree 2 or degree 6 insolation approximation. This means that the likelihood remains relatively flat as the obliquity decreases.

The relaxation to the mean model also exhibits pronounced differences between PDF0, PDF1, and PDF2 at high obliquities. For high values of the obliquity β, the likelihood of stable partial ice cover is lower for PDF2. The difference between PDF2 and other tested PDFs is due to the differences in α_crit between the relaxation to the mean model and the diffusion model. Since PDF2 changes the distribution of $\overline{\alpha }$ from a uniform distribution to a parabolic beta distribution, the shape of α_crit results in a more pronounced difference for the relaxation to the mean model than for the diffusion model. The resulting gap between the likelihood curves conveys decreased certainty about the likelihood of stable partial ice cover on high-obliquity planets. Note also that for PDF3 and PDF4, the relaxation to the mean model exhibits behavior similar to that for PDF1.

We quantify the effects of albedo contrast and efficiency of heat transport on the presence of hysteresis loops in radiative forcing. We find that the severity of the snowball bifurcation increases as the albedo contrast $\overline{\alpha }$ and the efficiency of heat transport (δ or μ) increase.

The above behavior can be explained by the effect of albedo contrast and efficiency of heat transport on the planetary climate mechanism. When the albedo contrast $\overline{\alpha }$ is low, ice is not much more reflective than the nonfrozen regions (either ground or water), resulting in suppressed ice–albedo feedback. When δ or μ is low, the near absence of heat transport across latitudes limits the interaction between the ice regions and the water regions of the planet, thus reducing the likelihood of the snowball catastrophe. When $\overline{\alpha }$ is high, the reflectivity of ice is much higher than that of the nonfrozen regions, expediting the ice–albedo processes. When μ is high, the heat transport across latitudes makes it difficult to maintain a difference in temperatures between ice regions and water regions, leading to an ice-free or a snowball planet.

In the range of obliquities where both ice caps and ice belts may be stable, it is not possible to transition continuously between the ice-free and snowball states, as q is varied; i.e., there will always be a hysteresis loop when varying q. The hysteresis will contain either a saddle node bifurcation or the most severe snowball bifurcation from ice-free to completely ice-covered. A future study might explore whether the lack of a region without hysteresis in the parameter space is a contributing factor for the decrease in the likelihood of stable partial ice cover for these obliquities in Figure 3.

It should be noted that dealing with the discontinuity at the ice line affects the definitions of q_free and q_snow. For ice caps, physically, the definition q_free = q₁ can be interpreted as vanishingly small ice caps, and q_snow = q₀ can be interpreted as a vanishingly small equatorial strip of water. This choice is consistent with the literature on the relaxation to the mean model (Widiasih 2013; McGehee & Widiasih 2014). This definition ensures that the bifurcation diagram for the relaxation to the mean model is continuous at η = 0 and 1. However, one could consider an alternative definition of q_free and q_snow; instead of taking the average of the two sides of the discontinuity at the ice line, one could use only the warm branch (Equation (8)) of the temperature equilibrium to define q_free and only the cold branch (Equation (9)) to define q_snow. This choice yields ${q}_{\mathrm{free}}=\tfrac{(1+\mu )}{{\sigma }_{2N}(\eta ,\zeta )+\mu }$ and ${q}_{\mathrm{snow}}=\tfrac{(1+\mu )}{(1-\overline{\alpha })({\sigma }_{2N}(\eta ,\zeta )+\mu )}$. Such a definition might be more physically intuitive, representing the relevant thresholds for the first appearance of freezing temperatures at the poles for a warm planet undergoing a gradual cooling process and the first appearance of thawing temperatures for a planet undergoing a gradual thawing process. Removing the dependence of q_free on $\overline{\alpha }$ might lead to a higher estimate for the likelihood of ice-free solutions. However, analyzing the implications of these alternatives would require extensive mathematical analysis of the bistability of ice-free and partial ice states that we defer to future work.

The robustness of the snowball catastrophe and the parameter regimes where an energy balance model might be applicable has been debated. In the global circulation model (GCM) simulations conducted by Ferreira et al. (2014), the snowball catastrophe occurs only for particular ocean regimes. Wagner & Eisenman (2015) showed that meridional heat transfer may increase ice cover stability. Rose & Marshall (2009) extended the energy balance models to include ocean heat transport and meridional structure and found that the snowball catastrophe is possible in those models.

If this work were to be applied to an observed planet, the obliquity ζ and albedo contrast $\overline{\alpha }$ could perhaps be measured, and albedo signatures could be compared to the predictions of our model. The parameter q will be more difficult to estimate. While the mean annual insolation Q can be derived from information about the star and orbit of the planet, the dependence of q on the atmospheric parameters A and B makes determining q challenging. The parameters δ and μ, the efficiency of heat transport, would also be difficult to measure directly because of our limited knowledge of the rates of heat transport on different planets.

In this paper, we have analyzed a one-dimensional energy balance model with two different methods of modeling the heat transport. The models have an explicit dependence on the planet's obliquity through the annual average insolation distribution, which we approximate with degree 2 and degree 6 polynomials. We pay particular attention to the planet's obliquity, radiative forcing, albedo contrast, and efficiency of heat transport and find the following.

• 1.  
  With an improved approximation to the insolation distribution function, planets at all values of obliquity exhibit a nonzero likelihood of stable partial ice cover regardless of the mode of heat transport. The minimum likelihood ranges from 7% to 32% relative to Earth's likelihood of stable partial ice cover, depending on the model and PDFs for the parameters. The maximum likelihood ranges from 110% to 140% relative to Earth.
• 2.  
  Several results from the earlier study by Rose et al. (2017) are seen here in both the diffusion and relaxation to the mean models.

  • (a)  
    Models with high obliquity are less likely to have stable partial ice cover than ones with low obliquity but more likely than models with moderate obliquity.
  • (b)  
    Low albedo contrast and low efficiency of heat transport favor stable partial ice cover.
  • (c)  
    High albedo contrast and high efficiency of heat transport favor a severe snowball catastrophe in a hysteresis loop caused by changes in radiative forcing.
• 3.  
  In the relaxation to the mean model, a larger range of heat transport efficiency supports partial stable ice cover at high obliquities, even for high albedo contrast. However, in the diffusion model, stable partial ice cover a high obliquities and high albedo contrast is possible only for very low heat transport efficiency.

Both models discussed in this work are highly simplified, so direct implications for real physical systems are tenuous when taken individually. We may interpret the results with higher confidence in regions of parameter space where their solutions give similar results. Places where they differ indicate a need for further investigations with more complex models.

For example, due to the disagreement of the models shown in Figure 5, a future study may explore the transition to the snowball state for high-obliquity planets in a general circulation model. Other extensions of this work would be to look at the climatic histories of planets by introducing obliquity or eccentricity variations in time into the energy balance model used in this study. Variations in a planet's orbital parameters change the behavior of the planetary ice cover over geological time and are easiest to model over long time periods with simple energy balance models like the ones presented here.

The authors thank Nikole Lewis, Tiffany Kataria, Toby Ault, and Steven Strogatz for their comments in preparing the manuscript. The authors also thank Brian Rose and an anonymous reviewer for their comments during the review process, which greatly improved the manuscript. A.N. was supported by a Mathematical Sciences Postdoctoral Research Fellowship (award No. DMS-1902887) during this project.

In this section, we discuss a brief comparison of the analytical equilibrium solutions. As mentioned in Section 2.1, North (1975a) commented that the divergence operator and the relaxation to the mean operator have the same effect on degree 2 polynomials when δ = μ/6. These operators are the same on the function subspace of degree 2 polynomials; however, solutions to the models live in larger spaces of functions, sometimes called solution spaces. In particular, the equilibrium temperature solutions must be in the solution space, so the relaxation model solution space contains piecewise polynomials of degree 2N, and the diffusion model solution space contains bounded differentiable hypergeometric functions on the interval [0, 1]. Degree 2 polynomials on [0, 1] are a subset of both spaces. Even though solutions to the models are not restricted to the function subspace of degree 2 polynomials, the relation between δ and μ of δ = μ/6 remains a remarkably good approximation for comparing solutions between the models. We find that in the small and large limits of the heat transport efficiency parameter with other parameters fixed, the equilibrium temperature solutions from both models approach the same functions.

As μ and δ approach zero, solutions of the relaxation to the mean model approach solutions to the diffusion model. This can be seen directly by considering Equations (4) and (5). As the parameters μ and δ approach zero, heat transport becomes negligible, and the models become an energy balance equation, where the latitude is simply a parameter. Convergence of the solutions across all latitudes as μ and δ approach zero appears to be rather slow due to the continuous diffusive solution approaching the jump discontinuity in the relaxation to the mean solution at the ice line.

This behavior is illustrated in Figure 6, where we compare the equilibrium temperature profile solutions of the two models with N = 3 and Earth values for all parameters except δ and μ.

Zoom In Zoom Out Reset image size

As μ and δ become large, solutions of the relaxation to the mean model again approach solutions to the diffusion model. Intuitively, this occurs because as the transport coefficients become large, heat is redistributed almost instantaneously, and the solutions approach the global mean temperature (see Equation (7)). Numerically, we see that the solutions approach each other faster than they approach the global mean. Convergence for large μ and δ appears to be much faster compared to when the parameters approach zero. For instance, for the parameters used in Figure 6, the solutions are within 3% of each other for all latitudes when δ = μ/6 = 1 and within 1% of each other for all latitudes when δ = μ/6 = 10. Further mathematical investigation of these properties is a natural direction for future work.