Mathematics of the Environment (Part 5)

We saw last time that the Earth’s temperature seems to have been getting colder but also more erratic for the last 30 million years or so. Here’s the last 5 million again:

People think these glacial cycles are due to variations in the Earth’s orbit, but as we’ll see later, those cause quite small changes in ‘insolation’—roughly, the amount of sunshine hitting the Earth (as a function of time and location). So, M. I. Budyko, an expert on the glacial cycles, wrote:

Thus, the present thermal regime and glaciations of the Earth prove to be characterized by high instability. Comparatively small changes of radiation—only by 1.0-1.5%—are sufficient for the development of ice cover on the land and oceans that reaches temperate latitudes.

How can small changes in the amount of sunlight hitting the Earth, or other parameters, create big changes in the Earth’s temperature? The obvious answer is positive feedback: some sort of amplifying effect.

But what could it be? Do we know feedback mechanisms that can amplify small changes in temperature? Yes. Here are a few obvious ones:

Water vapor feedback. When it gets warmer, more water evaporates, and the air becomes more humid. But water vapor is a greenhouse gas, which causes additional warming. Conversely, when the Earth cools down, the air becomes drier, so the greenhouse effect becomes weaker, which tends to cool things down.

Ice albedo feedback. Snow and ice reflect more light than liquid oceans or soil. When the Earth warms up, snow and ice melt, so the Earth becomes darker, absorbs more light, and tends to get get even warmer. Conversely, when the Earth cools down, more snow and ice form, so the Earth becomes lighter, absorbs less light, and tends to get even cooler.

Carbon dioxide solubility feedback. Cold water can hold more carbon dioxide than warm water: that’s why opening a warm can of soda can be so explosive. So, when the Earth’s oceans warm up, they release carbon dioxide. But carbon dioxide is a greenhouse gas, which causes additional warming. Conversely, when the oceaans cool down, they absorb more carbon dioxide, so the greenhouse effect becomes weaker, which tends to cool things down.

Of course, there are also negative feedbacks: otherwise the climate would be utterly unstable! There are also complicated feedbacks whose overall effect is harder to evaluate:

Planck feedback. A hotter world radiates more heat, which cools it down. This is the big negative feedback that keeps all the positive feedbacks from making the Earth insanely hot or insanely cold.

Cloud feedback. A warmer Earth has more clouds, which reflect more light but also increase the greenhouse effect.

Lapse rate feedback. An increased greenhouse effect changes the vertical temperature profile of the atmosphere, which has effects of its own—but this works differently near the poles and near the equator.

See also “week302″ of This Week’s Finds, where Nathan Urban tells us more about feedbacks and how big they’re likely to be.

Understanding all these feedbacks, and which ones are important for the glacial cycles we see, is a complicated business. Instead of diving straight into this, let’s try something much simpler. Let’s just think about how the ice albedo effect could, in theory, make the Earth bistable.

To do this, let’s look at the very simplest model in this great not-yet-published book:

• Gerald R. North, Simple Models of Global Climate.

This is a zero-dimensional energy balance model, meaning that it only involves the average temperature of the earth, the average solar radiation coming in, and the average infrared radiation going out.

The average temperature will be T, measured in Celsius. We’ll assume the Earth radiates power square meter equal to

\displaystyle{ A + B T }

where A = 218 watts/meter2 and B = 1.90 watts/meter2 per degree Celsius. This is a linear approximation taken from satellite data on our Earth. In reality, the power emitted grows faster than linearly with temperature.

We’ll assume the Earth absorbs solar energy power per square meter equal to

Q c(T)


Q is the average insolation: that is, the amount of solar power per square meter hitting the top of the Earth’s atmosphere, averaged over location and time of year. In reality Q is about 341.5 watts/meter2. This is one quarter of the solar constant, meaning the solar power per square meter that would hit a panel hovering in space above the Earth’s atmosphere and facing directly at the Sun. (Why a quarter? We’ve seen why: it’s because the area of a sphere is 4 \pi r^2 while the area of a circle is just \pi r^2.)

c(T) is the coalbedo: the fraction of solar power that gets absorbed. The coalbedo depends on the temperature; we’ll have to say how.

Given all this, we get

\displaystyle{ C \frac{d T}{d t} = - A - B T + Q c(T(t)) }

where C is Earth’s heat capacity in joules per degree per square meter. Of course this is a funny thing, because heat energy is stored not only at the surface but also in the air and/or water, and the details vary a lot depending on where we are. But if we consider a uniform planet with dry air and no ocean, North says we may roughly take C equal to about half the heat capacity at constant pressure of the column of dry air over a square meter, namely 5 million joules per degree Celsius.

The easiest thing to do is find equilibrium solutions, meaning solutions where \frac{d T}{d t} = 0, so that

