Climate and snowball earth (3 days)

( previous, home, next )


frozen planet earth
frozen planet earth

One of the great contemporary scientific efforts of the day is our study of the climate of planet Earth. Climate studies are somewhat unique, in that climate effects all of our lives daily. Our world today does 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 our own actions are changing the climate.

“Climate” is an intangible thing, representing the aggregation of chaotic weather conditions that vary in time and space, and which each one of us experiences only vary narrowly. With only limited personal experience, it is hard for us to develop intuition, experience, and foresight about changing climatic conditions. Modelling is one approach that can help us develop that intuition. In this chapter, we will develop two simple models of earths climate. The first model attempts to predict the average temperature of the earth using only geometry and one well-known 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. Along the way, we’ll come up with a new kind of mathematical model and some ideas for studying it.

Seller’s theory of 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 energy can enter and leave 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 Seller’s model
labeled diagram of earth-sun configuration for Seller’s 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. Astronomers have determined that \(R _ \oplus \approx 6.4 \times 10^6\) meters, \(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, which, for a planet that is small and far from the sun, should be approximate the cross-sectional area of the earth divided by the area of the spherical shell, \[\frac{\pi R _ \oplus^2}{4 \pi D^2}\]

\[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 applies. \[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 K = 1 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.

Budyko’s global temperature model

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.

Perhaps the most important component of 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 = \sin(\theta)\) as the normalized latitude (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 Seller’s 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 will 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*} R \frac{\partial \Theta}{\partial t} = 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 are albedo function at the ice line, and it’s not entirely clear yet how we should handle this ice line.

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

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 \[0 = Q s(y) \left(1 - \alpha{\left (y,y _ c \right )} \right) - (A + B \Theta) + C \left( \overline{\Theta} - \Theta{\left (y,t \right )} \right)\] Since the local temperature \(\Theta\) appears linearly on the right hand, we can rearrange to show \[\Theta(y) = \frac{ Q s(y) \left(1 - \alpha{\left (y,y _ c \right )} \right) - A + C \overline{\Theta} }{B+C} \] where \[\overline{\Theta} = \frac{Q}{B} \left( 1 - \overline{\alpha}(y _ c)\right) - \frac{A}{B}.\]

So, if we know where the ice line latitude \(y _ c\) is, we can calculate the steady-state global annual 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 what we’ve proposed about albedo \(\alpha(y,y _ c)\) and planetary geometry \(s(y)\) with this calculation of global temperature \(\overline{\Theta}^ * (y _ c)\), we can calculate the latitudinal temperature profile for different ice line latitudes.

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 a free parameter. 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? Let’s first consider the two extreme case of a Dune-line 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? Well, from the plots above, we see that it’s at the ice line?!?! Our model predicts a big jump in local average temperature at the ice line, and for the parameter values we’ve picked, the local temperature jumps from above the critical temperature for lower latitudes to below the critical temperature at higher latitudes.

Without some more information, we can not say precisely where the temperature profile will cross the critical temperature threshold. 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,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” after the first scientist to make use of one. Stefan conditions are a simple case of the more general free-boundary boundary conditions that can appear in the study of lake freezing, bubble movement, and blood circulation.

From our jump condition, we can plot the annual average temperature at the iceline as a function of the iceline latitude. At steady-state, we expect the Stefan condition to return the ice-line’s critical temperature: \[\begin{gather*} \Theta _ c = 1/2 \; \lim_{\epsilon \rightarrow 0} \left[ \Theta(y_c + \epsilon, y_c) + \Theta(y_c - \epsilon, y_c) \right]. \end{gather*}\]

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 degrees Celsius at two separate points. These two crossing points each correspond to a different steady-state solutions \((\overline{\Theta}^ * , \Theta^ * (y), y _ c^ * )\) where are two steady-state equations and the steady-state Stefan condition. In between these two points, our jump condition predicts the temperature will not be cold enough for ice year round. Let’s call the upper of these two crossing points Earth (since it seems closest to our current circumstances), and let us call the lower point 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. 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.

Stability of steady-states using perturbation analysis

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\). We can apply this same reasoning to analyze the stability of our 4 steady-states using the evolution equation for global annual average temperature

\[\begin{gather} \frac{\partial \Theta}{\partial t} = Q (1 - \overline{\alpha})(y _ c) ) - ( A + B \overline{\Theta} ).\end{gather}\]

Differentiating with respect to \(\overline{\Theta}\), we find the global temperature will be stable provided \[ - Q \frac{d \overline{\alpha}}{d y_c} \frac{d y_c}{d\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{d \overline{\alpha}}{d y_c} \frac{d y_c}{d\overline{\Theta}} \frac{d \overline{\Theta}}{d Q} - B \frac{d \overline{\Theta}}{d Q} \end{gather}\]

After substitution, our stability condition becomes \[\begin{gather} -(1-\overline{\alpha}) \left( \frac{d\overline{\Theta}}{dQ} \right)^{-1} < 0 \end{gather}\] which we can simplify to \[\begin{gather} \frac{d\overline{\Theta}}{dQ} > 0 \end{gather}\] for interior steady-states. This is a classic result 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. When we apply this condition, we see that Earth’s situation is stable, while Gethen’s situation is unstable!

Long-time analysis


Budyko’s climate theory is compelling, in that, beginning with the familiar idea of conservation of energy and the simply-state Stefan–Boltzmann law, we can derive climate predictions for entire planets. It’s even more compelling because there is geological evidence supporting the appearance of ice ages, and maybe even snowball earths.
supporting references?
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. See Steffen, 2018

But we have to 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.

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, including mountain ranges, ocean circulation, ice sheet melting, cloud formation and dispersal, daily weather patterns, and even 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.


  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. If Seller’s 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.

  4. 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.

  5. 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?

  6. Estimate the power entering the atmosphere from the earth’s crust.

( previous, home, next )