The Azimuth Project
Experiments with carbon cycle box model



I want to create a couple of box models to get to understand the global carbon cycle fluxes and reservoirs better. In this first experiment I look at the perspective from the atmosphere mostly and also restrict it to the the biochemical part of the cycle, not the longest rock cycle. Also I neglect volcanic eruptions/flows which are at a rate of 0.06 Gt carbon per year. My intention is to give an introduction to both the carbon cycle and box modelling.


I read some background material on the global carbon cycle and one (see the reference section in Post, Peng etal) quotes another paper published in 1957 by the late Roger Revell and Hans Suess:

Planetary scale experiment in which mankind is returning to the atmosphere and oceans the concentrated organic carbon stored in sedimentary rocks over hundreds of millions of years.

That is almost 60 years since that was published. So let us get some updated data on the reservoirs and flows in the carbon cycle. From our own page on biomass we have an estimate on 560 Gt, so the estimate used by Rudy Slingerland and Lee Kump on 600 Gt of carbon in the biophere is pretty good and easier to remember!

One box models

Just some assumptions and notation which I will use. When you use a box model you are assuming that whatever the boxes represent is uniformly spread . FF will be used to represent flows of mass into or out of a box and mm is the mass of the box in question. As there is a lot of carbon I’ll use the customary unit Giga-tonne (Gt) which is one (US) billion tonnes or 10 910^9 tons and a tonne is 1000 kg. Flows are measured in Gt per year usually.

If you want to convert to CO 2CO_2 Gt, just multiply the Gt carbon quantity with 3.66. When I need to convert to volumetric parts per million (ppmv) you can use the fact that 2.1 Gt CC in the atmosphere is 1 ppm - See Experiments in carbon cycle with Sage or Carbon emissions for further explanations.

The other concepts I’ll introduce now is residence time which is roughly how long a flow stays in a particular box, which I’ll denote with τ\tau. Above we used kk as the loss or growth rate and it is defined as k=1τk = \frac{1}{\tau}. Steady state, is setting the setting the equations to 00 and solving them for values of mass. m˙\dot m denotes the time derivative of mm and m\vec{m} signifies the vector m.

The mass balance equation for one box models is:

m˙=sourcessinks \dot m = \sum{sources} -\sum{sinks}

Where the sources are what flows into the reservoir (box) and any internal sources (production and emission). The sinks is what flows out plus internal sinks (loss and deposition) in the reservoir. Here is a picture of this:

general one box model

Carbon growth rates in the atmosphere

Assuming mm is the amount of carbon in the atmosphere, we know from the Experiments in carbon cycle with Sage measurements that in 2008 the annual growth is about 1.8 CO 2CO_2 ppm (v) or 3.9 Gt carbon/year. You can see this in the lower diagram on that page looking at the yellow moving average for 2008.

m co2˙=sourcessinks \dot{m_co_2}=\sum\sources-\sum{sinks}

For the atmosphere as a one box model, the major contributor to the emission sources is fossil fuel emission, see Carbon emissions and for 2008 its growth rate was around 8.7 Gt CC or 3.5 ppm CO 2CO_2 per year . Compare this to the volcanic flow which is less than 1% of this - so it is reasonable to neglect it here.

These values are our own estimates which uses the Mauna Loa measurements as a time series for CO 2CO_2 pppm and I also have verified them from the book by David Neelin in the references below. Here is how CO 2CO_2 is distributed in the atmosphere:

co2 distribuiton measured by HIPPO

Images from the HIPPO project who aims towards:

The HIAPER Pole-to-Pole Observations (HIPPO) project is investigating the carbon Cycle and greenhouse gases throughout various altitudes of the western hemisphere through the annual cycle.

So the main objective of HIPPO is to add a vertical dimension to all GHG measurements.

Forest burning is a source of CO 2CO_2 and for 2008 it was estimated to 1.6 Gt/y or 0.7 ppm CO 2CO_2, so

sources=3.5+0.7=4.2ppm/y \sum\sources = 3.5 + 0.7 = 4.2 ppm/y

Using the mass balance equation above:

sinks=4.21.8=2.4ppm/y \sum{sinks} = 4.2 - 1.8 = 2.4 ppm/y

So we have less than half of the carbon in the in the atmosphere than the fossil fuel emission ! So what is going on here?

Let us look at the Carbon cycle diagrams we see that the Ocean is a candidate to look at. So let us do that. We got 92 Gt carbon per year flowing into the oceans and 90 Gt/year so we have about 2 Gt annual uptake in the oceans (and slightly more in 2008 according to David Neeling; at 2.2 Gt). So we have an oceanic sink of 2.2 Gt/year carbon or 1.1 ppm/year of CO 2CO_2 which explains most of this, but not all. so let me explore some other sinks and sources in the CO 2CO_2 global box model.

