A simple energy balance model of global warming

(relaxation to steady state response temperature)

General problem

One important feature of the anthropogenic climate change (global warming) is that the rise of temperature lags the rise of carbon dioxide concentration in the atmosphere. It is not a simple delay (time shift), but relationship between stock and flow (CO₂ concentration makes the flow of energy, and temperature is related to the stock of energy). Thus, it involves an integral in time.

The delay can schematically (not very realistically) be expressed in a zero-dimensional energy balance model, where the stock of energy (per unit area of the earth's surface) is equal to surface temperature T multiplied by constant heat capacity (per unit area) C.

If we tentatively assume that energy input to the climate system (per unit area per unit time) is F, the system can be written in the following equation.

C dT(t)/dt = F(t)

If integrated in time,

T(t) = T(0) + (1/C) �ç_t'=0^t F(t') dt' .

Steady-state response

We assume that the pre-industrial climate is in a steady state, with CO₂ concentration q_preind and temperature T_preind.

[Note: I just tentatively use "q" for CO₂ concentration. It does not mean specific humidity (concentration of water vapor) in this context.]

We assume that the steady-state response of temperature to CO₂ doubling is ƒÂT_ssr,double. [If CO₂ concentration would be 2 q_preind for an indefinitely long time, the average temperature would be T_preind + ƒÂT_ssr,double.

From experience of modeling studies, steady-state response is approximately proportional to logarithm of CO₂ concentration. Thus, it can be written like this.

T_ssr(t) = T_preind + ƒÂT_ssr,double �E{ (log(q(t)) - log(q_preind) }/ log(2)

(Either log-base-10 (common logarithm) or log-base-e (natural logarithm) can be used as far as the same base is used consistently.)

Simple formulation of energy input: relaxation towards steady-state

As for a simple example, we assume that F is proportional to the difference between the actual temperature and the steady-state response temperature.

F = k �E { T_ssr(q(t)) - T(t) }

If T_ssr is constant, T approaches T_ssr according to an exponential function of t (with a negative coefficient). In applied physics, this type of behavior is called "relaxation" or "damping".

Specification of a problem with some CO₂ concentration scenarios

The RCP (reference concentration pathways) scenarios have carbon dioxide concentration from year 2010 to 2100 by 10 year intervals. It is shown in a table whose value is taken from the Table A II.4.1 of IPCC AR5 (WG1). For convenience, CO₂ concentration at year 2000, known based on observations (Table A II.1.1a), is also included in our table.

• ar5_co2conc_rcp.txt
• (Supplementary info) ar5_co2conc_obs.txt

Pre-industrial level of CO₂ concentration q_preind is assumed to be 280 ppm. (According to IPCC AR5 Table A II.1.1a, CO₂ concentration took this value in year 1790.)

Pre-industrial level of temperature is not precisely known. IPCC AR5 Table A II.3 contains historical decadal mean global surface air temperature. According to "HadCRUT4_Median" (compiled by UKMO Hadley Centre and Univ. East Anglia's Climatic Research Unit), the temperature in 1850 was -0.320 K and that in 2000 was +0.453 K, relative to 1961-1990 average. We (tentatively) assume that T in 1850 as T_preind, and show all temperature values as deviation from this pre-industrial level. Temperature in 2000 is 0.773 K, in this convention.

• (Supplementary info) ar5_temp_global_hadcrut4.txt

Steady-state response to CO₂ doubling ƒÂT_ssr,double is (tentatively) assumed to be 3.0 K.

Heat capacity is assumed to be that of a homogeneous layer of water. Heat capacity per unit area C is the product of the depth of layer (tentatively assumed to be 500 m), density of water (1 g/cm³ = 1000 kg/m³), and specific heat capacity of water (approx. 4.2 J/(g K) = 4200 J/(kg K) ).

The coefficient of relaxation k cannot be determined by physics. I tentatively give the value 1.5 W/(m² K), based on empirical matching of results of test runs reproducing the history of the 20th century with observations.

As temperature is expressed in deviation from the pre-industrial level, the model equations are as follows:

C dT(t)/dt = k { T_ssr(q(t)) - T(t) }

T_ssr(q(t)) = ƒÂT_{ssr, double}�E{ log (q(t)) - log (q_preind) } / log(2)

And the differential equation is approximated by a finite difference formula like this.

T(t + ƒ¢t) = T(t) + ƒ¢t�Ek�E{ T_ssr(q(t)) - T(t) } / C

Programing in R

I have made an example program in R language.

R software is designed for statistical data analysis. It is also good for making graphs. It is not specifically designed for time step simulations, but it can be used.

The implementation software of R is available free from CRAN (https://cran.r-project.org ) or its mirror sites. There are packages for MS Windows, Mac, and Linux.

Download the following file of R script in the same directory as the above data. Open it by R, and run.

• co2warming.R
• Note: Do not try to browse the source program and save the text. Choose the file name link, show menu of actions (perhaps by clicking the right button of your mouse), and choose "save link" with a file name ending with ".R".

By one run, it makes one time series of temperature according to one of the RCP scenarios chosen by the variable "kscenario".

The part for output is unfinished. It just dumps values of array variables. (Also note that the value of "tssr" for the year 2110 is invalid.)

[This paragraph is added on 2019-09-01.] A version of R script which displays the results in a graph is added. (This version does not show numbers on the console or files, however.)

• co2warming_graph.R (revised on 2019-09-16)

[This part is added on 2019-09-16.] A version of R script which writes the results into a file is added.

• co2warming_write.R

As the scripts is run without modification, the results will be stored in a file like this.

• co2warming_rcp45_depth500_tssrdouble3_relax1_5.txt

Another R script which reads this file and makes the same graph as "co2warming_graph.R" is here.

• co2warming_read_graph.R [end of addition 2019-09-16.]

Formatted output will be added later (but I do not promise I finish by what time).

Programing in spreadsheets

The calculation can be done in spreadsheets of MS Excel or Libreoffice Calc. It is done by putting formulas of calculation in the cells. Users must care "absolute reference" and "relative reference" of cells. When a formula is copied, absolute reference (denoted with "$") does not change, but relative reference (denoted without "$") shifts.

(I do not recommend using spreadsheets if you can use R, Python or Ruby.)

Cells with "..." should be filled with data values from RCP (shown above). (Perhaps you can open another sheet and take data from file there, and then copy the values and paste onto the worksheet for simulation.)

"�«" and "�¨" means that the formula shown at the upper left corner should be copied to all cells in 4 columns and 11 rows.

Exercises

• Compare the results between concentration scenarios (Change the number given to "kscenario" in R).
• Try different values of steady state response (SSR) to CO₂ doubling ("tssrdouble" in R), e.g. 1.5 K or 4.5 K.
• Try different values of depth of water layer ("depth" in R), e.g. 4000 m (whole ocean) or 100 m (mixed layer near surface).
• If interested, try different values of the coefficient of relaxation ("relax" in R), which is the most uncertain.