Climate and snowball earth (3 days)

( previous, home, next )


Hypothetical snowball earth

Hypothetical snowball earth

One of the great contemporary scientific efforts of today is the study of the climate of planet Earth. Climate fundamentally effects how we live. Life in Pennsylvania is different from life in Los Angeles or life in Miami. And the climate today is not seem to be the same as it was just a few thousand years ago. Where people were buried in boats in Tibet 4,000 years ago and paint pictures of swimming in Egypt 10,000 years ago, there is now only dessert. Geologists have discovered that before that, the earth has been through sporadic ice ages where much larger portions of the norther and southern latitudes have been covered in glaciers. And over the last two hundred years, our population and industry have grown so large that we ourselves are changing earth's climate.

The term "climate" is used to represent the aggregation of all the earth's chaotic weather conditions that are varying in time and space. This is not the same as the weather we experience each day -- no one person can see the entire state of the earth's climate. Each of us can observe and experience only a thin, insignificant slice of the climate. Without extensive personal experience, it is hard for us to develop intuition and foresight about changing climatic conditions. Modelling is one approach that can help us develop that intuition. In this chapter, we will explore two simple models of earth's climate. The first model attempts to predict the average temperature of the earth using only geometry and one well-studied law from thermodynamics. The second model is a generalization that supports a theory claiming the earth once froze solid and stayed that way for millions of years -- a snowball earth. Along the way, we'll come up with a new kind of mathematical model and some ideas for studying it.

Energy-balance and global temperature

One approach to modelling climate is to make use of the principle of conservation of energy. The idea is that the climate of Earth depends on the amount of energy in our atmosphere. One way to frame this is to say the temperature we experience here at the earth's surface is the consequence of energy entering and leaving the atmosphere through four sources -- the sun, the sky, the oceans, and the earth's crust. Intuition and order-of-magnitude analysis tells us that the majority of our atmosphere's energy enters the atmosphere from the sun than from the crust or the rest of the sky. Just think how quickly the temperature can drop at night when the sun goes down, and how much colder it gets on long winter nights than short summer nights. The earth's core is heated by radioactive decay, and cools by releasing energy into the atmosphere directly through the crust and indirectly through the oceans. Sometimes this release of energy takes the spectacular form of a volcanic eruption. However, the presence of permafrost, glaciers, and ice sheets at the earth's poles is a clear indicator that this is a relatively small contribution over all. The oceans act mostly as large reservoirs for heat, which they absorb, redistribute, and re-emit. The cooling of the atmosphere then occurs mainly through emission of energy out into the rest of the sky.

If we hypothesize that the earth's average temperature is roughly at steady-state, we can study the atmosphere's energy using an energy-balance model. The conservation of energy is a fundamental principle of the theory of thermodynamics. If, at steady-state, our atmosphere's energy is conserved, then the rate of change of the earth's atmosphere's energy must be equal to the power (energy per time) input (\(I\)) from the sun minus the power emission (\(E\)) by the earth back to the sky. Put as an equation,

\[I - E = 0.\]

To determine the power input and emission, we must apply some knowledge from thermodynamics and geometry. Naively, the toughest part is determining the power outputs. However, there is a very interesting and useful observation from the thermodynamics called the Stefan--Boltzmann relation. The Stefan--Boltzmann relation says that the power (P) emitted from an ideal body in the form of light is proportional to fourth power of the body's absolute temperature (\(\theta\)), \(P \propto \theta^4\). So, the hotter a body is, the brighter it shines. If we know the temperatures of the earth and sun, we can determine their power outputs. The rest, then, is geometry.

Labeled diagram of earth-sun configuration for our energy-balance model

Labeled diagram of earth-sun configuration for our energy-balance model

To work out the geometry, we can draw a simple diagram of the solar system consisting of the earth (represented with astronomical symbol \(\oplus\)) and the sun (represented with symbol \(\odot\)). Let \(r _ \odot\) be the radius of the sun, \(d\) be the average distance from the sun to the earth, and \(r _ \oplus\) be the radius of the earth. Through study, we have determined that \(r _ \oplus \approx 6.4 \times 10^6\) meters (using half the Karman line), \(r _ \odot \approx 7 \times 10^8\) meters, and \(d \approx 1.5 \times 10^{11}\) meters.

