Routine to compute coefficients for anthropogenic CO2 perturbation approach (GLO or MED approach)

The perturbation approach developped in Sarmiento et al.,(1992) identify a linear relation between the Ocean anthropogenic change in pCO2 and its ratio with the corresponding ocean change in total carbon :

δpCO2oδCT=z0+z1δpCO2o.
where z0 and z1 are quadratic functions of temperature, i.e :
z0=a0+a1T+a2T2,
and
z1=b0+b1T+b2T2.

Below we demonstrate how the coefficients needed to calculate z0 and z1 (i.e : ai and bi) have been calculated.

Start by loading Rmagic, which is part of Rpy2 for iPython.

In [1]:
%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.

In [2]:
%%R
#  remove.packages('seacarb')
#  install.packages('seacarb')
#  install.packages('abind')
 
NULL

Continue in the same R session.

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

In [3]:
%%R
# Load library to compute carbonate chemistry
  library(seacarb)
 
#Mole fraction (ppm) measured in dry air
 xCO2o = 280.0  
 
## Global ocean:
#  Temp <- seq( 0, 30, by=1)
#  Salt = 35
#  TA = 2300e-6
#  delpco2MAX = 200. 
 
## Mediterranean Sea:
  Temp <- seq(13, 30, by=1)
  Salt = 38
  TA = 2600e-6
  delpco2MAX = 280.
 

Correct for atmospheric humidity (100%) at equilibrium, see Sarmiento et al. (1992 - Eq. 9) Compute initial field pCO2,o, which varies with temperature (due to humidity correction)

In [4]:
%%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
 

Compute corresponding (pre-industrial) DIC field with DIC0 = f(T)

In [5]:
%%R
carb0 <- carb(flag=24, var1=pCO2o, var2=TA, S=Salt, T=Temp, P=0, Pt=0, Sit=0, k1k2="l", kf="dr", ks="d", pHscale="T")
dic0 <- carb0$DIC
 

Increasing DIC incrementally and compute pCO2 for each increment until we get δpCO2 from 0 to 280 ppm and the corresponding change in surface DIC.

In [6]:
%%R
DIC  <- numeric(0) ; pCO2 <- numeric(0)
delDIC  <- numeric(0) ; delpCO2 <- numeric(0)
 
In [7]:
%%R
DIC  <- numeric(0) ; pCO2 <- numeric(0)
delDIC  <- numeric(0) ; delpCO2 <- numeric(0)
 
  for (i in seq(5,delpco2MAX,5)) {
    #print(c("delpCO2 =", i))
    carbdel <- carb(flag=24, var1=pCO2o+i, var2=TA, S=Salt, T=Temp, P=0, Pt=0, Sit=0, k1k2="l", kf="dr", ks="d", pHscale="T")
    dic <- carbdel$DIC
    pco2 <- carbdel$pCO2
    DIC  <- rbind(DIC, dic)
    pCO2 <- rbind(pCO2, pco2)
    delDIC  <- rbind(delDIC, dic - dic0)
    delpCO2 <- rbind(delpCO2, pco2 - pCO2o)
   }
 

Compute z0 and z1, i.e., the intercept and slope relating the ratio δpCO2δCT and δCT. Do so by a linear regression.

In [8]:
%%R
# 5) Compute linear regression coefficients (slope and intercept) as a function of temperature
  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)
  }
 

Fit each linear regression coefficient to a quadratic equation in Temperature to deduce the ai and bi coefficients.

In [9]:
%%R
# 6) Fit each linear regression coefficient to a quadratic equation in Temperature
    t <- Temp
    t2 <- t^2
    s0 <- summary(lm(Z0 ~ t + t2))
    s1 <- summary(lm(Z1 ~ t + t2))
    print(s0)
    print(s1)
 
Call:
lm(formula = Z0 ~ t + t2)

Residuals:
       Min         1Q     Median         3Q        Max 
-9.030e-04 -4.207e-04  1.154e-05  4.051e-04  9.382e-04 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.419e+00  2.345e-03  605.01   <2e-16 ***
t           -2.135e-02  2.268e-04  -94.13   <2e-16 ***
t2           2.416e-04  5.244e-06   46.07   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0005331 on 15 degrees of freedom
Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999 
F-statistic: 1.034e+05 on 2 and 15 DF,  p-value: < 2.2e-16


Call:
lm(formula = Z1 ~ t + t2)

Residuals:
       Min         1Q     Median         3Q        Max 
-9.253e-07 -3.428e-07 -1.964e-08  3.956e-07  7.532e-07 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.301e-03  2.188e-06  1508.8   <2e-16 ***
t           -6.663e-05  2.116e-07  -314.9   <2e-16 ***
t2           6.521e-07  4.892e-09   133.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.974e-07 on 15 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 1.467e+06 on 2 and 15 DF,  p-value: < 2.2e-16

In [10]:
%%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 ai and bi coefficients, and their respective standard error (aei and bei):

In [11]:
%%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
 
 
        a0 X.....         ae0          a1 X......1          ae1           a2
1 1.418634    +/- 0.002344821 -0.02134721      +/- 0.0002267815 0.0002415878
  X......2          ae2
1      +/- 5.243825e-06

Show results in a nice table (using Pandas in python).

In [12]:
import numpy as np
import pandas as pd
import pandas.rpy.common as com
 
In [13]:
ai=com.load_data('partA0')
bi=com.load_data('partB0')
 
def CreateDataSet(Number=1): 
    Output = []
    Output.extend(zip(ai, bi))
    return Output
 
 
In [14]:
#print ai.shape
ai
 
Out[14]:
a0 X..... ae0 a1 X......1 ae1 a2 X......2 ae2
1 1.418634 +/- 0.002345 -0.021347 +/- 0.000227 0.000242 +/- 0.000005
In [15]:
bi
 
Out[15]:
b0 X..... be0 b1 X......1 be1 b2 X......2 be2
1 0.003301 +/- 0.000002 -0.000067 +/- 0 0.000001 +/- 4.892404e-09

Assess the error in the δpCO2,o calculated from the perturbation approach relative to "true" values from the seacarb carbonate chemistry package.

For a doubling of atmospheric xCO2 (280 to 560 ppm), compute relative error in perturbation approach.

In [16]:
%%R
# 8) Verify on Fitted vs. True results
    x0 <- a0 + a1*t +a2*t2
    x1 <- b0 + b1*t +b2*t2
#   For a doubling of atmospheric xCO2 (280 to 560 ppm), corresponding delta DIC will vary (due to changes in temperature)
    dxCO2true = 280
    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, Pt=0, Sit=0, k1k2="l", kf="dr", ks="d", pHscale="T")
    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.5% when delxCO2 < 280 (i.e., atm xCO2 ~ 560 ppm)
 
Figure 1: Relative error as a function of temperature for the δpCO2 computed by the perturbation approach relative to reference values from seacarb. Relative errors remain less than 0.5% across the observed temperature range of the Mediterranean Sea.
In [ ]: