# 40. Playing with a diffusive energy balance model

Posted on August 8th, 2013 in Isaac Held's Blog Latitude of ice margin as a function of a non-dimensional total solar irradiance $q$ in the diffusive energy balance climate model described by North 1975, for different values of the non-dimensional diffusion $d$.  Stable states are indicated by a thicker line.

When we were first starting out as graduate students, Max Suarez and I became interested in ice age theories and found it very helpful as a starting point to think about energy balance models for the latitudinal structure of the surface temperature.  At about the same time, Jerry North had simplified this kind of model to its bare essence: linear diffusion on the sphere with constant diffusivity, outgoing infrared flux that is a linear function of surface temperature, and absorbed solar flux equal to a specified function of latitude multiplied by a co-albedo that is itself a function of temperature to capture the different planetary albedos for ice-free and ice-covered areas.  Playing with this kind of “toy” model is valuable pedagogically —  I certainly learned a lot by building and elaborating this kind of model — and can even lead to some nuggets of insight about the climate system.

Using the same notation as in post #36, $C \partial T/\partial t = \nabla \cdot C\mathcal{D} \nabla T - (A + B (T-T_0)) + \mathcal{S}(\theta) \mathcal{A}(T)$. $\mathcal{S}(\theta)\mathcal{A(T)}$ is the absorbed solar flux, where $\mathcal{A(T)}$ is the co-albedo.  We set $S = S_0 \mathcal{\sigma}(\theta)$, where $S_0$ is the global mean incident flux = total solar irradiance/4 so that $\mathcal{\sigma}(\theta)$ averages to unity over the sphere.  The form $\sigma(\theta) = 1 - 0.5P_2(sin(\theta))$, with $P_2(x)$ the second Legendre polynomial $(3x^2-1)/2$ is a pretty good fit to the annual mean incident solar flux as a function of latitude $\theta$. The seasonal cycle is ignored here. $A + B(T-T_0)$ is the outgoing infrared flux linearized about the reference temperature $T_0$, $C$ is a heat capacity, and $\mathcal{D}$ is a kinematic diffusivity.   We can just set $T_0 = 0$ (and think of it as the freezing temperature when we set the albedo). We can define a non-dimensional temperature, $\mathcal{T} \equiv BT/A$, diffusivity $d \equiv C D/Ba^2$, and mean incident solar flux or total solar irradiance, $q \equiv S_0/A$.  Finally, we can use the simplest possible albedo formulation that provides some ice-albedo feedback: $\mathcal{A}(T) = \beta$ for $\mathcal{T} > 0$ and $\mathcal{A}(\mathcal{T}) = \alpha$ for $\mathcal{T} < 0$ .  I use $(\alpha,\beta) = (0.4, 0.7)$ for the results described here.

Our equation for steady state solutions independent of longitude, writing out the divergence of the diffusive flux in spherical coordinates, is now $\mathcal{T} -\frac{1}{\cos(\theta)}\frac{\partial}{\partial \theta}( \cos(\theta) d(\theta) \frac{\partial \mathcal{T}}{\partial \theta}) = q\sigma(\theta) \mathcal{A}(\mathcal{T}(\theta))-1$.

(Aug 12: Fixed a couple of typos in the last few days; hopefully OK now.) If both the diffusivity and the albedo are chosen to be spatially uniform one can solve this equation analytically for this specific choice of $\sigma(\theta)$ because $\mathcal{T}$ is then a constant plus a term proportional to the second Legendre polynomial $P_2$ which is an eigenfunction of the Laplacian.  Comparing the result with $d=0$ with that for non-zero $d$, one finds that the presence of diffusion reduces the equator-to-pole temperature gradient by the factor $1/(1+6d)$.  One needs this reduction to be about a factor of 2-3 to get a reasonable temperature gradient, which translates into a value of $d$ of 0.2-0.3.