The total power output from the sun is the output per unit area times the area of the sun. Approximating the sun as a sphere, we know this area is 4 times the cross-sectional area. Thus, the total power output will be given by \[P _ \odot A _ \odot = ( \sigma \theta _ \odot^4)(4 \pi r _ \odot^2) \] where \(\sigma = 5.67 \times 10^{-8}\) watts per quartic degrees is the Stefan--Boltzmann constant. Scientists have studied the sun closely and found it's surface temperature \(\theta _ \odot = 5800^{\circ}\) K.

The fraction of that power that hits the earth will be the fraction of the spherical shell occluded by the earth, divided by the area of the spherical shell with a radius equal to the the earth's average distance from the sun, \(d\). For a planet that is small and far from the sun, the area occluded is closely approximated by the cross-sectional area of the earth, so the fraction occluded is \((\pi r _ \oplus^2)/(4 \pi d^2)\). Thus, the earth's power input

\[I = P _ \odot A _ \odot \frac{\pi r _ \oplus^2}{4 \pi d^2}.\]

Even though we don't usually think of the earth as an incandescent object emitting energy, the exact same reasoning we used for the sun's power output applies for the earth. \[P _ \oplus = \sigma \theta _ \oplus^4\] So, the earth should emit energy at a rate \[E = P _ \oplus A _ \oplus = (\sigma \theta _ \oplus^4) (4 \pi r _ \oplus^2).\] Applying conservation of energy, the power absorbed by the earth must balance the power emitted back into space by the earth: \(I = E\). So, by substitution, \[\begin{gather} P _ \odot A _ \odot \left(\frac{\pi r _ \oplus^2}{4 \pi d^2}\right) = P _ \oplus A _ \oplus \\ \left( \sigma \theta _ \odot^4 \right) \left( 4 \pi r _ \odot^2 \right) \left( \frac{\pi r _ \oplus^2}{4 \pi d^2} \right) = \left( \sigma \theta _ \oplus^4 \right) \left(4 \pi r _ \oplus^2 \right) \end{gather}\] Solving for the earth's temperature, we get the simpler formula \[\theta _ \oplus = \theta _ \odot \sqrt[2]{\frac{r _ \odot}{2 d}} = 280^{\circ} K = 44^{\circ} F = 7^{\circ} C\] So, this geometric model predicts that the earth's temperature is about 44 degrees Fahrenheit. We should really caveate this statement, as the prediction is a long-term prediction for the whole earth, rather than the temperature at any one place and time. Our model predicts that the earth's global average annual temperature is about 44 degrees. Using satellite technology and other resources, climatologists estimated the earth's average temperature in 2013 was \(58^{\circ} F = 15^{\circ} C\). So our estimate is definitely in the right ballpark. If fact it seems very good!

It should be noted, however, that we've left out one particularly important component in our calculation -- albedo. The albedo of a body is its tendency to reflect incident light at the same wavelengths it arrives, rather than to absorb and re-emit that energy as black-body radiation. This is important because the reflected light typically has higher frequency than the emitted light, and thus has a different impact on temperature. From satellite measurements, we know the earth's albedo \(\alpha \approx 0.3\). When we incorporate albedo into our energy balance equation, \[ \left( \sigma \theta _ \odot^4 \right) \left( 4 \pi r _ \odot^2 \right) \frac{\pi r _ \oplus^2}{4 \pi d^2} (1-\alpha) = 4 \pi r _ \oplus^2 \sigma \theta _ \oplus^4\] Solving for the global average annual temperature again, we find \[\theta _ \oplus = \theta _ \odot \sqrt[4]{\frac{ (1-\alpha) r _ \odot^2}{4 d^2}} = 256^{\circ} K = 1^{\circ} F\]

It was a happy accident that our first calculation came so close to the earth's actual average temperature. This formula that includes albedo is a more accurate formula, and predicts an earth that is significantly colder than we experience. Part of the explanation for the difference between observed and predicted temperatures is the earth's atmospheric chemistry and dynamics, which reduce the rate of energy emissions back into space.

Sellers--Budyko theory

You may have wondered why we're modelling the global average annual temperature. You may have even asked yourself what use is the "global average annual temperature" anyway? As a person at one place and time, we never experience this average temperature. We experience the local temperature with it's regular variation from latitude, season, time of day, and weather. And you would be right. Global average annual temperature is a useful summary statistic for tracking certain trends, but a better model will capture some of the variation in temperature we know is important.

