AN ENERGY-BALANCE CLIMATE MODEL[1]

Aad van Ulden 7 june 2001, adapted by Boeker 6 june 2002

The energy balance climate model is based on Hoffert et al[2]. It is a 1-D model, which describes the global temperature response to prescribed radiative forcings. This response is described by temperature anomalies relative to an equilibrium state. The anomalies are described by a greek D, but at the end we will omit that D . In this adaptation we ignore atmosphere and land.

OCEAN

The ocean is represented by a box with an area A_m equal to 65% of the earth surface A and a depth of 4100 m. The upper ocean is a well-mixed layer with depth h_m and uniform temperature anomaly DT_m.

Vertical transport below the mixed layer is described by down-welling, upwelling and diffusion. Down-welling occurs in small regions at high latitudes with an area A_d, where cold water sinks to the bottom of the ocean with a fixed velocity w_d. The down-welling mass flux is r_dA_dw_d, where r_d is the density of the down-welling water.

Upwelling occurs in a very large fraction A_u of the ocean at a fixed mean rate w_u. The upwelling mass flux is r_uA_uw_u. Continuity[3] requires that:

r_dA_dw_d + r_uA_uw_u = 0 (1)

In the down-welling regions the temperature anomaly DT_d is assumed vertically uniform[4]. Thus the down-welling heat flux anomaly equals:

DH_d = c r_d A_d w_d DT_d (2)

where c is the heat capacity of water. Using (1) yields:

DH_d = - c r_u A_u w_u DT_d (3)

In the upwelling area the vertical distribution of the temperature anomaly is non-uniform and described by an advection-diffusion equation[5]:

c r_u A_u ¶DT_u/¶t = - ¶DH_u /¶z (4)

where[6]

DH_u = c r_u A_u( - k ¶DT_u /¶z + w_u DT_u) (5)

is the heat flux anomaly in the upwelling area and k a diffusion coefficient.

The following boundary conditions are assumed:

At the bottom of the mixed layer[7]:

z = - h_m : DT_d = DT_u = DT_m (6)

At the bottom of the ocean:

z = - h_d : DH_u + DH_d = 0 (7)

Using r_u » r_m , A_u » A_m , (3), (5) and (6) it follows that at z = - h_m:

DH_m = DH_u + DH_d = - c r_m A_m k [¶DT_u /¶z]_z=-hm (8)

The heat budget of the oceanic mixed-layer can be written as the sum of the flux, which enters at the bottom of the layer minus the flux, which leaves at the top:

C_m dDT_m/dt = DH_m - DH_s (9)

where

C_m = c r_m A_m h_m (10)

is the total heat capacity of the oceanic mixed layer and DH_s the heat flux into the atmosphere.

Note (EB): In order to judge the signs of variables, it is always handy to look at a simple example. Let us therefore take a heating up of the mixed layer because of IR down from the atmosphere. Ignore for a moment the DH_m from exchange with lower ocean layers. Then it is clear from (9) that DH_s must be negative, which agrees with the last statement above this note.

RADIATIVE FORCING

The flux (-DH_s ) going down, can be written as:

- DH_s = ADF - ADT_m/S_eq (11)

where DF is the radiative forcing per unit area[8] and A DT_m/ S_eq the increase in the outgoing radiative flux due to the temperature response of the climate system to this forcing. S_eq is the equilibrium climate sensitivity[9] [K/W/m²].

From (10) it follows that (9) can be written as

c r_m A_m h_m dDT_m/dt = - DH_s + DH_m (12)

TOTAL

From (8) - (12) it follows that:

c r_m A_m h_m dDT_m/dt = ADF_t - ADT_m/S_eq - c r_m A_m k [¶DT_u /¶z]_z=-hm (13)

Dividing (13) by c r_m A_m we obtain:

h_m dDT_m/dt = DF_t/c_w - DT_m/c_wS_eq - k [¶DT_u /¶z]_z=-hm (14)

where

c_w = c r_m A_m/A (15)

NUMERICAL PROCEDURE[10]

In order to solve the equations, the problem is made discrete. Below the mixed layer N ocean layers are taken. They are indicated by an index i, with temperatures in the middle represented by T_i. The thickness of a layer is taken as s_i. The subscripts u are omitted as well as the D. Of course they should be added in the final solution. We give the equations for a discrete number of levels and calculate the coupled differential equations with respect to time in Mathematica. One also may discretize the time and calculate everything in Excel, which was done by van Ulden in his original paper and used in our code ‘oceans-1’.