We look for solutions that are symmetric about the equator and have temperatures below freezing poleward of a given latitude, the” icecap edge”, and above freezing equatorward of this latitude.  The resulting ice edge as a function of $q$ for different values of $d$ is shown in the figure at the top of the post.  Ice-free states are indicated by the horizontal line at $\theta = 90$ and ice-covered “snowball Earth” states by the horizontal line at $\theta = 0$.  Partially glaciated steady states also exist, some of which are unstable.  The branches on which the ice edge moves equatorward with increasing $q$ are unstable, not surprisingly. There is a small ice cap instability, with ice caps smaller than this threshold receding unstably to the ice-free state.  And there is a large ice cap instability beyond which the ice grows unstably due to a runaway albedo feedback, until one reaches the snowball state.  This kind of model attracted considerable attention because of the rather small range of solar flux for which partially glaciated states exist and the proximity of these state to the large ice cap instability, for plausible values of the diffusivity.  As a point of comparison, Voigt and Marotske 2010, using a modern comprehensive coupled atmosphere-ocean GCM, find that a reduction of 6-9% in the solar flux is sufficient to generate the large-icecap instability. This is obviously an interesting number.

The small ice cap instability in this simple model captures the basic idea that an icecap has to have a certain size to protect itself from diffusion of warm air from lower latitudes.  The critical size increases with increasing diffusivity. But this small ice cap instability (unlike the large ice cap instability) turns out to be fragile to modifications to the model such as smoothing the temperature dependence of the albedo or adding a seasonal cycle or adding some noise.  To decide if the Arctic ice possesses a small icecap instability requires much more realistic atmospheric and sea ice models.  But this simple model does get you thinking about the importance of heat transport from lower latitudes for this issue.

In post #36 there are some plots indicating that the effective diffusivity for heat in the atmosphere should be thought of as having a maximum in midlatitudes.  To mimic this schematically, I have set $d = 0.4$ for $40^\circ < \theta < 60^\circ$ and $d = 0.2$ elsewhere.  The solution is shown by the black line in the figure below. (The results for uniform values of $d = 0.2$ with $d =0.4$ are copied over from the figure at the top for comparison.) A new instability has been created, with no stable ice caps ending within the region in which the diffusivity has been given the larger value of 0.4.  Comparison with the uniform $d = 0.4$ case indicates that it is not simply the magnitude of the diffusivity that creates this instability but rather its horizontal structure. David Linder, Max, and I touched on this kind of behavior in an old paper Held et al 1981, Some related results are discussed in Rose and Marshall, 2009, coming from the direction of trying to include the effects of ocean heat transport in an energy balance model. This kind of stability diagram could generate some interesting hysteresis loops from a time-dependent parameter like the obliquity in Milankovitch ice age theories.

Irrespective of any imagined relevance for the ice age problem, I am interested in closing the gap between this kind of diffusive model and GCMs.  In particular, it would be interesting to take an aqua planet atmospheric GCM over a slab ocean with some heat capacity but no horizontal oceanic heat fluxes, no seasonal cycle, no sea ice, but with a simple specified surface albedo as a function of temperature, and map out its behavior as a function of the incident solar flux (with and without cloud feedback).  You could then try to mimic this idealized  but still chaotic and turbulent GCM’s behavior with steady state energy balance models incorporating theories for the atmospheric heat flux. Points of interest include how an advancing ice edge affects the effective diffusivity and how best to represent cross-equatorial influences.

[I think it is an excellent project for beginning students to generate these solutions on their own.  You can reintroduce the time-dependence and integrate forward in time, but this will only give you the stable states.  For this simple case of a step function albedo, there is an easy way of getting all of the states shown: specify the ice edge and, therefore, the albedo, then solve the boundary value problem directly for the steady state by inverting the tridiagonal matrix that you get from simple finite-differencing.  The solutions will not have the $T = 0$ point coincide with the ice edge, so you then need to iterate and find the value of $q$ that gives you consistency between temperature and albedo.]