One important component of climate temperature variation is the variation of temperature with latitude. The differences in climate and weather between polar and tropical latitudes is dramatic, and global weather patterns are largely attributable to the atmospheric mixing driven by these temperature differences. One model of temperature variation with latitude and time is Budyko's model.

Model set-up

The basic setup of Budyko's climate model is based on the earth's energy balance -- energy comes in from the sun and is radiated back out, and over decade time scales, will be in equilibrium when input rates match the output rates. However, energy input and output rates are influenced by some important things that we need to take into account. For one thing, the rate of energy input depends on latitude -- because of the earth's spherical shape, the north and south poles get less light than the equator. Latitudinal temperature differences are further effected by snow and ice that alter the earth's albedo and reflect light back into space. The north and south poles have lots of ice, while the equator has very little. Weather works to undo the differences in heating between the poles and equator -- cold air gets moved south, while warm air gets moved north. We can summarize this energy-balance dynamic as \[\frac{\partial \theta}{\partial t} \propto I - E + D\] where \(\theta(y,t)\) is the typical temperature at latitude \(y\) at time \(t\), \(I\) is the rate of energy input from the sun, \(E\) is the rate of energy emission from the earth back into space, and \(D\) is the rate of atmospheric mixing evening out the temperature. We will use \(y\) as the rectilinear latitude (scaled linearly from 0 at the equator to 1 at the north pole), and \(t\) is time (in years).

The rate of energy input \(I\) from the sun can be treated geometrically, as we did with our energy-balance model, with additional terms that account latitude dependences on geometry and albedo. The average annual input rate per unit of latitude \[\begin{gather*} I = Q s(y) [1- \alpha(y, y _ c)] \end{gather*}\] where \(Q\) is the per unit area, \(s(y)\) is the effective area at latitude \(y\), and \(\alpha(y, y_c)\) is the effective albedo at latitude \(y\) when the iceline is at latitude \(y_c\). Using integration with observations of the earth's tilt and orbit, the geometry term can be reasonably approximated by the parabola \(s(y) = 1 - 0.241 (3 y^2 - 1)\). The albedo is relatively constant except for the jump it takes at poles because of ice, so it may be approximated by the piecewise constant function \(\alpha(y, y _ c) = 0.32 + 0.3 H(y - y _ c)\), where \(H(x)\) is the Heaviside function (which we leave undefined at zero for the moment).

The rate of energy emission \(E\) could be approximated with the Stefan--Boltzmann relation, but for the range of temperatures we are interested in, it is simpler and just as accurate to use a linearized version of this relation, so \[\begin{gather*} E = A + B \theta. \end{gather*}\] Now, we must include something for the term \(D\) that accounts for movement of heat between latitudes. Since mixing just moves energy around, we still expect the total energy to be conserved, so \[\int _ {-1}^{1} D dy = 0.\] The idea Budyko used was that global mixing moves local temperatures towards the global average temperature according to Newton's law of cooling: \[\begin{gather*} D = C (\overline{\theta} - \theta ) \end{gather*}\] where \(C\) is a constant and \(\overline{\theta}\) is the global average temperature.

Putting all of these pieces together, \[\begin{gather} \frac{\partial \theta}{\partial t} \propto Q s(y) (1-\alpha(y,y_c)) - (A+B\theta) + C \left( \overline{\theta} - \theta \right), \\ R \frac{\partial \overline{\theta}}{\partial t} = Q (1-\overline{\alpha}(y_c)) - (A+B\overline{\theta}), \\ \overline{\theta}(t) = \frac{1}{2} \int _ {-1}^1 \theta(y,t) dy = \int_0^1 \theta(y,t) dy, \\ \overline{\alpha}(y_c) = \int _ {-1} ^{1} s(y) \alpha(y,y _ c) dy, \end{gather}\]

where \(Q = 340\) watts per square meter, \(A + B \theta = 202 + 1.9 \theta\), \(C = 1.6\), \(\alpha(y, y _ c) = 0.32 + 0.3 H(y - y _ c)\), \(H(x)\) is the Heaviside function, and \(s(y) = 1 - 0.241 (3 y^2 - 1)\). This is the dynamic Budyko model for global temperature dynamics. The Budyko model can be rewritten as the single partial differential equation \[R \frac{\partial \theta}{\partial t} = Q s(y) \left(1 - \alpha{\left (y,y _ c \right )} \right) - (A + B \theta) + C \left( \int_{0}^{1} \theta{\left (y,t \right )}\, dy - \theta{\left (y,t \right )} \right)\]

