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_{2} 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' .

We assume that the pre-industrial climate is in a steady state,
with CO_{2} concentration *q*_{preind} and temperature *T*_{preind}.

[Note: I just tentatively use "q" for CO_{2} 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_{2} doubling
is *ƒÂT*_{ssr,double}.
[If CO_{2} 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_{2} 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.)

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".

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_{2} 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_{2} concentration *q*_{preind} is assumed to be 280 ppm.
(According to IPCC AR5 Table A II.1.1a, CO_{2} 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_{2} 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^{3} = 1000 kg/m^{3}), 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^{2} 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) + ƒ¢tEkE{ T

_{ssr}(q(t)) - T(t) } / C

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.

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

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).

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.)

A | B | C | D | E | F | G | H | I | J | K | L | M | N | O | P | Q | R | |

1 | A simple energy balance model of global warming (relaxation to steady state response tempearture) 2019-07-16 Masuda | |||||||||||||||||

2 | constants | year | CO2 conc. [ppm] | Steady state response [K] | Temp. (base year 1850) [K] | |||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

3 | CO2 pre-industrial [ppm] | 280 | RCP2.6 | RCP4.5 | RCP6.0 | RCP8.5 | RCP2.6 | RCP4.5 | RCP6.0 | RCP8.5 | RCP2.6 | RCP4.5 | RCP6.0 | RCP8.5 | ||||

4 | SSR to CO2 doubling [K] | 3 | 2000 | 368.7 | =$E4 | =$E4 | =$E4 | =$B$4*(LOG(E4)-LOG($B$3))/LOG(2) | 0.773 | =$O4 | =$O4 | =$O4 | ||||||

5 | 2010 | 389.3 | 389.1 | 389.1 | 389.3 | « | ¨ | ¨ | ¨ | =O4+(J4-O4)*$B$6*$B$16/$B$11 | ||||||||

6 | relaxation coeff. [W/(m2 K)] | 1.5 | 2020 | 412.1 | 411.1 | 409.4 | 415.8 | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

7 | 2030 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||||

8 | depth of water [m] | 500 | 2040 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

9 | density of water [kg/m3] | 1000 | 2050 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

10 | specific heat of water [J/(kg K)] | 4200 | 2060 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

11 | heat capacity / area [J/(m2 K)] | =$B$8*$B$9*$B$10 | 2070 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

12 | 2080 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||||

13 | time step [year] | 10 | 2090 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

14 | day / year | 365.25 | 2100 | ... | ... | ... | ... | « | ¨ | ¨ | ¨ | « | ¨ | ¨ | ¨ | |||

15 | second / day | 86400 | 2110 | « | ¨ | ¨ | ¨ | |||||||||||

16 | time step [s] | =$B$13*$B$14*$B$15 |

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.

- Compare the results between concentration scenarios (Change the number given to "kscenario" in R).
- Try different values of steady state response (SSR) to CO
_{2}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.

2019-Jul-16, 2019-Jul-17, 2019-Sep-01, 2019-Sep-16

MASUDA Kooiti