A + B T = Q c(T)

Now C doesn’t matter anymore! We’d like to solve for T as a function of the insolation Q, but it’s easier to solve for Q as a function of T:

\displaystyle{ Q = \frac{ A + B T } {c(T)} }

To go further, we need to guess some formula for the coalbedo c(T). The coalbedo, remember, is the fraction of sunlight that gets absorbed when it hits the Earth. It’s 1 minus the albedo, which is the fraction that gets reflected. Here’s a little chart of albedos:

If you get mixed up between albedo and coalbedo, just remember: coal has a high coalbedo.

Since we’re trying to keep things very simple right not, not model nature in all its glorious complexity, let’s just say the average albedo of the Earth is 0.65 when it’s very cold and there’s lots of snow. So, let

c_i = 1  - 0.65 =  0.35

be the ‘icy’ coalbedo, good for very low temperatures. Similarly, let’s say the average albedo drops to 0.3 when its very hot and the Earth is darker. So, let

c_f = 1 - 0.3 = 0.7

be the ‘ice-free’ coalbedo, good for high temperatures when the Earth is darker.

Then, we need a function of temperature that interpolates between c_i and c_f. Let’s try this:

c(T) = c_i + \frac{1}{2} (c_f-c_i) (1 + \tanh(\gamma T))

If you’re not a fan of the hyperbolic tangent function \tanh, this may seem scary. But don’t be intimidated!

The function \frac{1}{2}(1 + \tanh(\gamma T)) is just a function that goes smoothly from 0 at low temperatures to 1 at high temperatures. This ensures that the coalbedo is near its icy value c_i at low temperatures, and near its ice-free value c_f at high temperatures. But the fun part here is \gamma, a parameter that says how rapidly the coalbedo rises as the Earth gets warmer. Depending on this, we’ll get different effects!

The function c(T) rises fastest at T = 0, since that’s where \tanh (\gamma T) has the biggest slope. We’re just lucky that in Celsius T = 0 is the melting point of ice, so this makes a bit of sense.

Now Allan Erskine‘s programming magic comes into play! I’m very fortunate that the Azimuth Project has attracted some programmers who can make nice software for me to show you. Unfortunately his software doesn’t work on this blog—yet!—so please hop over here to see it in action:

Temperature versus insolation.

You can slide a slider to adjust the parameter \gamma to various values between 0 and 1.

In the little graph at right, you can see how the coalbedo c(T) changes as a function of the temperature T. In this graph the temperature ranges from -50 °C and 50 °C; the graph depends on what value of \gamma you choose with slider.

In the big graph at left, you can see how the insolation Q required to yield a given temperature T between -50 °C and 50 °C. As we’ve seen, it’s easiest to graph Q as a function of T:

\displaystyle{ Q = \frac{ A + B T } {c_i + \frac{1}{2} (c_f-c_i) (1 + \tanh(\gamma T))} }

Solving for T here is hard, but we can just flip the graph over to see what equilibrium temperatures T are allowed for a given insolation Q between 200 and 500 watts per square mater.

The exciting thing is that when \gamma gets big enough, three different temperatures are compatible with the same amount of insolation! This means the Earth can be hot, cold or something intermediate even when the amount of sunlight hitting it is fixed. The intermediate state is unstable, it turns out—we’ll see why later. Only the hot and cold states are stable. So, we say the Earth is bistable in this simplified model.

Can you see how big \gamma needs to be for this bistability to kick in? It’s certainly there when \gamma = 0.05, since then we get a graph like this:

When the insolation is less than about 385 W/m2 there’s only a cold state. When it hits 385 W/m2, as shown by the green line, suddenly there are two possible temperatures: a cold one and a much hotter one. When the insolation is higher, as shown by the black line, there are three possible temperatures: a cold one, and unstable intermediate one, and a hot one. And when the insolation gets above 465 W/m2, as shown by the red line, there’s only a hot state!

Mathematically, this model illustrates catastrophe theory. As we slowly turn up \gamma, we get different curves showing how temperature is a function of insolation… until suddenly the curve isn’t the graph of a function anymore: it becomes infinitely steep at one point! After that, we get bistability:

\gamma = 0.00

\gamma = 0.01

\gamma = 0.02

\gamma = 0.03

\gamma = 0.04

\gamma = 0.05

This is called a cusp catastrophe, and you can visualize these curves as slices of a surface in 3d, which looks roughly like this picture:

from here:

• Wolfram Mathworld, Cusp catastrophe. (Includes Mathematica package.)

The cusp catastrophe is ‘structurally stable’, meaning that small perturbations don’t change its qualitative behavior. In other words, whenever you have a smooth graph of a function that gets steeper and steeper until it ‘overhangs’ and ceases to be the graph of a function, it looks like this cusp catastrophe. This statement is quite vague as I’ve just said it— but it’s made 100% precise in catastrophe theory.