Model analysis

Budyko's theory is puzzling at first. It's equations look and feel very different from the models we have previously studied. At first, it seems like a partial differential equation for \(\theta(y,t)\), but there are no partial derivatives with respect to \(y\), only a definite integral and some functional dependences. There is also a discontinuity in our albedo function at the ice line latitude \(y_c\), and it's not entirely clear yet how we should handle this ice line.

Budyko's theory is an example of a Stefan problem -- a problem where an unknown internal moving boundary plays an important role in the solution. Stefan problems were originally studied in the context of freezing lakes, where we were interested using temperature observations to estimate ice thickness. Today, generalized versions of Stefan problems called free boundary problems are of wide interest. In our case, the importance of the ice line is not evident a priori, but will become clear momentarily.

To make sense of Budyko's model, let us start by hoping the dynamics will approach a steady-state solution after sufficient time, and attempt to find this steady-state.

At steady-state, \(\partial \theta/\partial t = 0\), so the steady-state temperature \(\theta^*\) solves \[0 = Q s(y) \left(1 - \alpha{\left (y,y _ c \right )} \right) - (A + B \theta^*) + C \left( \overline{\theta}^* - \theta^* \right)\] Since the local temperature \(\theta\) appears linearly on the right hand, we can rearrange to show \[\theta^*(y,y_c) := \frac{ Q s(y) \left(1 - \alpha{\left (y,y _ c \right )} \right) - A + C \overline{\theta}^* }{B+C} \] where the steady-state global annual average temperature \[\overline{\theta}^ * (y _ c) := \frac{Q}{B} \left( 1 - \overline{\alpha}(y _ c)\right) - \frac{A}{B}.\] For the moment, let's treat the ice-line latitude as an independent variable, even though its location should depend on the temperature. When we know where the ice line latitude \(y _ c\) is, we can calculate the steady-state average temperature \(\overline{\theta}^ * (y _ c)\).

[Show code]
#!/usr/bin/env python3
We are working with the dynamic Budyko model and governing law ...

>>> Eq(R*T(y,t).diff(t), Q*s(y)*(1-alpha(y,y_s)) - (A+B*T) + C*(Integral(T(y,t),(y,0,1)) - T(y,t)))

                                                  / 1                     \
      dT                                          | /\                    |
    R -- = Q (1 - a(y, y_s)) s(y) - (A + B T) + C | | T(y, t) dy - T(y, t)|
      dt                                          |\/                     |
                                                  \ 0                     /

from pylab import *

def heaviside(x):
    return (0. if x < 0. else ( 1. if x > 0. else 0.5))

def zap(x):
    # 1 if x >= 0 otherwise NaN
    return ( float('nan') if x < 0. else 1. )

Q_real = 343.
#Q = 330.
A = 202.
B = 1.90
C = 1.6*B
TIce = -10. # average temperature at ice line

def s(y):
    # Latitudinal distribution approximation
    return  1 - 0.482*(3*y**2-1)*0.5

def alpha(y, ys):
    # albedo as a function of latitude and ice-line latitude
    return 0.32 + heaviside(y-ys)*0.3

def mid_alpha(ys):
    # this works only because alpha is a flat Heaviside
    return (alpha(0., ys) + alpha(1., ys))*0.5

def net_alpha(ys):
    # average global albedo, which is Integral(s(y)*alpha(y,ys),(y,0,1))
    return 0.62 - 0.3*ys*(1-0.241*(ys**2-1))

def average_T(ys, Q):
    return (Q*(1 - net_alpha(ys)) - A )/B

def T(y, ys, Q):
    # local steady-state temperature
    return (Q*s(y)*(1 - alpha(y, ys)) - A + C*average_T(ys, Q))/(B + C)

def midpoint_temperature(ys, Q):
    # at the midpoint of the ice line
    return ( Q*s(ys)*(1 - mid_alpha(ys) ) -A + C*average_T(ys, Q))/(B+C)

