Chapter 7: Nonlinear Phenomena in EBMs – Energy Balance Climate Models

Chapter 7
Nonlinear Phenomena in EBMs

Most of the energy balance models (EBMs) we have considered so far have been essentially linear. This allowed us the luxury of using the well-known methods of classical theoretical physics that have been developed over the last few centuries. Linear systems have many great advantages including the decomposition of fields into modes that in many cases are orthogonal, but even systems with non-orthogonal modes can be handled. Linear models can also be analyzed in very elegant forms when statistical noise drives the system, as we will see in Chapter 9. Linear systems described by partial differential equations in space and time often have symmetries that can be exploited as well. Large classes of these linear systems can be classified and their properties can be studied in great detail and generality. In particular, if the homogeneous problem is translationally invariant in time (e.g., no time-dependent coefficients), we can use Fourier methods. If in space, the geometry is simple such as translationally invariant along the line, or in the infinite plane, or rotationally on a spherical surface, we can make use of many established results. In many of these high-symmetry cases, our forebears have developed an encyclopedic archive of special functions that can be used. If no such symmetries are evident in the problem, we still can use numerical methods to find the modes, their shapes, time constants, and so on. For example, the type of equation most commonly encountered in this book is the damped diffusion equation with nonhomogeneous driving terms. A prototype is


Similar equations occur in electricity and magnetism, heat transfer, particle and molecular diffusion, as well as in many other fields. In this equation, is a linear-decay timescale, is a length scale, and so on. When the drivers are set to zero, the initial anomaly field decays to zero, spreading and smoothing as it decays. Typically, large spatial scales have longer decay times than small spatial scales. If the coefficients are time independent, the dynamical decay modes can be found and usually they form a basis set into which the solution field and the drivers can beexpanded. Modes found in the drivers will excite response modes in the solution field. In particular, frequency components appearing in the drivers will appear as energy (in our case, variance) in the responses at the very same frequencies.

Another class of problems encountered in the climate system involves a second time derivative on the left-hand side of the last equation. In this case, the response modes are wavelike. If , the waves are not damped; if , the waves are damped. There can also be both a second and a first time derivative (an example is waves on a string, set in an air-resistant medium), the first derivative being a damper as well as the term. Other than an occasional mention, we will not deal with wave equations in this book, although there are some excellent sources for that fascinating subject (e.g., Pedlosky, 2003) .

There are a few nonlinear problems that are amenable to analytical solutions in the EBM world. One obvious example is the form in blackbody radiation rules. This one turns out to be a mild nonlinearity not of much interest here, but it can be handled easily in the context of the potential function discussed in Chapter 2. This nonlinearity does not produce new solutions or alter the stability characteristics of the solutions we have already found. Another nonlinear issue arises when we consider the ice–albedo feedback mechanism. This particular mechanism leads to the same multiple-solution structure studied in Chapter 2, but there are a couple of new twists that we would like to elaborate on here. Another nonlinear problem of interest to EBM modelers (especially those working in vertical atmospheric profiles and other planets) is that of the runaway greenhouse which we considered in Chapter 4.

When more than one steady-state solution can exist for the same values of the parameters, we are always interested in the stability of such steady states. We will solve this problem for the one dimensional zonally symmetric planet in this chapter.

Figure 7.1 Shapes of the coalbedo as a function of sine of latitude, , in a solid line for the ice-cap edge at and in a dashed line at .

7.1 Formulation of the Nonlinear Feedback Model

Consider a north–south symmetric planet with no zonal features. We let the coalbedo depend on (cosine of the polar angle) along with a discontinuity at the ice-cap edge, (see Figure 7.1):


The smooth dependence multiplying the step function mimics a zenith angle dependence and comes from satellite data (Graves et al., 1993). The steady-state temperature field, , is governed by the heat conduction equation

where is a thermal diffusion coefficient; other symbols are as before. We require application of the boundary conditions:


The equatorial boundary condition ensures symmetry, the polar condition states that no heat flux enters the poles. The necessity for the latter is, of course, just a consequence of our using the polar coordinate system. Basically, it forces only the regular solution of the energy balance equation at the pole. The function is the mean annual distribution of sunlight at the top of the atmosphere which is given approximately by


It is normalized, so that


Legendre polynomials are the eigenfunctions of the diffusion operator,


with eigenvalues Note that only evenly indexed modes are retained because of the north–south symmetry. Since they form a complete basis set, we can use them to express the temperature field





Also define

After inserting (7.8) into (7.3), multiplying through by and integrating from 0 to 1, we find


We may now solve for and reconstruct by use of (7.8).