Below the mixed layer the heat budget for layer i follows from (4) and (5) and can be written as:

s_i ¶T_i/¶t = [w T - k ¶T /¶z]_i+1/2 - [w T - k ¶T /¶z]_i-1/2 (16)

where the subscript i+1/2 denotes the half level at the bottom of the layer and i-1/2 the half level at the top of the layer. The following boundary conditions hold:

- at the top of the first layer below the mixed layer (z_1/2 = -h_m):

T _1/2 = T_m (17a)

[¶T /¶z]_1/2 = (T_m - T₁)/( s₁ /2) (17b)

- at the bottom of the lowest layer ( z_N+1/2 = -h_d):

[w T - k ¶T /¶z]_N+1/2 = w T_d =w T_m (17c)

At the other half levels, values of T and ¶T /¶z are obtained by interpolation:

T_i+½=(s_i+1 T_i + s_i T _i+1)/( s_i + s_i+1) (18)

¶T /¶z _i+1/2 = (T_i – T _i+1)/(½ s_i + ½ s_i+1) (19)

Note that (17b)is somewhat different from (19) as it takes into account that the mixed layer has the same temperature everywhere. For completeness we also write (14) for the mixed layer in our simplified notation:

h_m dT_m/dt = F_t/c_w - T_m/c_wS_eq - k [¶T /¶z]_z=-hm (20)

UNITS

In climate studies the year is a more convenient time scale than the second. Therefore we express all our variables and parameters in years. For the parameters this is done by van Ulden in the next subsection. Eq. (20) shows us that we find the correct variables if we take 2 precautions:

- The forcing F_t is by convention expressed in Wm^-2 = J s^-1 m^-2. In our new units the forcing should be expressed in J yr ^-1 m^-2. Therefore F_t in Wm^-2 should be multiplied by the length of the year in seconds.

- the second term on the right in (20) contains S_eq, which is expressed in K/ Wm^-2 = K J^-1 s m^-2. When we change the s into yr we have to divide by the length l_yr of the year in seconds. Equivalently, the second term on the right of (20) should be multiplied by the length of the year l_yr. Concluding, the first two terms on the right in (20) should be multiplied by the length of the year l_yr. Then everything is OK. With the right parameters, (16) then is correct as well.

CHECKS

For a constant forcing F_t = 1 Wm^-2 all ocean layers should go to a temperature increase of S_eq. For the parameters given below this happens after a few 1000 years. Incidentally, one may check that this homogeneous temperature increase makes all derivatives zero and therefore gives a solution of eqs (16) and (20).

ACCURACIES

In the numerical procedure the number of layers is fixed. If one computes eqs (16) to (20) with Excel one has to choose a time step and define the forcing F(t) at each time point to be considered. If one has a strongly varying forcing as is the case with the description of volcanic eruptions, it is not straightforward in Excel to decrease time steps in solving the time dependent equations. Also it is not straightforward to define the resulting accuracies.

In using Mathematica as programming language the time steps in solving the equations are reduced automatically until the desired accuracy is achieved.

We find that for a constant forcing and a blockpulse, Excel with time steps of 0.25 years give a good approximation to the Mathematica results. For a strongly varying forcing Excel is not reliable. On the web (www.nat.vu.nl/envphysexp) we show Excel results, which for volcanic forcing only may be used ‘to guide the eye’.

Another aspect concerns the space steps to take. When one has got two adjoining layers of very different thickness, eq. (18) will put its weight on the thinnest layer. This may lead to erroneous results. When for example the thinner layer has the higher temperature and the advection points into the direction of that thinner layer, the advection term in eq. (16) will send too much heat from the colder layer into the warmer layer. This could even result in negative temperature anomalies in the thicker layer. The conclusion here is that one needs physical insight in deciding what steps in space one will take.

PARAMETERS

Numerical values given by Hoffert et al. (1980).

A_m/A = 0.65 r_m = 1030 kg/m³ c = 4000 J/kgK

c_w = c r_m A_m/A = 2.678 x10⁶ J/m³K

l_yr = 3.158 x10⁷ s