[The views expressed on this blog are in no sense official positions of the Geophysical Fluid Dynamics Laboratory, the National Oceanic and Atmospheric Administration, or the Department of Commerce.]

## 5 thoughts on “40. Playing with a diffusive energy balance model”

1. Steve Fitzpatrick says:

A toy model can be interesting (and even informative). But in this case, I wonder how informative it is. If I remember correctly, there is strong correlation between albedo in the tropical Pacific and ENSO, with La Nina conditions (cooler tropical pacific surface water) corresponding to significantly lower tropical Pacific albedo. In my mind, this argues that a cooing tropical ocean would substantially lower Earth’s net albedo, and so help to avoid the unstable growth of ice caps. (A warmer tropical ocean would presumably do the opposite.) The long term non-snowball-Earth conditions which have existed over vast geological periods (500+ million years) suggests to me that unstable ice cap growth is at best improbable.

1. Isaac Held says:

If the snowball earth state is stable at the present value of TSI (and CO2) then from the stability diagram there must be an unstable branch between the present climate and the snowball state. For a recent paper on simulations of the snowball state, see Abbott et al 2013. The climate of the snowball state is so extreme, with almost no water vapor in the atmosphere, that I can’t imagine that present day interannual variability has much to say about the stability of this state. The amount by which you have to reduce TSI to get into this state is going to be dependent on cloud feedbacks. A model like this is just meant to raise questions; it is obviously not the framework in which to talk about cloud feedbacks. if you think cloud variability during ENSO has something to say about the large icecap instability, a logical approach would be to look at GCMs that try to simulate both.

2. Chris Colose says:

Steve,

I don’t agree with your concern here.

The interannual/multidecadal variability in Earth’s albedo is quite small compared to the stuff that is trying to be teased out of this type of model. If something like the tendency for a more El Nino/La Nino-like state was a critical factor in determining whether a snowball (or hothouse) climate was set in motion, then our planet would not be habitable as we know it. In my view, paleoclimate (including the Holocene) provides compelling evidence that the Earth is not susceptible to such large bifurcations in the absence of an external, applied forcing (for example, just due to large internal variability passing some threshold and locking the climate into some new, persistent state).

Such behavior can’t be excluded on first principles, but there is no example of an unforced coupled model (between atmosphere-ocean-sea ice, etc) that spontaneously produces a change as large as a doubling of CO2, let alone a Snowball Earth, particularly in a Holocene-type background setting. Thus, there is little reason to believe that the bifurcation structure in the above figures is remarkable sensitive to the detailed representation of internal variability, or e.g., the diurnal and seasonal cycle. In fact, if the results were highly dependent on such features, then I think that would be good evidence of a model that is not too useful.

3. Aptoide says:

I agree with you Chris.

There is little reason to believe that the bifurcation structure in the above figures is remarkable sensitive to the detailed representation of internal variability. The Snowball Earth bifurcation point as well as the transition times are well reproduced by a zero-dimensional energy balance model of the mean ocean potential temperature. During the transition, the asymmetric distribution of continents between the Northern and Southern Hemisphere causes heat transports toward the more water-covered Southern Hemisphere. This is accompanied by an intensification of the southern Hadley cell and the wind-driven subtropical ocean cells by a factor of 4. If we set back TSI to 100% shortly before the transition to a modern Snowball Earth is completed, a narrow band of open equatorial water is sufficient for rapid melting.

1. Isaac Held says:

I think Steve’s point was not that internal variability directly affects the large icecap instability but that it might serve as a probe of the sign and strength of cloud feedbacks that could affect the instability. My concern was that the snowball climate, or a climate with very extensive icecaps trying to decide whether to freeze over or not, is so different from today’s that expecting ENSO to be a useful probe of the relevant feedbacks is unjustified. (Actually, I question whether they are relevant for estimating cloud feedbacks to small perturbations in CO2 about the present climate, but that is the subject of a different post.)