The problem is not yet solved, however, as we have not found the value of . This is done by enforcing the Budyko condition that the temperature at the ice-cap edge is 10 C. In other words,


By evaluating (7.13) at , we get an identity which constitutes a relation between and . We can plot as in Figure 7.2.

In the calculation that went into Figure 7.2, we included only the first two terms in the sum in the denominator of (7.15). Values of the parameters for the calculations were (from Chapter 5), W m (from Graves et al., 1993), = 2.00 W m K .

Figure 7.2 The solid curve is the sine of the latitude of the ice-cap edge, as a function of the solar irradiance . The vertical thin line indicates the present value of the total solar irradiance (), 340 W m. Note that there are three roots for this operating curve, I near the present climate, , II at , and III for an ice-covered planet at .

7.2 Stürm–Liouville Modes

The differential equations we have encountered so far are amenable to solution by expansion into space–time modes. Usually, the time modes are simply the Fourier Series amplitudes (harmonics). These amplitudes are often spatial modes whose graphs have more zero crossings the higher their index. There is a general class of these modes is called the Stürm-Liouville system. Consider the eigenvalue equation:

where is a continuous function on the interval. The interval of the domain can be taken to be symmetric: with boundary conditions:


where indicates the derivative evaluated at the argument in parentheses. The eigenvalues are bounded from below. To show this, multiply (7.16) by and integrate by parts. The first term of the LHS becomes


The first term after the equal sign vanishes because of the boundary conditions and the second term is positive definite. We may write


and , and we assume (as is the infrared radiative damping function in the EBMs, usually taken to be the positive constant coefficient ). All three integrals are positive, thus we have . If the function is negative over part of the interval, the will still be bounded from below, but not necessarily by a positive bound.1

7.2.1 Orthogonality of SL Modes

To prove orthogonality, begin with (7.16) and multiply through by . Rewrite the result with interchanged with . Now integrate both equations by parts from to . The LHS of the two are identical. This means we can subtract them to obtain


where the coefficient is in case the functions are not normalized. They can be normalized by dividing by . Once the modes are normalized, they form an orthonormal set. For the Legendre polynomials, . Note that the Legendre polynomials are not normalized with respect to integration from1 to 1, but rather their scale is set by the condition for historical reasons.

Other important special functions such as the Hermite polynomials and Bessel functions utilize a weight function, . We will encounter such a situation in Chapter 9 where the weight function is related to the land–sea distribution on the planet.

We will soon encounter distinct modes that have the same eigenvalue. In that case, we say there is an -fold degeneracy, where is the number of distinct modes with that same eigenvalue. Two or more modes with the same eigenvalue have the property (easily shown by insertion in the eigen-equation) that is also an eigenfunction with eigenvalue where is the eigenvalue common to the two eigenfunctions. Knowing this, we can form two linear combinations (choices of and ) that are mutually orthogonal and properly normalized. In this case, we can choose and . Or


where and are the new orthonormal modes that have been orthonormalized.

7.3 Linear Stability Analysis

Once we have found that there are multiple solutions for the climate for a given value of the solar constant, we are obliged to examine the stability2 of the solutions.3 In this chapter, we first examine the linear stability of the solutions and find a slope-stability theorem analogous but not identical to that in Chapter 2. Next we construct and examine the potential function[al] for the one-dimensional problem.Consider now the EBM defined by

where we have set to simplify the notation. Also, we have used as our dependent variable, which simplifies the algebra in what follows. Note that this change means that the diffusion coefficient is different from the one where is an dependent variable, that is, , where the subscript T indicates the case where is the independent variable. This difference will have no effect on our stability analysis, as is the case where we have set . We have the ice-line constraint from Budyko,


where C) and subject to the usual north–south symmetric boundary condition


We are including a -dependence in to make the results as general as possible. As in Chapters 5 and 6 and the previous section of this chapter, we expand into the eigenfunctions of the total convergence-of-flux operator of (7.24). It is an example of a Stürm–Liouville System, whose properties were outlined in the previous section.


where the are the Stürm–Liouville orthonormal eigenfunctions of degree , and the eigenvalues are real and positive. We define


The linear stability analysis begins by assuming the temperature field is perturbed by a small amount :


and the ice cap is similarly perturbed.


The equilibrium values satisfy


Inserting the small departure definitions and retaining only the linear terms, we have

where indicates the derivative of evaluated at where is the value of at equilibrium. Next, we use the relation

with the constraint to find


where we have now suppressed the argument of to mean its value at and , indicating the derivative of evaluated at .

The last expression becomes