Structural stability is a useful concept, because it focuses our attention on robust features of models: features that don’t go away if the model is slightly wrong, as it always is.

There are lots more things to say, but the most urgent question to answer is this: why is the intermediate state unstable when it exists? Why are the other two equilibria stable? We’ll talk about that next time!

16 Responses to Mathematics of the Environment (Part 5)

  1. I found it interesting that B is two orders of magnitude less than A. A quick and dirty numerical simulation shows that setting Q = A + B T + C T^2 with a C two orders of magnitude less than B yields the same qualitative behavior. This is no surprise if one has structural stability. However, if C is just one order of magnitude less than B, e.g. C = 0.1 than things seem to change. I am not quite sure whether that means something important, but maybe for real world problems structural stability is not enough. One might also have to argue that the linearization is indeed ‘good’ enough.

    One a second thought. C might well be zero for symmetry reasons and the physically ‘correct’ form should then be more like Q = A + B T + D T^3.

  2. amarashiki says:

    Dear John. I dedicated you this post, the 50th, in my young blog. I hope you will like it… My LOG 50th available here

    • John Baez says:

      That’s very nice! Thank you! That’s a very beautiful, detailed post, so I’m flattered that you dedicated it to me…

      … especially because I always I wish I had the time and patience to understand more about the Riemann zeta function. I’ve always loved the ‘primon gas’ idea… but nowadays, I’m busy trying to be a bit more practical, so I don’t have time to think about this. Someday, when I feel in the mood to write another pure math paper, I’ll try to finish and publish this:

      • John Baez and James Dolan, Zeta functions.

      This emphasizes the combinatorics behind the Hasse–Weil zeta function.

  3. amarashiki says:

    No, John…You…I am in debt to you more…Your This Week Finds was an incredible blog, like this Azimuth blog! Believe me, I am sure many people started blogs reading your stuff. Some time ago, I had envy…And finally, I couldn’t stand anymore…I decided to begin my blog!!! And about this post…It is not as detailed as I wanted, but more stuff is coming soon…My draft dashboard is going to explode right now, but I must focus to create good logs or posts…And it carries time.

    Well, … I think that the zeta stuff is practical stuff!!!!! It is like the fiber bundle approach behind the gauge theories…It is real math behind our (ultra)quantum world… In my opinion, we have to understand the zeta stuff better in order to understand QFT, the vacuum and, likely, quantum gravity! Of course, this is a crazy idea!

    Don’t be flattered, John…You make people to be flattered as long as we read your posts, LOL… If you read my log 51, I also wrote a basic zeta function”zoological park” for beginners. I would like to extend that in the future…But I am so disturbed by the RH than sometimes I am seduced by the own formulae you, mathematicians, invent…Hehehe…

    Best wishes!

  4. amarashiki says:

    By the way… Zeta functions are everywhere, even they arise in graph theory, hypergraph theory and it can also appear, I presume, in network theory. How? I am not sure, since I began to learn about zeta seriously only some time ago…And it is hard to keep you updated with the zeta theory…

  5. Renato Iturriaga says:

    Hi John

    I have a couple of questions and comments. About the values of A= 218 \frac{W}{m^2} and B=1.9 \frac{W}{m^2 K} , the value of A is the emission of a body at -23, very close to the radiate equilibrium of -18, I guess this is not a coincidence. Or is it? The derivative of \sigma T^4 at this value is 3.5, B is in this order of magnitude, but not as close one could imagine. Is there some reason for that?

    To see the stability or instability of a particular equilibrium is convenient to look at the graph of \frac {dT}{dt}. In this case -A-BT+Qc(T). In particular the zeros give the equilibrium points. For large negative values of T , the function is positive, and for large positive temperatures the function is negative so there is a zero somewhere, an equilibrium. Almost surely there is an odd number of them: The first time you cross the zero (the coldest equilibrium) you cross it downwards, and if you go up and reach zero you will cross it upwards, afterwards you will have to come down again. So we have either one or three, the only exception is if the zero is double, meaning that the derivative is also zero. The equilibrium points where the derivative is negative, there is local negative feedback effect: if your are slightly above the equilibrium you tend to go down if you are slightly down you tend to go up, so the point is stable, in the points where you have positive derivative -when you cross upwards- there is a positive feedback: if you are a little bit up you tend to go upwards and if you are a little bit down you tend to downwards, hence unstable. So if you have three points, the coldest and the hottest are stable and the middle one is unstable.

    • John Baez says:

      Renato wrote:

      I have a couple of questions and comments. About the values of A= 218 \frac{W}{m^2} and B=1.9 \frac{W}{m^2 K} , the value of A is the emission of a body at -23, very close to the radiate equilibrium of -18, I guess this is not a coincidence. Or is it? The derivative of \sigma T^4 at this value is 3.5, B is in this order of magnitude, but not as close one could imagine. Is there some reason for that?

      Good question! It would be nice to understand this in detail.

      The main reason the observed infrared power emission of the Earth doesn’t match the power emission of a blackbody is the greenhouse effect. Quite a bit of infrared doesn’t make it through the atmosphere. So, the simplest thing we can say is that we should have

      A + B T < \sigma T^4

      for temperatures in the range we see on Earth. But it would be nice to predict the (linearized) observed average power emission per square meter, A + B T, starting from a simple model.

      Could Simpson’s model in Part 3 be good enough to approximately calculate A and B? Or is it too simple?

      In its simplest form, this model amounts to saying that the fraction of infrared radiation that escapes to space is a function

      0 \le  F(\lambda) \le 1

      of the wavelength \lambda, given as follows. To a crude approximation, the atmosphere is:

      • completely opaque from 5.5 to 7 micrometers (due to water vapor),

      • partly transparent from 7 to 8.5 micrometers (interpolating between opaque and transparent),

      • completely transparent from 8.5 to 11 micrometers,

      • partly transparent from 11 to 14 micrometers (interpolating between transparent and opaque),

      • completely opaque above 14 micrometers (due to carbon dioxide and water vapor).

      Simpson’s function F is 0 where the atmosphere is ‘completely opaque’, 1 where it’s ‘completely transparent’, and it linearly interpolates between these values for wavelengths where the atmosphere is ‘partly transparent’.

      One could multiply the Planck distribution at temperature T by this function of wavelength, then integrate it over wavelengths, to estimate the power emitted per square meter as a function of T.

      This neglects many effects like ‘water vapor feedback’ (there’s more water vapor when it’s hotter), ‘lapse rate feedback’ and ‘cloud feedback’, as mentioned in this post. But it would still be worth doing!

      • Renato Iturriaga says:

        Hi John

        In the homepage of Kerry Emanuel ( there are a lot of things. In particular there are presentations of a lot of courses. In one of them (Tropical meteorology) there is a model for the equilibrium radiation. The simplest model has two layers, the atmosphere and the surface. The atmosphere is completely transparent to visible light and completely opaque to infrared. Let T_e be the temperature of equilibrium (250 K), then the equilibrium is achieved at at T_e in the atmosphere and 2^{\frac{1}{4}} T_e in the surface.

        A 3 layer model predicts T_e in the outer part and $3^{\frac{1}{4}} T_e$ in the surface. This is way to hot and the explanation is partially because it is not completely opaque and because of convection. Something else should be wrong because if we do the same thing with n layers, we would get a surface temperature of n^{\frac{1}{4}} T_e. But it is funny that the outer part is always the equilibrium temperature. Finally there is a complete calculation, that unfortunately is not clear how it is done, that gives a vertical profile of the temperature if only radiation is taken into account. The profile is very hot in the surface, again the explanation is the lack of the convection in the model.

    • John Baez says:

      Renato wrote:

      To see the stability or instability of a particular equilibrium is convenient to look at the graph of \frac {dT}{dt}.

      Thanks for explaining this! I hope all the students in my class read your explanation, which is better than mine.

  6. Last time we saw a ‘bistable’ climate model, where the temperatures compatible with a given amount of sunshine can form an S-shaped curve like this: [...]

  7. In Part 5 and Part 6, we studied a model where the Earth can be bistable. [...]

  8. gymnosperm says:

    Very nice. When one extends the timeline back to our horizons, basically the Proterozoic for semi reliable temperature, it is clear that bistability applies not only to glacial and interglacial periods within a snowball earth, but to the alternation between snowball and greenhouse climates. At the larger scale the warm state appears more stable.

  9. John Roe says:

    I set up a version of this model on InsightMaker, you can find it here. If you want to play around, note that A, B are as above; C is rescaled to measure time in days; G is twice gamma. Please, let me know if I messed something up.

    The time scale is quicker than I had expected, this is because the heat capacity C is quite low. In particular, the latent heat of fusion of a 1cm layer of ice appears to be of the same order of magnitude as C… one might try to incorporate this in some way, e.g. with an “ice quantity” variable driving c, and the “ice quantity” evolving according to temperature.

  10. [...] Part 5 – A model showing bistability of the Earth’s climate due to the ice albedo effect: statics. [...]

You can use HTML in your comments. You can also use LaTeX, like this: $latex E = m c^2 $. The word 'latex' comes right after the first dollar sign, with a space after it.

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Get every new post delivered to your Inbox.

Join 2,711 other followers