Compute coefficients for anthropogenic CO$_2$ perturbation approach¶
James Orr
LSCE/IPSL
26 March 2018
Introduction and set up to use R¶
The perturbation approach developped in Sarmiento et al.,(1992) identifies a linear relation between the Ocean anthropogenic change in $p$CO$_2$ and its ratio with the corresponding ocean change in total carbon :
where $z_0$ and $z_1$ are quadratic functions of temperature, i.e : \begin{eqnarray} z_0 = a_0 + a_1 T + a_2 T^2 \hbox{,} \end{eqnarray} and \begin{eqnarray} z_1 = b_0 + b_1 T + b_2 T^2 \hbox{.} \end{eqnarray}
Below we demonstrate how the coefficients needed to calculate $z_0$ and $z_1$ (i.e : $a_i$ and $b_i$) have been calculated considering two different reference states (1765 and 1870), with one set of coefficients for each.
The seacarb package is used to make these calculations (using the Rmagic environment in a Jupyter notebook).
0. Set up to use R¶
Start by loading Rmagic, which is part of Rpy2 for iPython.
%reload_ext rpy2.ipython
Then launch R and install the needed R libraries (packages), if they are not already installed. After installing, comment out corresponding line as below.
%%R
# remove.packages('seacarb')
# install.packages('seacarb')
# install.packages('abind')
Continue in the same R session.
1. Initialise the basic information needed to compute the ocean carbonate chemistry.¶
In the next cell, chose between the Global conditions (compute the GLO approach coefficients) and Mediterranean Sea conditions (compute the MED coefficients). That is, uncomment the preferred choice, and comment out the other (precede lines by #).
%%R
# Load library to compute carbonate chemistry
library(seacarb)
# REFERENCE atm xCO2: Mole fraction (ppm) measured in dry air
# xCO2o = 280.0 #(Sarmiento et al., 1992)
xCO2o = 277.860 #(xCO2atm at time 1765.0 in record from Meinshausen et al. (2017, GMD, Supplement)
# xCO2o = 287.290 #(xCO2atm at time 1870.0 in combined record from M17 & J. Simeon's forcing for ORCA05)
## Global ocean:
Temp <- seq( -1.5, 30, by=0.5)
Salt = 35
TA = 2300e-6
# Maximum perturbation
delpco2MAX = 280.
## Mediterranean Sea:
# Temp <- seq(13, 30, by=1)
# Salt = 38
# TA = 2600e-6
# delpco2MAX = 280.
2. Correct for atmospheric humidity (100%) at equilibrium, see Sarmiento et al. (1992 - Eq. 9)¶
Compute initial field pCO$_{2,o}$, which varies with temperature (due to humidity correction)
%%R
tk <- Temp + 273.15 #Absolute temperature (K)
LNespa <- 20.1050 - 0.0097982*tk - 6163.10/tk
espa <- exp(LNespa)
pCO2o <- xCO2o * (1 - espa) #Convert from mole fraction (ppm) to partial pressure (uatm), assumes Patm = 1 atm
%%R
help(carb)
3. Compute corresponding (pre-industrial) DIC field with DIC0 = f(T)¶
%%R
carb0 <- carb(flag=24, var1=pCO2o, var2=TA, S=Salt, T=Temp, P=0, Patm=1, Pt=0, Sit=0,
k1k2="l", kf="dr", ks="d", pHscale="T", b="u74", gas="standard")
dic0 <- carb0$DIC
4. Increase ocean pCO$_2$ incrementally & compute corresponding surface DIC each time for all $T$'s¶
That is, do so by specifying $\delta$pCO$_2$ from 0 to 280 ppm and computing the corresponding change in surface DIC.
%%R
DIC <- numeric(0) ; pCO2 <- numeric(0)
delDIC <- numeric(0) ; delpCO2 <- numeric(0)
%%R
DIC <- numeric(0) ; pCO2 <- numeric(0)
delDIC <- numeric(0) ; delpCO2 <- numeric(0)
for (i in seq(1,delpco2MAX,1)) {
#print(c("delpCO2 =", i))
carbdel <- carb(flag=24, var1=pCO2o+i, var2=TA, S=Salt, T=Temp, P=0, Patm=1, Pt=0, Sit=0,
k1k2="l", kf="dr", ks="d", pHscale="T", b="u74", gas="standard")
dic <- carbdel$DIC
pco2 <- carbdel$pCO2
DIC <- rbind(DIC, dic)
pCO2 <- rbind(pCO2, pco2)
delDIC <- rbind(delDIC, dic - dic0)
delpCO2 <- rbind(delpCO2, pco2 - pCO2o)
}
5. Compute linear regression coefficients (slope and intercept) as a function of temperature¶
Compute $z_0$ and $z_1$, i.e., the intercept and slope relating the ratio $\frac{\delta pCO_2}{\delta C_T}$ and $\delta C_T$. Do so by a linear regression.
%%R
delDIC <- delDIC * 1e+6 #Convert from mol/kg to umol/kg
ratio <- delpCO2 / delDIC
Z0 <- numeric(0) ; Z1 <- numeric(0)
for (i in 1:length(Temp)) {
z0 <- lm(ratio[,i] ~ delpCO2[,i])$coeff[[1]] #Intercept
z1 <- lm(ratio[,i] ~ delpCO2[,i])$coeff[[2]] #Slope
Z0 <- append(Z0, z0)
Z1 <- append(Z1, z1)
}
6. Fit each linear regression coefficient to a quadratic equation in Temperature¶
Fit each linear regression coefficient to a quadratic equation in Temperature to deduce the $a_i$ and $b_i$ coefficients.
%%R
t <- Temp
t2 <- t^2
s0 <- summary(lm(Z0 ~ t + t2))
s1 <- summary(lm(Z1 ~ t + t2))
print(s0)
print(s1)
7. Retrieve coefficients from quadratic fits, and write them to output file¶
%%R
# 7) Retrieve coefficients from quadratic fits, and write them to output file
# Coefficients, associated errors, and R^2 for Z0 = a0 + a1*t + a2*t^2
a0 <- (lm(Z0 ~ t + t2))$coeff[[1]]
a1 <- (lm(Z0 ~ t + t2))$coeff[[2]]
a2 <- (lm(Z0 ~ t + t2))$coeff[[3]]
ae0 <- s0$coefficients[[1,2]] #Std. Error of a0
ae1 <- s0$coefficients[[2,2]] #Std. Error of a1
ae2 <- s0$coefficients[[3,2]] #Std. Error of a2
aR2 <- s0$adj.r.squared #Adjusted R^2
# Coefficients, associated errors, and R^2 for Z1 = a0 + a1*t + a2*t^2
b0 <- (lm(Z1 ~ t + t2))$coeff[[1]]
b1 <- (lm(Z1 ~ t + t2))$coeff[[2]]
b2 <- (lm(Z1 ~ t + t2))$coeff[[3]]
be0 <- s1$coefficients[[1,2]] #Std. Error of b0
be1 <- s1$coefficients[[2,2]] #Std. Error of b1
be2 <- s1$coefficients[[3,2]] #Std. Error of b2
bR2 <- s1$adj.r.squared #Adjusted R^2
The calculation is now completed. Below are the $a_i$ and $b_i$ coefficients, and their respective standard error ($ae_i$ and $be_i$):
%%R
## present the coefficients :
## coefficients ai, bi, and ci (i index varies from 0 to 9) used to compute z1, z2 and z3 (equation 13 of the manuscript):
part00 <- data.frame("i", "ai", "bi")
part0 <- data.frame("0", a0, b0)
part1 <- data.frame("1", a1, b1)
part2 <- data.frame("2", a2, b2)
partA0 <- data.frame( a0,'+/-',ae0, a1,'+/-',ae1, a2,'+/-',ae2)
partB0 <- data.frame( b0,'+/-',be0, b1,'+/-',be1, b2,'+/-',be2)
partA0
#cat(partB0)
%%R
partB0
8. Show results in a nice table (using Pandas in python).¶
#import numpy as np
#import pandas as pd
from rpy2.robjects import r, pandas2ri
pandas2ri.activate()
r.data('partAO')
ai = r['partA0']
ai
r.data('partBO')
bi = r['partB0']
bi
9. Assess error in the $\delta pCO_{2,o}$ calculated from the perturbation approach relative to "true" values from seacarb package¶
For a 200-ppm change in atm xCO2 (280 to 480 ppm), compute relative error in perturbation approach.
%%R
# 8) Verify Fitted vs. True results
x0 <- a0 + a1*t +a2*t2
x1 <- b0 + b1*t +b2*t2
# For a 200 ppm change of atmospheric xCO2 (280 to 480 ppm),
# corresponding delta DIC will vary (from changes in T)
dxCO2true = 200
xCO2true = 280 + dxCO2true
# Convert xCO2 to pCO2 (and calculate corresponding dpCO2)
pCO2true <- xCO2true * (1-espa)
dpCO2true <- pCO2true - pCO2o
# Compute "true" reference value from seacarb
carbd <- carb(flag=24, var1=pCO2true, var2=TA, S=Salt, T=Temp, P=0, Patm=1, Pt=0, Sit=0,
k1k2="l", kf="dr", ks="d", pHscale="T", b="u74", gas="standard")
ddic <- (carbd$DIC - dic0) * 1e+6
# Compute delta pCO2 from perturbation approach (using true ddic from seacarb)
dpCO2fit <- (x0 * ddic) / (1 - x1*ddic)
# Plot results
plot(x=Temp, y=100*(dpCO2fit - dpCO2true)/dpCO2true, col="blue")
# Rough conclusion (from above plot): relative errors always < +/-0.3% when delxCO2 < 200 (i.e., atm xCO2 ~ 560 ppm)
Figure 1: Relative error as a function of temperature for the $\delta$pCO$_2$ computed by the perturbation approach relative to reference values from seacarb.¶
Conclusion¶
Relative errors remain less than +/-0.3% across the observed temperature range of the Global Ocean.