In the last formula and in what follows, we suppress the argument . We may substitute the last equation into (7.33) and after making the small-departure assumption that


as in Chapter 2, we find the eigenvalue problem




To obtain the last equation, we assumed that is a step function with discontinuity at ,


The stability of the system is determined by the sign of the eigenvalues . Because is real and symmetric, all eigenvalues are real and bounded from below.4 If the lowest eigenvalue is negative, the system is unstable; grows exponentially in time. By casting the eigenvalue problem into a different form, we can determine the sign of the lowest root. From this point in the proof, all arguments are suppressed with the understanding that the function in question is to be evaluated at .

To determine this sign, we rearrange (7.37) and use (7.38) to obtain


Dividing this expression by , multiplying by , and summing over we obtain


This relation is a transcendental equation that is satisfied for certain discrete values of , the stability eigenvalues. By further rearrangement, we arrive at the sign of the lowest eigenvalue. Using the definition of , we find


which, with the equilibrium condition, becomes



this leads to

Finally, substituting (7.46) into (7.44) we obtain

Figure 7.3 shows a plot of the LHS and the function on the RHS as functions of the stability parameter . Intersections of the curve and the horizontal line indicate roots of the equation above and are therefore the eigenvalues corresponding to the stability problem. If all roots are positive, the climate corresponding to is stable. However, if even one root is negative, the climate at that point on the operating curve will be unstable. We see from the graph that as long as the slope of the operating curve ( vs ) is positive, the climate will be stable. But if the dashed line falls below the horizontal axis (meaning there is a negative slope of the operating curve), the left-most root in the figure will become negative, indicating instability. This is referred to as the slope-stability theorem. Note that it differs from the case treated in Chapter 2 in which theoperating curve versus was used to calculate the slope. In the one-dimensional model, the slope must be calculated from the versus curve.

Figure 7.3 Schematic graph of the left-hand side of the stability equation versus the right right-hand side, where are taken to be positive quantities for all . The are taken to be the eigenvalues for a constant diffusion coefficient model (). Intersections of the horizontal line (nominally ) with the continuous curves represent roots of the system and therefore eigenvalues . The horizontal line will lie above the abscissa so long as the slope of the operating curve is positive. All these positive roots imply stability. When the slope of the operating curve is negative (horizontal dashed line below the -axis), there will be one root which is negative, leading to instability.

Another proof of the slope-stability theorem somewhat simpler in complexity than that presented here is given by Shen and North (1999). The most general proof is by Cahalan and North (1990). Drazin and Griffel (1977) examined the solutions in the case where the hemispheres might not be constrained to be symmetric. Since the length scale is short compared to the circumference of the planet, an ice cap in one hemisphere cannot feel the effect of an ice cap in the other. Indeed, they found this to be the case by careful numerical solution of the problem (using the “shooting method” in which one fits the polar condition in one hemisphere and adjusts the ice cap in the other by trial and error until a match is made). Our proof of the slope-stability theorem would not hold in such a situation, as for some solutions the ice-cap size in one hemisphere is not the same size as that in the other.

One wonders what the stability conditions would be if there is nontrivial geography as in the next chapter. This problem has not been addressed to our knowledge except incidentally during numerical solutions in paleoclimatology. We do not even have a proof as yet for the case of a perturbation of the ice-cap edge that has wave number above zero for the homogeneous planet case. Such proofs await the creative efforts of others.

7.4 Finite Perturbation Analysis and Potential Function

In analogy to the potential function analysis of Chapter 2, we present here its generalization to the case of one-dimensional climate models (North et al., 1990). In this case, we must deal with a functional instead of a function, as the potential will depend on the function at each point in the interval . We can best describe what is meant by a functional by directly proceeding with our example. Consider the integral


where for the sake of compactness we have introduced subscripts for partial derivatives ,


Here we have chosen to write as a function of and , for example,


where is the unit step function,


Note that as


Consider the functional for an arbitrary small variation ; that is, is replaced by to form . After subtracting , we have


where we used only the first-order terms in Taylor expansions of and in . Now noting that


we integrate the first term above by parts (endpoint contributions vanish) to obtain


where the prime denotes differentiation with respect to . If the functional is to be stationary for an arbitrary but small variation , it must vanish. Since the functional form of is arbitrary, the quantity in brackets in the integrand above must vanish. The latter is just the expression for the energy balance equation. In other words, our functional is a quantity which is a local extremum when the energy balance equation is satisfied.

It is easier to visualize the functional in spectral form, . In the spectral form, we may think of as an ordinary function of the variables ; an infinite number of variables. An extremum of may be expressed as for all .