def critPower(ys):
    return ((B+C)*TIce + A - C*average_T(ys, Q_real))/(s(ys)*(1-mid_alpha(ys)))

def midTemp(ys):
    return average_T(ys, critPower(ys) )

def topTemperature(Q):
    # assume an ice-free world
    ys = 1.
    y = 1.
    return zap(T(y, ys, Q)-TIce)*average_T(ys, Q)

def bottomTemperature(Q):
    # assume an ice-covered world
    ys = 0.
    y = 0.+1e-5
    return zap(TIce-T(y, ys, Q))*average_T(ys, Q)

def main():
    figure(1, figsize=(12, 10))
    y = linspace(1e-4, 1-1e-4, 1000)
    print(T(1e-3, 0., Q_real), T(0.9, 0., Q_real))
    print(T(0., 1., Q_real), T(1., 1., Q_real))

    subplot(2, 2, 4)
    for ys in linspace(0, 1, 10):
        plot(y, T(y, ys, Q_real), 'k-')
    #plot(y, -10+0*y, 'r:', y, midpoint_temperature(y, Q_real), 'g--')
    plot(y, -10+0*y, 'r:')
    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Local temperature for each $y_s$', fontsize=16)
    text(0.1, -50, '$T(y, y_s)$ for different ice lines $y_s$', fontsize=14)
    ylim(-60, 30)

    subplot(2, 2, 3)
    plot(y, s(y), 'b-')

    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Latitudinal distribution $s(y)$', fontsize=16)

    subplot(2, 2, 1)
    plot(y, alpha(y, 0.5), 'g-')
    xlabel('Latitude (normalized)', fontsize=16)
    ylabel('Albedo if ice line $y_s= 0.5$', fontsize=16)

    subplot(2, 2, 2)
    plot(y, net_alpha(y), 'g-')
    xlabel('Latitude of ice line $y_s$ (normalized)', fontsize=16)
    ylabel(r'Global average albedo $\overline{\alpha}(y_s)$', fontsize=16)

    savefig('snowball.png', transparent=True)

    for ys in linspace(0, 1, 10):
        plot(y, T(y, ys, Q_real), 'k-')
    plot(y, -10+0*y, 'r:', y, midpoint_temperature(y, Q_real), 'g--')
    xlabel('Latitude (normalized)', fontsize=18)
    ylabel('Local temperature for each $y_s$', fontsize=18)
    ylim(-60, 30)


    ys = linspace(1e-8, 1-1e-8, 100)
    Q = array([ critPower(i) for i in ys ])
    plot(Q, ys, 'b-', linewidth=3)
    plot([343, 343], [0, 1], 'g--', linewidth=2)
    QL = linspace(-250,550,100)
    plot(QL, zap(T(1., 1., QL)-TIce)*1, 'k-', linewidth=3)
    plot(QL, zap(TIce - T(1e-5, 0., QL))*0, 'k-', linewidth=3)
    xlim(270, 470)
    xlabel('Solar power (Q watts per meter square)', fontsize=18)
    ylabel('Stationary ice line lattitude ($y_s$)', fontsize=18)

    mT = array([ midTemp(i) for i in ys])
    plot(Q, mT, 'b-', linewidth=1)
    Q = linspace(270, 470, 1000)
    #hT = array([ highTemp(i) for i in Q])
    #lT = array([ lowTemp(i) for i in Q])
    #plot(Q, lT, 'k-', Q, hT, 'k-')
    plot(Q, topTemperature(Q), 'k-', Q, bottomTemperature(Q), 'k-')
    #plot(Q, T(1e-3, 0., Q), 'g-', Q, T(1.-1e-3, 0., Q), 'g-')
    #plot(Q, T(1e-3, 1., Q), 'r-', Q, T(1.-1e-3, 1., Q), 'r-')
    #plot([Q_real]*2, [-80, 100], 'b:')
    xlabel('Solar input rate (Q W/m$^2$)', fontsize=18)
    ylabel('Global temperature ($^\circ$C)', fontsize=18)
    ylim(-50, 60)
    xlim(Q[0], Q[-1])

    ys = linspace(0, 1-1e-5, 100)
    aT = average_T(ys, Q_real)
    plot(ys, aT, 'k--')
    xlabel('Ice line lattitude $y_s$(normalized)', fontsize=18)
    ylabel('Global average temperature (degrees C)', fontsize=18)