k = 2000 m²/y = 6.34 x10^-5 m²/s = 0.634 cm²/s

w_u = 4m/y = 1.27 x10^-7 m/s

Further we took h_m = 52.8 m , including a layer of 2.8 m of water representing the heat capacity of the troposphere-land system. And S_eq= 0.6 K/W m^-2.

Explanation of the simulation ocean-1 or ocean-2 (on the web)

We put in 4 standard radiative forcings, described below and the possibility of editing your own forcing.

Volcanic forcing

Explosive volcanic eruptions lead to the formation of stratospheric dust layers consisting of small aerosols. These dust layers remain in the stratosphere for 1 or 2 year, reflect sunlight and produce a negative radiative forcing which cools the earth surface. Van Ulden and van Dorland reconstructed the radiative forcing due to major volcanic eruptions from various observations and historic data. Van Ulden and van Dorland used volcanic eruptions where SO₂ was liberated, reducing the radiative forcing in the atmosphere. In the period from 1850 to 1999, one may recognize the eruption of the Krakatao in 1883 and of Pinatubo in 1991 . The negative forcing was calculated using historic data, and shown in the top graph. The resulting negative temperature pulses are shown in the bottom graph. One observes that the pulses propagate slowly to the lower layers and smoothen out as well.

Block pulse

A simple pulse in order to observe the behaviour of the water temperatures

Single pulse

A single Gaussian pulse

Predicting

When the radiative forcing increases like an exponential one may observe the subsequent temperatures

Editing forcing

The most simple forcing is the constant straight line (type in A=0 and B=1). This forcing should eventually lead to a temperature increase of 0.6 K. On the given time scale of 100 years this does not happen. One may however type in an own time scale as e.g. 10 000 years. Then it indeed converges. One also observes that at first the top and bottom layer heat up (because of the down-welling connecting both), while eventually all layers heat up in the order from top to bottom.

---

[1] The model described here is a summary by van Ulden of the equations used by A. P. van Ulden and R. van Dorland in their calculations (see Aad P. van Ulden and Rob van Dorland, Natural Variability of Global Mean Temperatures: Contributions from solar irradiance changes, volcanic eruptions and El Nino, Proc. 1^st Solar and Space Euroconference, ‘The Solar Cycle and Terrestrial Climate’, Santa Cruz de Tenerife, Tenerife, Spain, 25-29 September 2000, (ESA SP-463, December 2000). The footnotes were added by E. Boeker in order to clarify the connection between this paper and those used in ‘Environmental Physics’ of Egbert Boeker and Rienk van Grondelle (Wiley, 1999. This book is indicated by EP). The authors acknowledge Dr van Ulden’s permission to use his paper and his codes. The text after eq. (15) is added by E. Boeker in order to clarify the numerical procedure.

[2] Martin I. Hoffert, Andrew J. Callegari and Ching-Tzong Hsieh, The Role of Deep Sea Heat Storage in the Secular Response to Climatic Forcing, J Geophys. Res. 85, C11, November 20,1980, pp 6667-6679).

[3] Conservation of mass

[4] This means that on the timescales, which are considered here (years) the down-welling is described as a small column with the same temperature anomaly from the bottom of the top layer down to the bottom layer

[5] EP eq. (4.14) without heat production

[6] EP eq. (4.1) with at the right added an advection term, written down in EP eq. (5A.7) or EP eq. (5.4).

[7] At the bottom of the mixed layer the temperature anomaly of upwelling and down-welling are the same, as the temperature will change continuously at the interface.

[8] According to convention the radiative forcing is counted positive, when downwards. All other fluxes are counted positive if upwards.

[9] Assume for a moment that the anomalous heat flux at the tropopause is zero: DH_t=0. Then eq. (11) is precisely EP eq. (3.6). Consequently the equilibrium climate sensitivity S_eq is equal to the gain factor G of EP eq. (3.6). From EP eq. (3.17) it follows that the gain factor G is the temperature response after long times (equilibrium, say) for a radiative forcing equaling 1 W m^-2. This indeed is the way van Ulden et al have defined equilibrium climate sensitivity. This is more straightforward than the older IPCC definition, where the climate sensitivity is the response to doubling of CO₂ concentration. For the corresponding forcing will change all the time with improved models.

[10] This section on procedure is adapted by E. Boeker, the section on parameters is again completely by van Ulden.