Substitution of the spectral form into the definition of the functional leads to



For simplicity, we have taken and , the values of coalbedo over ice-free and ice-covered surfaces, to be constants. It is understood that is to be substituted for in the last expression.

The condition that the vanish simultaneously leads to


which takes us back to the slope-stability theorem.

7.4.1 Neighborhood of an Extremum

Consider a particular extremum (steady state) centered at the point

Nearby (in function space), we may write , or the deviation may be written in terms of its spectral components . Expanding about the local extremum, we obtain


where the subscript zero denotes evaluation at the extremum. The terms linear in vanish because vanishes at the extremum. Up to the terms considered, is locally a quadratic in . The matrix elements,


are the structure constants for the quadratic geometric surface, . If all eigenvalues of are positive, the surface is concave upward; if one or more of the eigenvalues are negative, the surface is locally a saddle point. We proceed to show that these eigenvalues are the stability eigenvalues studied earlier.

First note that if the temperature field is allowed to be a function of time, then by following the approach for the zero-dimensional models, we have


that is, the time derivative is given by the gradient in this multidimensional space. For infinitesimal departures from steady state, we set above, expand about , and obtain


The latter equation concludes the proof that the local geometrical structure constants of yield the stability eigenvalues for that particular steady state. The last equation is the analog of the simple zero-dimensional equation (2.90) in Chapter 2.

Finally, as a conclusion to this section, consider the time behavior of the value of when the point is governed by the time-dependent energy balance equation


where we inserted the equation of motion. This latter result is the multidimensional analog of (2.91) in Chapter 2. It has a corresponding interpretation: initial departures of the state from a local extremum of lead to a trajectory of the system point such that decreases. The point will continue down the gradient of until a local extremum is found. Clearly, saddle points are unstable—if perturbed infinitesimally, the system point can leak out into some neighboring basin.

7.4.2 Relation to Gibbs Energy or Entropy

Solutions of the energy balance equation satisfy a minimum principle reminiscent of Gibbs energy in thermodynamics (when pressure and temperature are held fixed). If we take the negative for , it might be interpreted as the maximum of entropy production. Then the extremum condition would be somewhat like Prigogine's theory of nonequilibrium thermodynamic states (Prigogine, 1968). Golitsyn and Mokhov (1978) showed that the linear climate models satisfy Prigogine's conditions. The idea is tantalizing, but has so far not been satisfactorily implemented in general terms.

7.4.3 Attractor Basins—Numerical Example

It is possible to work out in some detail an example illustrating the concept of the potential surface for a two-mode model. In this case, the functional becomes the truncated version of (7.57). The function is


where, in the two-mode model, the ice edge can be expressed in terms of :


Together with the other terms in (7.57), we have enough information to plot the surface corresponding to as a contour diagram in the plane. We choose the parameters such that there is a cusp: W m, = 220 W m, = 1.90 W m K, = 0.68, = 0.38, = 0.30. The operating curve for this system is shown in Figure 7.4a. This choice of parameters leads to five solutions labeled , and . The potential surface is mapped in (b) with the steady-state points indicated. It can be clearly seen that is a saddle point corresponding to an unstable solution. Figure 7.4c and d shows successively higher resolution focusing on the right-hand basin. The higher resolution shows that there are two distinct minima: and with a saddle point . The ice-free solution corresponds to a very shallow minimum and a very small but finite agitation would push it over the hill () into the deeper (more-stable) basin .

The relative minima in a purely dissipative system such as those we are studying here are isolated points called attractors. It is possible to map out the attractor basins for this class of problems as illustrated by the two-mode approximation with the basins shown in Figure 7.4.

Figure 7.4 Potential function for a two-mode, one-dimensional nonlinear climate model. (a) The operating curve indicating five solutions (labeled ) for . (b) A coarse resolution plot of the potential function in the plane. (c) and (d) Successively higher resolution around the right-hand basin which turns out to have two local minima. The relative minima are stable climate states, the saddle points correspond to unstable states. Further discussion in the text.

7.5 Small Ice Cap Instability

Consider the operating curve in Figure 7.5. This two-mode operating curve has five equilibrium solutions for W m. Because of the slope-stability theorem, only the ones with positive or zero slope are stable. The cusp at the top of the graph is magnified in the right-hand panel as the solid line. To illustrate the sensitivity of the cusp, we can vary the coalbedo of ice cover. According to Figure 7.6 (take the solid curve parameters), the operating curve has a negative slope for values of , which corresponds to a latitude of about 63. The negative slope of the operating curve implies that polar ice caps whose radius is smaller than about 25–37 on a great circle will be unstable. This point is worth examining more closely as it may have implications for paleoclimate and for future climates. For example, the Antarctic ice cap is of about this size and the seasonal sea ice in the North Pole area is also of roughly this size. Is it on the verge of unstable collapse if the planet is heated slightly by the greenhouse effect? Moreover, did it form suddenly? What parameters control the size of the least stable ice cap?