Global annual average temperature treated as a function of the ice-line latitude y_c.

Global annual average temperature treated as a function of the ice-line latitude \(y_c\).

Combining our descriptions of albedo \(\alpha(y,y _ c)\) and solar geometry \(s(y)\) with the average temperature \(\overline{\theta}^ * (y _ c)\), we can calculate the steady-state latitudinal temperature profile for different ice lines.

Further example climate model calculations

Further example climate model calculations

Note that our steady-state latitudinal temperature profiles never overlap -- each point between the maximum and minimum solutions falls on exactly one steady-state solution. This is not aways the case, and a special consequence of piecewise continuity and monotonicity of albedo. We also see that the annual average temperature decreases monotonely with latitude.

The catch with our analysis so far is that the ice line's location is not an independent variable. The location of the ice line should be a function of the earth's temperature profile. The ice line will be close to latitudes where the annual average temperature is \(0\) degrees Celsius. It won't be at exactly 0, since seasonal variation in temperature will melt ice during the summer at that latitude. Let's call the critical annual average temperature below which ice persists year round \(\theta _ c\). Evidence suggests the critical annual average temperature \(\theta _ c = -10^{\circ}\) C.

Where is the average temperature equal to the critical temperature? There are two extreme cases --- a Dune-like planet without ice (\(y _ c = 1\)) and a Hoth-like planet encased in ice (\(y _ c = 0\)). In the Dune case, we see \(\theta^*(y,1) > \theta _ c\) everywhere, suggesting the planet is plenty warm enough to stay ice-free. But this will only hold as long as \(\theta^*(1,1) > \theta _ c\). On the other hand, in the Hoth case, we see \(\theta^*(y,0) < \theta _ c\) everywhere, suggesting the planet would be cold enough to stay ice-covered everywhere as long as \(\theta^*(0,0) < \theta _ c\). For our estimated parameter values, both the Dune and Hoth scenarios are valid steady-state solutions.

Now, what about intermediate ice-line latitudes? Are there latitudes where our model predicts the temperature at the ice line will be the critical temperature \(\theta_c\)? Well, from the plots above, we see the latitudes where the iceline reaches critical temperaure are at the ice line?!?! Our model predicts a big jump in local average temperature whereever the ice line occurs, and for the parameter values we are considering, the local temperature jumps actually straddle the critical temperature!

Without some more information, we can not say precisely where the temperature profile will cross the critical temperature. We need some other condition. Perhaps what we need is to define the annual average temperature at the ice line? A natural (but not necessarily optimal) idea is to define the temperature at the ice-line latitude \(y _ c\) as the average of the temperatures at the immediately higher and lower latitudes: \[\begin{gather} \theta^*(y_c,y_c) = 1/2 \; \lim_{\epsilon \rightarrow 0} \left[ \theta^*(y_c + \epsilon, y_c) + \theta^*(y_c - \epsilon, y_c) \right]. \end{gather}\]

We call this extra equation a "jump condition" or "Stefan condition". Stefan conditions are a simple case of the more general free-boundary boundary conditions that can appear in the study of metal melting, bubble movement, and blood circulation. At steady-state, we expect the Stefan condition to give the ice-line's critical latitude \(y_c^*\) implicitly by \(\theta _ c = \theta^*(y_c^*, y_c^*)\). We can plot the average temperature at the iceline as a function of the iceline latitude to determine solutions.

Steady-state ice line temperature (green) as a function of the ice-line location y_c.

Steady-state ice line temperature (green) as a function of the ice-line location \(y_c\).

What we discover is that the ice-line's temperature crosses the critical threshold of \(-10^{\circ}\) Celsius at two separate points. These two crossing points each correspond to different steady-state solutions \(\theta^ * (y, y _ c^ * )\). In between these two points, our jump condition predicts the temperature will not be cold enough for ice year round. Let's call the larger of these two steady-state solutions Earth (since it seems closest to our current circumstances), and let us call the smaller solution Gethen, in honor of Ursula K. Leguin's speculative fiction world.

If we plot the iceline latitudes of the 4 steady-state solutions (Dune, Earth, Gethen, and Hoth) as a function of the incoming solar power, we get a more complete picture of the situation. For our realistic parameter estimates, all four steady-states exist. The Hoth and Dune steady-states only exist for temperatures where there temperature profiles don't intersect the critical temperature. If there is much less incoming power, then only the Hoth steady-state exists. If there was much more power coming in, then only the Dune steady-state would exist.

