CO2 data and trends

This post was inspired by a figure I saw recently, meant to highlight the critical state of atmospheric CO2 and what kinds of predictions we can make from the historical data. The original article, from the National Observer, can be found here, but what struck me as interesting was this figure:

atm c02 trends

NOAA Mauna Loa C02 data

We can get the monthly mean data from:

ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt

d <- read.table('ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt',
                na.strings = -99.99)
names(d) <- c('year', 'month', 'decimalYear', 'average', 'interpolated', 'trend', 'ndays')
d$time <- with(d, as.Date(paste(year, sprintf('%02d', month), '01', sep='-'), format='%Y-%m-%d', tz='UTC'))
time <- d$time
co2 <- d$interpolated

We can smooth out the seasonal variations with a moving average (i.e. boxcar) filter:

library(zoo)

Attaching package: 㤼㸱zoo㤼㸲

The following objects are masked from 㤼㸱package:base㤼㸲:

    as.Date, as.Date.numeric

co2sm <- rollmean(co2, 12, NA)
plot(time, co2, type='l')
lines(time, co2sm, col=2)

Fits and predictions

Now lets go through the excercise shown from the above figure to produce a series of linear fits based on the decadal values.

starts <- as.Date(c('1960-01-01', '1970-01-01', '1980-01-01', '1990-01-01',
                    '2000-01-01', '2010-01-01'))
plot(time, co2sm, type='l', ylab='C02', xlab='', lwd=3,
     xlim=c(min(time), as.Date('2050-01-01')),
     ylim=c(min(co2), 500))
abline(h=450, lwd=2, col=2)

i <- 1
for (s in starts) {
  II <- s <= time & time <= (s + 10*365)
  m <- lm(c ~ t, data=list(c=co2sm[II], t=time[II]))
  tp <- as.Date(seq(s, as.numeric(as.Date('2050-01-01')), 365))
  lines(tp, predict(m, newdata=list(t=tp)), col=i)
  i <- i + 1
}
legend('topleft', legend=format(starts, '%Y'), lwd=1, col=seq_along(starts),
       title='Start of decade')

But wait – fitting a bunch of straight lines to a curve seems like an odd way of making predictions. Why don’t we try to fit the underlying functional form and use that to predict into the future.

The CO2 time series looks to be exponential in nature, e.g.

\[ CO_2 = a + b e^{t/\tau} \]

Let’s try to fit that form using nls(), and compare to the linear fit from data from 2010 to present:

t <- as.numeric(time)
nm <- nls(c ~ a + b*exp(t/tau), data=list(c=co2sm, t=as.numeric(time)),
            start=list(a=320, b=10, tau=8000))
plot(time, co2sm, type='l', ylab='C02', xlab='', lwd=4,
     xlim=c(min(time), as.Date('2050-01-01')),
     ylim=c(min(co2), 500))
abline(h=450, lwd=2, col=2)

tt <- seq(min(t), 30000, 365)
lines(tt, predict(nm, list(t=tt)), col=2, lwd=2)
abline(m)

expdate <- tt[which.min(abs(predict(nm, list(t=tt))-450))]
abline(v=expdate, col=2)
mtext(format(as.Date(expdate)), at=expdate, adj=1, col=2)

lindate <- tt[which.min(abs(predict(m, list(t=as.Date(tt)))-450))]
abline(v=lindate)
mtext(format(as.Date(lindate)), at=lindate, adj=0)

Looks like the 450 ppm mark will hit three years sooner than the current decadal trend predicts.