Figure 7.5 An operating curve with five equilibrium solutions for W m. In this illustration, the , = 210 W m; B =1.90 W m K, and the coalbedo is taken to be flat as function of latitude with and .

Figure 7.6 Illustration of the sensitivity of the cusp's shape to small changes in the coalbedo of ice: The solid curve is for the same values as in Figure 7.5 (). The long-dashed (leftmost) curve is for dark ice (), while the rightmost curve is for bright ice ().

7.5.1 Perturbation of an Exact Ice-Free Solution

Let us imagine a situation where we are on the ice-free branch of the solution and slowly lower the solar constant until we are at . Our solution will constitute a point on the upper branch of Figure 7.6. The model is linear in the (infinitesimal) neighborhood of this solution as there is no ice cap. Using the parameters of Figures 7.5 and 7.6, we find that the temperature is given by


At every point, the temperature is above 10 C (Figure 7.7). At the pole, the temperature is 9.5 C, which is close to the critical value. Consider next the effect of adding a heat source centered at the pole whose density is (). Let the departure from the ice-free solution be . It satisfies


This equation can be solved by use of Legendre polynomial series on the sphere, but it is more instructive to look at the solution on the plane tangent to the pole where solutions in closed form can be obtained. Let be the plane polar coordinate corresponding to the distance from the pole. And let be the response function to the thermal perturbation . Then satisfies


where the length scale is defined as .

For a small patch of ice (diameter ), the solution5 is


where is the radius of the ice patch and is the zero-order modified Bessel function which is pictured as the solid curve of Figure 7.8. This particular Bessel function looks like a decaying exponential and in fact has the asymptotic form for large . The important point is that the influence of the small patch extends from its edge a distance the order of . It is hardly surprising that stable ice caps smaller than about in radius did not exist in model solutions.

Figure 7.7 Plot of the temperature near the pole as function of , the cosine of the polar angle in the case where the temperature at the pole is just above the critical value of 10 C. The value of the curve ; at the pole, this yields C. What happens if a small patch of ice is placed at the pole?

Figure 7.8 Graph of the modified Bessel function (solid line) as a function of and an asymptotic approximation, (dashed line). The Bessel function is the response to a point source in plane polar coordinates. The figure shows that it is nearly proportional to for large . It partially answers the question posed in the previous figure.

Imagine the no-icecap solution to be in place; this state's temperature (C) is above the critical value (C). Now “by hand,” we add a small ice patch at thepole. The patch as a source of cooling has an influence a distance away from the pole. If the albedo of the ice is high enough, the patch can pull the temperature down below critical over a distance comparable to . This means that two solutions can exist, one with no ice and another with a patch no smaller than approximately . This is the underlying meaning of the small ice cap instability. More details can be found in North (1984).

Figure 7.9 Solid line: The magnitude of the complex modified Bessel function as a function of , where for and (the low frequency limit). The dashed line indicates how the length scale is shortened for higher frequency forcing. In this latter case, . The length scale tends to zero as the frequency increases further.

7.5.2 Frequency Dependence of the Length Scale

Next consider the solution to a small scale source in the heating problem, when the strength of heating has a sinusoidal time dependence: . Then we should use the time-dependent version of the energy balance equation:


By the assumption that has the same sinusoidal dependence as , we find the same result except that is replaced by


or its magnitude is


To see the dependence of the solution, we plot the magnitude of in Figure 7.9. We see that the characteristic length scale is strongly dependent on the frequency of the forcing as the forcing frequency is raised to . Over this interval of frequency, the length scale shrinks by a factor of 2. It continues to shrink toward zero as the frequency is increased. The magnitude of the Bessel function decreases similarly to an exponential as a function of distance from the origin, but the real and imaginary parts of the function oscillate as dampened pulses diffuse away from the origin.

An application of the shortened length scale for oscillating point sources is in the seasonal cycle. In this case, over land is about 1/12 year and year; hence with respect to a land surface, the seasonal cycle is very low frequency and the low frequency limit is a reasonable approximation. On the other hand, over ocean, is about 5 years and year, hence, we are far from the low frequency limit over oceans. The consequence is that in seasonal cycle applications the reach or fetch of a seasonally oscillating source is only a fraction of its lowfrequency or static value. During the seasonal cycle, the oceans and land surfaces have different lags behind the heating. The above considerations mean that a neighboring lagging ocean surface has an influence far inland into a continent. But the influence of the land temperatures on neighboring ocean surfaces extends only a fraction of that distance off the shore. A larger study of the frequency dependence of the length scale can be found in North et al. (2011).

