One-Dimensional Energy Balance Model

So far, we have developed a zero-dimensional model of the earth. If desired, we could continue to refine this model by including more complicated functions for some of the parameters (such as time and temperature dependent terms in the greenhouse effect parameters or the albedo

The next logical step is to allow spatial resolution into our model. We will do this by breaking the northern hemisphere into 9 latitude regions of 10o each. The general geometry is shown in the following diagram:

Global Climate Model
Reprinted with Permission From: A Climate Modeling Primer, A. Henderson-Sellers and K. McGuffie, Wiley, pg. 58, (1987)

In the model we will consider, the solar flux and albedo both have a latitude dependence (Si and α i, where i runs from 1 to 9 indicating the latitude band). Also, if a latitudinal band is colder or warmer than the average global temperature, heat will flow into/out-of the region. We will assume that this heat flow depends linearly on the temperature difference between the region and the global temperature, in other words it is: F×(Ti-TAve). Putting this additional term into Equation 3 from the Global Climate Modeling Project results in an energy balance equation of:

(Eqn. 4): ΔT = (PGain-PLoss)Δt/ CE
In this case, for each of the i=1..9 regions
PGain = Si(1-αi)/4
PLoss =A + B×Ti + F×(Ti-TAve)

The equilibrium limit of this equation (setting ΔT=0) is:

(Eqn. 5): Ti = [Si(1- α i) + FTAve - A]/(B+F)

The Matlab script: ebm1.m for this model is available. For this model, we will use the following parameterizations:
SiSolar Flux in Latitude Band i
This is the product of S/4 (the average global solar constant) times the insolation si
siSolar Insolation
The fraction of solar flux incident on each latitudinal band.
α iAlbedo in Latitude Band i
The albedo of ice is much larger than the albedo of land/water. We can do a crude model of the temperature dependence of the albedo by using α i=0.3 for Ti > Tc or α i=0.6 for Ti <= Tc
TcBelow this Temperature, we assume a Permanent Ice Pack (Tc = -10oC)
FHeat Transport Coefficient (F=3.80 W m-2 oC-1)
TAveGlobal Weighted Average Temperature
This temperature is the weighted average of the temperature in all of the latitude zones on the previous iteration. The weighting factors fi are just the relative fraction of the surface area of the sphere in each latitude zone.
A & BCoefficients expressing Infrared Radiation Loss (A=204 W m-2 and B=2.17 W m-2 oC-1)
CEHeat Capacity (CE = 2.08 × 108 J/m2 oC)
If we want to find a steady state solution using Eqn 5, we will need to iterate the program from an initial temperature distribution until it finds the final state. This is required because the heat flow in each latitude zone depends on TAvg, but TAvg depends on the temperature distribution. By having the program cycle through several times, it should converge to a steady-state solution.

The Rough Flow of the Matlab script ebm2.m (which iterates Eqn. 5) is:

  1. Initialize the physical parameters for the system
  2. Call the keyboard command to allow you to easily change any variables
  3. Initialization finishes up
  4. Set MaxTempDiff=1E6 - we will be checking for a small value in the next step
  5. Start a while loop - this is terminated if the Maximum Temperature difference (MaxTempDiff) is smaller than the tolerance specified by TolTempDiff
  6. Calculate the albedo using the desired algorithm
  7. Calculate the new latitudinal temperature distribution
  8. Calculate the Global Average Temperature (using appropriate weighting factors)
  9. Calculate the maximum of all of the temperature changes.
  10. Repeat the "while" loop until complete
At the end of this process, you will have a list of the temperature as a function of time.

EXERCISE 4: Run the Matlab script using the parameters given in the program and this handout. Try adjusting the factor SX (which is multiplied by the Solar Constant) to find out how sensitive the model is to variations in the solar constant. By what factor does the solar constant need to decrease before the earth is completely glaciated? Before all the permanent ice pack melts? Modify your program so instead of finding a steady-state solution it will find a time dependent solution (using Eqn. 4).

There are numerous extensions to the basic model presented which you can explore. In particular, some questions you may want to answer include
This table [from S.G. Warren & S.T. Schneider, J. Atmos. Sci. 36, 1377-91 (1979)] contains some actual measurements - this may be useful in extending your models:
Zone Annual Mean
Temp (oC)a
Solar Insolation
(Fraction of
Solar Constant) b
α i
Heat Transfer
Into/Out of
Zone (W/m2)c
80-90o N-16.9 0.500 0.589 -103
70-80 -12.3 0.531 0.544 -94
60-70 -5.1 0.624 0.452 -72
50-60 2.2 0.770 0.407 -47
40-50 8.8 0.892 0.357 -21
30-40 16.2 1.021 0.309 1
20-30 22.9 1.120 0.272 18
10-20 26.1 1.189 0.248 46
0-10o N 26.4 1.219 0.254 59
0-10o S 26.1 1.219 0.241 56
10-20 24.6 1.189 0.236 41
20-30 21.4 1.120 0.251 22
30-40 16.5 1.021 0.296 0
40-50 9.9 0.892 0.358 -27
50-60 2.9 0.770 0.426 -57
60-70 -6.9 0.624 0.513 -86
70-80 -29.5 0.531 0.602 -90
80-90o S -42.3 0.500 0.617 -88

a Data from S.G. Warren & S.T. Schneider, J. Atmos. Sci. 36, 1377-91 (1979)
b Data from P. Chylek & J.A. Coakley, Jr., J. Atmos. Sci. 32, 675-9 (1975)
c Data from J.S. Ellis & T.H. Vonder Haar, Atoms. Sci. Paper 240, Colorado State University, 50pp. (1976)

Electronic Copy:
Revised: 14-JUL-97 by Tom Huber, Physics Department, Gustavus Adolphus College.