Steady-state ice-line latitudes as a function of incoming solar power

Steady-state ice-line latitudes as a function of incoming solar power

Steady-state global annual average temperature depending on incoming solar power Q. Depending on the power input, there may be 1, 2, 3, or 4 steady-state solutions.

Steady-state global annual average temperature depending on incoming solar power Q. Depending on the power input, there may be 1, 2, 3, or 4 steady-state solutions.

Analyzing stability

So far, our model gives us multiple predictions about the steady-state of earth's climate. But just because it is a steady-state doesn't mean it is realistic. A steady-state also has to be stable, in some sense, in order for it to correspond to something we are likely to observe in real life. Unfortunately, we'll have to develop a new method to assess stability of a system with a discontinuity.

Recall from our elementary differential equations studies that for a one-dimensional autonomous ordinary differential equation \(\dot{u} = f(u)\), any steady-state solution \(u^*\) solving \(f(u^*)=0\) is locally stable when \(f'(u^ * ) < 0\). Let's try to apply this same reasoning to analyze the stability of our 4 steady-states using the evolution equation for global annual average temperature. Dynamically, \[\begin{gather} \frac{d \overline{\theta}}{dt} = Q \left(1 - \overline{\alpha}(y _ c(\overline{\theta})) \right) - ( A + B \overline{\theta} ). \end{gather}\] Here, we've used the implicit function theorem to treat the ice-line location as a function of the average temperature. Differentiating the right hand side with respect to \(\overline{\theta}\), we find the average temperature will be stable provided \[ - Q \frac{\partial \overline{\alpha}}{\partial y_c} \frac{\partial y_c}{\partial \overline{\theta}} - B < 0.\] But the inequality is inconclusive because of the partial derivatives we have not yet defined. But, now, let's do something tricky -- let's differentiation our steady-state condition \[0 = Q (1 - \overline{\alpha}(y _ c) ) - ( A + B \overline{\theta} )\] with respect to the power \(Q\). \[\begin{gather} 0 = (1 - \overline{\alpha}(y _ c) ) - Q \frac{\partial \overline{\alpha}}{\partial y_c} \frac{\partial y_c}{\partial\overline{\theta}} \frac{\partial \overline{\theta}}{\partial Q} - B \frac{\partial \overline{\theta}}{\partial Q}. \end{gather}\] After substitution, our stability condition becomes \[\begin{gather} -(1-\overline{\alpha}) \left( \frac{\partial \overline{\theta}}{\partial Q} \right)^{-1} < 0 \end{gather}\] which we can simplify to \[\begin{gather} \frac{\partial \overline{\theta}}{\partial Q} > 0 \end{gather}\]

for interior steady-states. This is a classic result known as the slope-stability theorem by (Cahalan and North, 1979) with an elementary interpretation -- any equilibrium global average annual temperature with an interior ice-line will be a stable equilibrium if the global temperature increases as power increases. Looking at our plot, this tells us that Earth's situation is stable, while Gethen's situation is unstable. These results in-turn imply stability results for the extreme scenarios of Hoth and Dune. Hoth will be stable for the range of low power inputs and inputs where Gethen exists, while it will be unstable for larger power inputs. Dune will be unstable for low power inputs and when Earth exists, while Dune will be stable for higher power inputs. So our model says that there are two stable steady-state climates we could have -- the "Earth-like" climate we have today, or a snow-ball earth "Hoth-like" climate that some geologists have speculated about.

You shouldn't accept this stability analysis un-critically. It seems to suggest that the stability of our temperature profiles can be determined entire in terms of a single ordinary differential equation. That's probably not so. If our stability condition above fails, the steady-state will be unstable, since an appropriately picked perturbation is sure to grow. However, the inverse does not necessarily hold. There may be perturbations to our steady-state solutions that are not encompassed in our simple analysis and lead to instability. However, the presence of a discontinuity in our steady-state solution breaks the conventional methods we would usually try to use. We will have to leave the questions concerning the local stability of the steady-states in Sellers--Budyko theory incompletely answered.


Fortunately, our planet is not like Hoth right-now. But the geological record does indicate we've had ice-ages and maybe even snowball earths in the past. That means there are probably some slow but powerful effects that are moving our climate back and forth allong this spectrum. Pushed in one direction, our climate will steadily warm up until there are no more ice caps. This change would be reversible. Pushed the other way, climate will cool for a bit. But if it was cooled enough, it would hit a turning point of no return, where the spread of the ice sheet would accelerate and encase the whole planet in glaciers. That change could be irreversible on the time scale of human lives. Eventually, the transition back to our familiar climate might be equally rapid. This process of long periods of slow change alternating with short periods of rapid change is called a relaxation oscillation. Relaxation oscillations can be very dangerous because there way be very little warning or opportunity to prepare before the onset of rapid change.

Budyko's climate theory is provoking, in that, beginning with the familiar theories of temperature, conservation of energy, and the simply-stated Stefan--Boltzmann law, we can derive climate predictions for entire planets that are surprisingly good. It's even more provocative because there is geological evidence supporting the appearance of ice ages, and maybe even snowball earths. But we should view Budyko's model with a fair degree of skepticism as well, given our relative naivety about the complexities of climate science. The model has multiple equilibria because of the particular Stefan condition we've adopted. There are two interior equilibria because the Stefan condition for the ice-line creates a concave function with an interior peak. If we examine that condition closely, it begins to appear somewhat arbitrary, and we might wonder if there are other equally reasonable alternatives that will flip that concavity. Over the course of a year, we expect the ice line to have strong seasonal oscillations which would demand a more complicated stability analysis, perhaps based on Floque theory. At a higher level, our definition of global annual average temperature is suspect because it seems to weigh all latitudes equally. And our submodel for heat redistribution by atmospheric mixing seems little more than a phenomenological convenience, about as realistic and complete a representation of Rembrant the person as one of his self-portraits.

Overcoming these criticisms and building really good predictive models of our planet's global climate is an incredibly challenging research problem. There are many details that have to be including in the accounting, like mountain ranges, ocean circulation, ice sheet melting, cloud formation, daily weather patterns, and even salinety fresh-water inflows from large rivers. Even our largest, most expensive computers are barely able to manage the demanding task currently, and a great deal of research remains for those interested and up for the challenge.

there is not just one tipping point, but many -- melting of glaciers, collapse of ocean circulation, decaying of permifrost, ... We don't know where these tipping points are or how/if we can avoid them. Steffen, 2018


  1. What fraction of the sun's power output does the earth re-emit back to the sun?

  2. In our calculations above, we ignored power entering the atmosphere from the light bouncing off the moon. Estimate the average magnitude of power from moonlight relative to direct solar power.

  3. In our energy-balance model, we used an approximate formula for the fraction of power output from the sun occluded by the earth. Find the exact formula for the fraction of power output from a point source occluded by a sphere of radius \(r\) a distance \(d\) from the source (distance being measured to the nearest point on the sphere. Show that our approximate formula is the correct asymptotic approximation as \(r/d \rightarrow 0\). How much error does this approximation introduce to our calculations?

  4. If energy-balance theory of global climate is correct, it should be applicable not just to the earth but to other bodies in our solar system as well. Lets consider the special case of Mars.

    1. What data would you need to predict the average temperature of Mars?

    2. Where can you find this data on the internet? (Please choose a source or sources that are trusted and have a verifiable citation.)

    3. Using the methods from class, predict the average temperature on Mars.

    4. In 1976, the Viking 1 lander became the first spacecraft to reach mars. As part of its mission, Viking 1 recorder the temperature on Mars ( Plot the actual temperature over 10 sols (a Martian day), and discuss the relationship between the predicted and observed values.

    5. (Hard) Extend our climate model to predict daily temperature variation on Mars, and compare your results with Viking 1's data.

  5. If we want to build a ringworld in our solar system that has an average climate of 295 degrees Kelvin and an albedo 30 %, what should the ringworld’s radius be (in astronomical units)? You can assume the outer sides of the ringworld are perfectly insulated and do not lose heat.

  6. to finish and test After a sudden weather change, a pond that was uniformly 35 degrees Fahrenheit is subject to a cold snap of 15 degrees Farenheit for a week. How thick is the ice on the top of the pond at the end of the week?

  7. Estimate the power entering the atmosphere from the earth's crust.

( previous, home, next )