Climate and snowball earth (3 days)

( previous, home, next )


References


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 used to be buried in boats in Tibet, there is no river. Where people would paint pictures of themselves swimming in Egypt is now desert. Before that, geologic study has revealed that 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 sum 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 these lectures, 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.

labeld diagram of earth-sun configuration for Seller’s model
labeld 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, \[\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)\] 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 – albido. The albido 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 albido \(\alpha \approx 0.3\). When we incorporate albido 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 albido 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 lattitude, 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 lattitude. The differences in climate and weather between polar and tropical lattitudes is dramatic, and global weather patterns are largely attributable to the atmosphereic 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 decadal 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. 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).

Latitudinal temperature differences are further effected by snow and ice that alter the earth’s albido and reflect light back into space. The north and south poles have lots of ice, while the equator has very little.

\[\begin{gather*} I = ... \end{gather*}\] \[\begin{gather*} E = ... \end{gather*}\] Since mixing just moves energy around, we still expect the total energy to be conserved, so

\[\int D dy = 0.\]

\[\begin{gather*} D = C (\overline{\Theta} - \Theta ) \end{gather*}\]

Putting all of these pieces together, \[\begin{gather*} R \frac{\partial \Theta}{\partial t} = Q s(y) (1-\alpha(y,y_s)) - (A+B\Theta) + C \left[ \overline{\Theta} - \Theta \right], \\ R \frac{\partial \overline{\Theta}}{\partial t} = Q (1-\overline{\alpha}(y_s)) - (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_s) = \int _ {-1} ^{1} s(y) \alpha(y,y_s) dy, \end{gather*}\] where \(Q = 340\) watts per square meter, \(A + B \Theta = 202 + 1.9 \Theta\), \(C = 1.6\), \(\alpha(y, y _ s) = 0.32 + 0.3 H(y - y _ s)\), \(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 _ {s} \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 paritial derivatives with respect to \(y\), only a definite integral and some functiional dependences. There is also a discontinuity in are albido 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 _ {s} \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 since, we can rearrage to show \[\Theta(y) = \frac{ Q s(y) \left(1 - \alpha{\left (y,y _ {s} \right )} \right) - A + C \overline{\Theta} }{B+C} \] where \[\overline{\Theta} = \frac{Q}{B} \left( 1 - \overline{\alpha}(y _ s)\right) - \frac{A}{B}.\]

So, if we know where the ice line lattitude \(y _ s\) is, we can calculate the steady-state global annual average temperature \(\overline{\Theta}^ * (y _ s)\).

[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 *

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

@vectorize
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)


    figure(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--')
    xlabel('Latitude (normalized)', fontsize=18)
    ylabel('Local temperature for each $y_s$', fontsize=18)
    ylim(-60, 30)

    savefig('icelinesols.pdf')
    savefig('icelinesols.png',transparent=True)



    figure(2)
    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)
    text(425,0.95,'Dune')
    text(300,0.01,'Hoth')
    text(375,0.12,'Gethen')
    text(315,0.93,'Earth')
    savefig('power_to_lattitude.png',transparent=True)
    savefig('power_to_lattitude.pdf')


    figure(3)
    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])
    tight_layout()
    savefig('power_to_temp.png',transparent=True)
    savefig('power_to_temp.pdf')

    figure()
    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)
    savefig('average_temp.png',transparent=True)

    #show()


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

Combining what we’ve proposed about albido \(\alpha(y,y _ s)\) and planetary geometry \(s(y)\) with this calculation of global temperature \(\overline{\Theta}^ * (y _ s)\), we can calculate the latitudinal temperature distribution for different ice line lattitudes.

Further example climate model calculations
Further example climate model calculations

Note that our steady-state latitudinal temperature distributions are foliated – each point betten 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 albido. 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 distribution. The ice line will be close to lattitudes 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 _ s = 1\)) and a Hoth-like planet encased in ice (\(y _ s = 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 lattitudes to below the critical temperature at higher latitudes.

Without some more information, we cann’t say precisely where the temperature distribution 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 definte the temperature at the ice-line lattitude \(y _ s\) as the average of the temperatures at the immediately higher and lower lattitudes: \[\begin{gather*} \Theta(y,y_s) = 1/2 \; \lim_{\epsilon \rightarrow 0} \left[ \Theta(y_s + \epsilon, y_s) + \Theta(y_s - \epsilon, y_s) \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 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 lattitude. 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_s + \epsilon, y_s) + \Theta(y_s - \epsilon, y_s) \right]. \end{gather*}\]

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

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{T}^ * , T^ * (y), y _ s^ * )\) 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 if $u^* $ is a steady-state solution of a one-dimensional ordinary differential equation \(\dot{u} = f(u)\), then \(u^*\) is stable if \(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 _ s) ) - ( A + B \overline{T} ).\end{gather}\]

Differentiating with respect to \(\overline{T}\), we find the global temperature will be stable provided \[ - Q \frac{d \overline{\alpha}}{d y_s} \frac{d y_s}{d\overline{T}} - 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 _ s) ) - ( A + B \overline{T} )\] with respect to the power \(Q\). \[\begin{gather}0 = (1 - \overline{\alpha}(y _ s) ) - Q \frac{d \overline{\alpha}}{d y_s} \frac{d y_s}{d\overline{T}} \frac{d \overline{T}}{d Q} - B \frac{d \overline{T}}{d Q} \end{gather}\]

After substitution, our stability condition becomes \[\begin{gather} -(1-\overline{\alpha}) \left( \frac{d\overline{T}}{dQ} \right)^{-1} < 0 \end{gather}\] which we can simplify to \[\begin{gather} \frac{d\overline{T}}{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!

Discussion

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?)

But we have to view Budyko’s model with a fair degree of scepticism 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.

Exercises

  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 our theo1

    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. The Viking 1 lander recorder that actual temperature on Mars (http://www-k12.atmos.washington.edu/k12/mars/data/vl1/part1.html). 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 albido 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 )