The biomass sink is 2.5 Gt/y or 1.2 ppm CO 2CO_2. So currently the land and ocean are sinks for about half of the carbon emissions from fossil fuel. There is concern about this as both the ocean and land sink appears to shrink and its important to keep track of how fast it shrinks. This is from Le Quéré et al., Trends in the sources and sinks of carbon dioxide in the reference section of carbon cycle. So below negative values are sinks and positive sources and the first green curve shows land and the second blue ocean sinks:

shrinking sinks for co2

So when I sum all source and sinks it doesn’t add up to zero and this is typical, because of how good or bad the model reflects actual measurements. One example is the soil in the land sink - should I use one meter or two of soil or something in between to calculate that sink ? Recently the Global Carbon Project published its report for 2000 to 2010 and here it continues with the same trend upto 2010:

global carbon budget 2010

Monitoring of the shrinking sinks can also be done by observing the Airborne fraction.

Radiocarbon growth

I will use one example from Slingerland which also has an analytic solution. Its a linear model with an initially empty atmosphere which simulates the 13C^{13}C growth:

linear carbon two box model

Radio carbon is usually used for finding out the age of things and in earth and climate science it is mainly used to distinguish and also trace different Greenhouse gases. Several isotopes are used including isotopes of carbon and oxygen. See our page on isotope geochemistry and its application in isotopic signatures and environmental forensics.

More boxes

It is simple to state if I use vector notation:

m˙=sourcessinks \dot \vec m = \sum{sources} -\sum{sinks}

It looks very similar using vector notation! But let us look at a 2-box model to show how we go from two explicit equations towards vectors.

Here is the first naive linear box model from Slingerland and Kump, where M 1M_1 is the carbon content in the atmosphere and M 2M_2 is the biomass on Earth. The box model is described by:

M 1˙=k 2M 2k 1M 1 \dot{M_1} = k_2M_2 -k_1M_1
M 2˙=k 1M 1k 2M 2 \dot{M_2} = k_1M_1 -k_2M_2

Where k 1k_1 is the rate constant for photosynthesis and k 2k_2 for respiration and decay. So the assumption here is that the rates are proportional to their masses. For typical values; M 1=800M_1=800 Gt and M 2=600M_2=600 Gt and fluxes k 2M 2=k 1M1=60k_2M_2 = k_1M1=60 Gt/year we can calculate at steady state, (setting the differential equations to 00). This gives us k 1=0.075k_1=0.075 and k 2=0.1k_2=0.1 per year. The inverse of the k-values is known as the residence time or also removal time; how long the flux stays in atmosphere respective biomass.

I made a plot The that shows tentatively what happens when we suddenly move 300 Gt to the atmosphere, like an asteroid crashing into Earth:

linear carbon two box model

So if we introduce vector notation I can write M 1,M 2\M_1,M_2 as M=(M 1,M 2)\vec M = (M_1,M_2). The response time (residence time for the whole system) is about 6 years, a lot shorter than the individual boxes residence time for 10 years for the biosphere and more than 13 for the atmosphere.

Non-linearity in the model

Here is the radiocarbon example by Rudy and Lee, with 11-year sunspot cycles included in the production PP. I doubled the amplitude to see the cycles better Here we see the start (or transient) to the left and after steady-state to the right:

transient of non linear radio carbon one box model

add stationary pic above

Here is the 2-box example from Kump, but this time the photosynthesis rate is proportional to both the mass of carbon in the atmosphere and the biosphere:

response of non linear carbon two box model

So the response time of 23 years is approximately 4 times longer than in the linear case of 6 years.

In the second experiment I will scrutinize the effects on the oceans in Experiments with ocean carbon box model.

Sage code

It will just be here temporarily if you want to experiment yourself, with the program in progress. It is code to solve the ODE-system which is a bit of overkill for the linear part, but I needed it for the non-linear examples. This is just for the 2-box mode:

# Dynamical Systems (MMEDS)
# Slingerland et. al PUP 2011

var('t k1 k2')

M1 = function('M1', t)
M2 = function('M2', t)

# rate of removal ~ mass. there is a typo in eq 3.17 and onwards to 3.21 in MMEDS,
# but it does not have any effect on the model.

# photosynthesis rate
F12 = k1*M1

# decay and respiration rate
F21 = k2*M2

# initial transfer of 300 Gt carbon to the atmosphere

de1 = diff(M1,t) == F21 - F12 
de2 = diff(M2,t) == F12 - F21
sol = desolve_system([de1, de2], vars= [M1,M2],ivar=t,ics = [0,1100,300])

# extract the symbolic solutions
solm1,solm2 = sol[0].rhs(),sol[1].rhs()

# linear carbon cycle plot
m1p = plot(solm1.substitute({M1(0):800,M2(0):600,k1:0.075,k2:0.1}),(t,0,25),color="purple", legend_label='$M_1$')
m2p = plot(solm2.substitute({M1(0):800,M2(0):600,k1:0.075,k2:0.1}),(t,0,25), color="black", legend_label='$M_2$')
(m1p + m2p).show()


shows several interesting animations and here is an animation of the global carbon dioxide 2002-2008 as measiured by the satelite AIRS