One caveat about the small ice cap instability is its sensitivity to the parameters chosen. Even tiny changes in the parameters of the model can cause it to appear or not in the operating curve. One might conclude that it is an artifact of such a simple model. However, it might be indicative of phenomena that do actually occur in the climate system. The fact that it is so sensitive to the parameter values suggests that it might be there at some times in history and not at others depending on the environment at the time.

7.6 Snow Caps and the Seasonal Cycle

The latitude of the steady-state snow line varies through the seasons following roughly the freezing line. By steady state of course, we mean the ensemble average over many seasonal cycles as before. We postpone the issue of the snow-line realizations in the presence of climate fluctuations. The seasonal progression of the snow line and its modification under changed external forcing conditions has important consequences as we will see. A very important issue is whether snow lingers through the summer at a given location, say the pole. If this does happen, there is the possibility of growing an ice sheet if there is persistent snow over land.

First consider an all-ocean planet. In this situation, the seasonal cycle is small and the pole will remain below freezing all year round. The Arctic Ocean may be an analog to this case and it is interesting to note that there is perennial sea ice at the North Pole and that this may be suggestive of a small ice cap poised for abrupt removal under the influence of a warm forcing such as the greenhouse effect or low frequency changes in oceanic conditions. Removal of perennial sea ice at the North Pole would have negligible direct thermal effects on the planet but it might cause profound changes in the circulation of the Arctic Ocean and consequently the formation of deep water in the North Atlantic.

Next consider an all-land planet. It is easy to show that, in this case, with the present obliquity (tilt of the Earth's spin axis with respect to the ecliptic plane) the summers are very hot at the poles (exceeding 30 C). Hence, for the all-land planet there will be no persistent snow at the poles and the planet is much too hot for a small ice cap summer solution.

7.7 Mengel's Land-Cap Model

An interesting intermediate case is that of a symmetrical cap of land at the pole surrounded by open ocean (Mengel et al., 1988). The analog of this case is the continent of Antarctica whose size is about the size of an annual cycle length scale in the simple climate models.

As an introduction to problems with partial land and sea geographies, consider a planet with a disklike island centered at the pole. The rest of the planet is ocean covered. If the solar irradiance is in a certain range, the temperature will fall below freezing in winter in which case the albedo will become higher in the snow-covered areas. This is a model that still has only one horizontal dimension, the latitude. The way of implementing the land–sea contrast is to introduce an effective heat capacity that is small over land and large over ocean. This means the coefficient of the term has a strong discontinuity at the land cap's edge. Solving this boundary value problem analytically is not likely to be possible; hence, one must resort to numerical methods. Other problems include how to handle the heat capacity when an area of ocean surface is ice covered. In this initial study, this latter effect has been ignored.

Mengel et al. (1988) used an expansion of the temperature field into Legendre modes. Because of the discontinuity in albedo and heat capacity, many modes (29) had to be retained in the expansion, but experiments with half that many did not reveal different results. Details of the methodology can be found in the paper.

Figures 7.10 and 7.11 show results of numerical calculations for such a zonally symmetric planet. The figure shows that for warm values of the solar constant, the snow line retreats rapidly to the pole in spring and reappears in fall, moving well out into the ocean areas in winter. In other words, we have a snow-free summer over the land cap. However, as the solar constant is lowered, an abrupt transition occurs to a climate configuration with persistent snow throughout the summer. This transition is a form of the small ice cap instability but is more complicated than the simple static ice cap studied earlier in this chapter. It is clearly affected by the seasonality and the geography. The smallness of the land cap makes it possible for the oceanic influence to stretch all the way to the pole, keeping summers marginally cool enough to support ice in summer if the snow albedo is present. On the basis of this model finding (and it persists in general circulation model (GCM) simulations), one might expect that the transition to full glaciation on Antarctica might have been rather sudden as the parameters passed slowly over the critical point. Of course, in reality, it was not the changing solar constant that induced the Antarctic transition but some other parameter (continental drift) acting in a similar way. We will return to this paleoclimate question in Chapter 8.

Figure 7.10 Illustration of the extreme sensitivity of the seasonal cycle of the snow line to solar brightness when a bifurcation (tipping point is near). Shown is a plot of the seasonal variation of the snow line for a land-cap planet (analogous to Antarctica) for two neighboring values of normalized solar constant. The shaded area indicates land area. The ordinate is the latitude in degrees and the abscissa is time in years. The origin is at NH winter solstice.

(Mengel et al. (1988). Reproduced with permission of Springer.)

Figure 7.11 Similar to Figure 7.10, except that the variable plotted is the temperature at the pole. Shown is a plot of the seasonal variation of the polar temperature as a function of time of the year. The shaded area indicates land area. It shows the temperature region above freezing. The origin is at winter solstice. The values of total solar irradiance for this pair of curves coincides with those of Figure 7.10.

(Mengel et al. (1988). Reproduced with permission of Springer.)

An example of counterintuitive behavior can be found in Figure 7.12. Here the model was initialized on one side of the bifurcation but the parametric conditions were on the other side. The relaxation to the seasonal cycle appropriateto the parametric conditions did not occur with the linear timescale of a few years, but actually took tens of years. This is because the solution and parametric conditions were very close to a bifurcation where nonlinear effects are dominant. This is a problem often ignored in climate simulation experiments. For example, we may well have already passed the bifurcation for sea ice-free summers in the Arctic, but the planet is still adjusting to this condition. Similar long-delayed adjustments may be underway in Greenland and Antarctica.

Figure 7.12 Substantial modification of the relaxation time to the steady-state seasonal cycle by the nearness of a bifurcation (tipping point). Shown is a plot of the time series for two realizations of the model for initial conditions near a singularity or bifurcation in the solution. (a) Initially the model is set to be in the icy pole case, but suddenly the solar constant is switched to the ice-free pole case. The time constant for the transition in this case is of the order of 43 years, but it is sensitive to the initial condition. (b) The opposite case where the solution is initially for the ice-free pole, but the total solar irradiance is suddenly switched to the icy pole value. This time the adjustment takes about 120 years.

(Mengel et al. (1988). Reproduced with permission of Springer.)

The work of Mengel et al. (1988) (see Figure 7.12) was extended by Lin and North (1990) to several other geographies (see also, Huang and Bowman, 1992). The techniques used by these authors were also different and that might be interesting for some readers. Lin used the Fourier harmonic technique to solve for the steady seasonal cycle. This turned out to be more efficient. The nonlinearity was taken into account by iteration. First, Lin demonstrated that she could repeat Mengel's work for a disklike land mass centered at the pole. She investigated how the bifurcation occurs as it is expressed in an operating curve, shown here in Figure 7.13. Using the same method, Lin investigated the case of a planet with a polar ocean but several cases of zonally symmetric bands of land configurations. Figure 7.14 shows an example of a band of land configuration where the land surface lies between 50 and 75 latitude. In this case, as indicated in Figure 7.14, there are two bifurcations, one with an ice-free Arctic Ocean as well as snow-free land for warmest values of solar irradiance. As the solar irradiance is lowered to 0.9316 of the present value, the former ice-free polar ocean exhibits summer ice cover, and the land surface is nearly snow free for most of the summer. Finally, at solar irradiance 0.9255, the entire polar area is ice covered year round.

Figure 7.13 Graph of the operating curve for the north polar temperature (time = 0 at winter solstice). Operating curve for the polar temperature at summer solstice for the Mengel's land-cap model (analogous to Antarctica).

(Lin and North (1990). Reproduced with permission of Springer.)

Figure 7.14 The complicated bifurcation structure of a seasonal cycle of a “band of land” model. Similar to previous figures except for a different land/sea configuration. Snow-line edge versus time of year (in tenths). The shaded latitudes are land, unshaded are ocean. Winter solstice is at . Note that there are two bifurcations for this geography. The numbers beside the snow-line curves indicate the total solar irradiance in ratio to its present value for the particular case.

(Lin and North (1990). Reproduced with permission of Springer.)

The Mengel-like experiment has also been examined in a GCM (the Genesis Model of the early 1990s) by Crowley and Baum (1993). They used a geography from the Carboniferous period which was similar to the continent at the pole configuration in the Mengel et al. (1988) work (Figure 7.15).

Figure 7.15 Bifurcation in the operating curve for the Genesis GCM with a Carboniferous land/sea distribution that is roughly comparable to a land mass centered at the South Pole. Both EBM and GCM show summer snow lines that have a discontinuous break at a certain value of the solar constant (shown here as , the ratio to its value in the Carboniferous period.

(Crowley and Kim (1994). Reproduced with permission of AAAS.)

7.8 Chapter Summary

The chapter began with a discussion of the ease of studying linear systems compared to the difficulty of dealing with nonlinear systems. Hence, it is wise to work linear problems when possible. In many cases, when nonlinearity is small, one can get by with such an approximation. On the other hand, the great power garnered from experience based on linear methods fails spectacularly when the nonlinear effects are large enough to cause bifurcations in steady-state solutions especially when parameter conditions bring the solution very near to a bifurcation or other form of singularity.

The first problem taken up in the chapter was the nonlinear ice-cap problem where one finds three and sometimes five solutions depending sensitively on parameter values. An interesting problem consists in finding the stability of steady-state solutions. It turns out that with ice-cap models, the stability of a solution depends on the sign of the slope of the operating curve at the solution point. This is the point that depicts the polar ice-cap edge's latitude versus the solar irradiance. If the ice cap increases with increasing solar irradiance, the solution is unstable, a result anticipated heuristically by Budyko. The nonlinear ice cap's stability properties can also be studied by means of the potential function introduced in Chapter 2, only this time it must be multidimensional (usually, infinitely so).

The idea of a length scale suggests an explanation for the cusp-like behavior of the operating curve near small ice-cap sizes. It turns out that because a tiny patch of high-albedo material such as snow can induce cool surface air over large regions. If the albedo change is large enough even for a small patch, the surrounding region up to a radius of the length scale can become below freezing. Hence, ice caps can form.

The chapter closes with a seasonal model with a mix of land and sea surface. It requires a new addition to EBCMs, the position-dependent effective heat capacity. Position-dependent heat capacity is to be explored in more detail in Chapter 8. This addition keeps the model linear but it was introduced here to solve an interesting seasonal model, the Mengel model, which includes a strong nonlinear albedo feedback, owing to the high albedo of seasonal snow cover. The first model explored is one with a cap of land centered at the pole analogous to Antarctica. The model solutions are for a steadily repeating seasonal cycle that includes a snow line hooked to the freezing point. The model solutions indicate a strong bifurcation at a particular value of total solar irradiance. The seasonal cycle warmer than a critical value of the control parameter leads to snow-free summers over the polar land cap. Values lower than the critical value lead to snow cover large enough to cover the polar continent.

A second model by Lin and North (1990) expands the work to include more geographical configurations including zonally symmetric bands of land. Again interesting bifurcations are found with even more richness of behavior. This effect brings to mind a theory for how glaciers or ice sheets might initiate. When the control parameter is below critical, the summers are perpetually snow covered, allowing the accumulation of ice as the temperature never goes above freezing. Another effect discovered in these models is an explicit demonstration that relaxation times near one of these bifurcations can exceed decades, departing by an order of magnitude from the linear timescales of EBCMs. There were a number of experiments on the two-dimensional model (Chapter 8) with realistic geography with nonlinear snow/ice albedo feedback Hyde et al. (1990), which will be discussed in Chapter 12.

Notes for Further Reading

There are many books on nonlinear systems, but a particularly good one is that by Drazin (1992).


  1. 7.1 On the basis of (7.11), find the closed form for as a function of for and .
  2. 7.2 On the basis of your answer in Exercise 7.1, determine the insolation for which the ice-cap edge is at . Use W m, W mC), and .
  3. 7.3 Consider a time-dependent, one-dimensional energy balance model with unit heat capacity in the form

    where is the outgoing radiation flux density, which is a function of latitude and time, is the latitude-dependent diffusivity and is the coalbedo depending on both latitude and the latitude of the ice-cap edge. Assume further that

    namely, there exist orthonormal eigenfunctions and eigenvalues satisfying

    Then, the outgoing radiation field can be expanded in terms of eigenfunctions , that is,

    Solve the energy balance model to determine the radiation (and therefore the temperature) as a function of , in closed form. Find the equilibrium solution for a given ice-cap-edge location .

  4. 7.4 In Exercise 7.3, imagine that the equilibrium infrared radiation to space, , is perturbed by a small amount ; that is

    As a result of this perturbation, the latitude of the ice-cap edge changes slightly by an amount , that is,

    1. a. Show that the perturbed temperature, to a first-order approximation, satisfies
    2. b. Show that the perturbation of the ice-cap edge is determined to be
  5. 7.5 In Exercise 7.4, a linearized equation for a small temperature departure, , for a one-dimensional nonlinear model is shown to satisfy


    and the coalbedo is defined as a step function with a jump



    Assume a solution of the form

    show that

  6. 7.6 Show that
    1. a.
    2. b.
    3. c.
  7. 7.7 Plot the right-hand side of Exercise 7.6(c), and find the solutions as a function of the value and sign of the left-hand side . Interpret the result in terms of the stability of the solution .