Routine to compute coefficients for anthropogenic CO2 perturbation approach (VAR approach)

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

δpCO2,o=0+z1δCT+z2δC2T+z3δC3T.
where z0 and z1 are quadratic functions of temperature, i.e :
z1=a0+a1T+a2S+a3T2+a4S2+a5T3+a6S3+a7TS+a8T2S+a9TS2,
z2=b0+b1T+b2S+b3T2+b4S2+b5T3+b6S3+b7TS+b8T2S+b9TS2,
and
z3=c0+c1T+c2S+c3T2+c4S2+c5T3+c6S3+c7TS+c8T2S+c9TS2.

In the following, we will show how the coefficients needed to calculate z1, z2 and z3 (i.e : ai, bi and ci) have been calculated.

So let's use R to compute the partials. 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. You may need to install the Hmisc library with same command in an R session on X terminal. After installing, comment out corresponding line as below.

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

Continue the same R session.

First, we initialise the basic informations that are needed to compute the ocean carbonate chemistry. --> Here we compute the Mediterranean sea conditions of the VAR experiment (with varying temperature and salinity).

In [3]:
%%R
# Load libraries to compute carbonate chemistry and make nice table (in LaTeX)
  library(Hmisc)
  library(seacarb)
  #library(Hmisc)
 
#Mole fraction (ppm) measured in dry air
 xCO2o = 280.0  
 
# Mediterranean Sea:
  delpco2MAX = 280.
  salt <- seq(33.5, 40.0, by=0.5)
  temp <- seq(13, 30, by=1)
  hydro <- expand.grid(temp, salt)
  Temp <- hydro$Var1
  Salt <- hydro$Var2
  TA <- Salt
  TA <- (73.7*Salt - 285.7)*1e-6     #Schneider et al. (2007)
  TA2 <- (93.996*Salt - 1038.1)*1e-6 #Copin et Begovic (2002)
 
 
Le chargement a nécessité le package : grid
Le chargement a nécessité le package : lattice
Le chargement a nécessité le package : survival
Le chargement a nécessité le package : splines
Le chargement a nécessité le package : Formula

Attachement du package : ‘Hmisc’

L'objet suivant est masqué from ‘package:base’:

    format.pval, round.POSIXt, trunc.POSIXt, units

now we 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
 
 
 

And here is computed the corresponding (pre-industrial) DIC field with DIC0 = f(T,S)

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
 
 

We continu, increasing DIC incrementally, and computing pCO2 for each increment, until we get δpCO2 from 0 to 280 ppm and the corresponding change in surface dissolved inorganic carbon.

In [6]:
%%R
DIC  <- numeric(0) ; pCO2 <- numeric(0)
delDIC  <- numeric(0) ; delpCO2 <- numeric(0)
Tc  <- numeric(0) ; Sc <- numeric(0)
 
delpco2RANGE <- c(0.1, 1, seq(5, delpco2MAX, 5))
delpco2RANGE <- c(seq(5, delpco2MAX, 5))
  for (i in delpco2RANGE) {
    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)
    tc <- carbdel$T
    sc <- carbdel$S
    Tc  <- rbind(Tc, tc)
    Sc  <- rbind(Sc, sc)
   }
 
[1] "delpCO2 =" "5"        
[1] "delpCO2 =" "10"       
[1] "delpCO2 =" "15"       
[1] "delpCO2 =" "20"       
[1] "delpCO2 =" "25"       
[1] "delpCO2 =" "30"       
[1] "delpCO2 =" "35"       
[1] "delpCO2 =" "40"       
[1] "delpCO2 =" "45"       
[1] "delpCO2 =" "50"       
[1] "delpCO2 =" "55"       
[1] "delpCO2 =" "60"       
[1] "delpCO2 =" "65"       
[1] "delpCO2 =" "70"       
[1] "delpCO2 =" "75"       
[1] "delpCO2 =" "80"       
[1] "delpCO2 =" "85"       
[1] "delpCO2 =" "90"       
[1] "delpCO2 =" "95"       
[1] "delpCO2 =" "100"      
[1] "delpCO2 =" "105"      
[1] "delpCO2 =" "110"      
[1] "delpCO2 =" "115"      
[1] "delpCO2 =" "120"      
[1] "delpCO2 =" "125"      
[1] "delpCO2 =" "130"      
[1] "delpCO2 =" "135"      
[1] "delpCO2 =" "140"      
[1] "delpCO2 =" "145"      
[1] "delpCO2 =" "150"      
[1] "delpCO2 =" "155"      
[1] "delpCO2 =" "160"      
[1] "delpCO2 =" "165"      
[1] "delpCO2 =" "170"      
[1] "delpCO2 =" "175"      
[1] "delpCO2 =" "180"      
[1] "delpCO2 =" "185"      
[1] "delpCO2 =" "190"      
[1] "delpCO2 =" "195"      
[1] "delpCO2 =" "200"      
[1] "delpCO2 =" "205"      
[1] "delpCO2 =" "210"      
[1] "delpCO2 =" "215"      
[1] "delpCO2 =" "220"      
[1] "delpCO2 =" "225"      
[1] "delpCO2 =" "230"      
[1] "delpCO2 =" "235"      
[1] "delpCO2 =" "240"      
[1] "delpCO2 =" "245"      
[1] "delpCO2 =" "250"      
[1] "delpCO2 =" "255"      
[1] "delpCO2 =" "260"      
[1] "delpCO2 =" "265"      
[1] "delpCO2 =" "270"      
[1] "delpCO2 =" "275"      
[1] "delpCO2 =" "280"      

Now we compute z1, z2 and z3, respectively the intercept, the slope, and the second term coefficient of the cubic relation between Mediterranean δpCO2 and δCT, from a polynomial regression.

In [7]:
%%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
 
  Z1  <- numeric(0) ; Z2  <- numeric(0) ;  Z3  <- numeric(0)
  nsalt = length(salt)
  ntemp = length(temp)
# for (it in 1:length(Temp)) {
  for (is in 1:nsalt) {
    for (it in 1:ntemp) {
      n <- (is - 1) * ntemp + it
      dp <- delpCO2[,n]
      ddic <- delDIC[,n]
      ddic2 <- ddic*ddic
      ddic3 <- ddic*ddic2
#     f <- lm(ratio[,n] ~ dp)         #Fit to ratio = a0 + a1*dp
#     f2 <- lm(ratio[,n] ~ dp + dp2)  #Fit to ratio = a0 + a1*dp + a2*dp2
#
#     Switch to polynomial of dpCO2 = f(dDIC)
      f <- lm(dp ~ 0 + ddic + ddic2 + ddic3)  #Fit to dpCO2 = 0 + a1*ddic + a2*ddic^2 + a3*ddic^3
      z1 <- f$coeff[[1]]  #ddic
      z2 <- f$coeff[[2]]  #Slope
      z3 <- f$coeff[[3]]  #coefficient on dp2 term
      Z1 <- append(Z1, z1)
      Z2 <- append(Z2, z2)
      Z3 <- append(Z3, z3)
    }
  }
 
 
 

Then we fit each coefficient to a cubic equation in Temperature and Salinity, so we can deduce the ai, bi and ci coefficients.

In [8]:
%%R
# 6) Fit each linear regression coefficient to a quadratic equation in Temperature
# 6) Fit each linear regression coefficient to multivariate quadratic equations in T and S
    t <- Temp
    s <- Salt
    t2 <- t^2
    s2 <- s^2
    s3 <- s2 * s
    t3 <- t2 * t
    ts <- t*s
    t2s <- t2*s
    ts2 <- t*s2
#   s4 <- s^4
#   sinv <- 1/s
    fitZ1 <- lm(Z1 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)
    fitZ2 <- lm(Z2 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)
    fitZ3 <- lm(Z3 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)
    sum1 <- summary(fitZ1)
    sum2 <- summary(fitZ2)
    sum3 <- summary(fitZ3)
    print(sum1)
    print(sum2)
    print(sum3)
 
 
 
Call:
lm(formula = Z1 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)

Residuals:
       Min         1Q     Median         3Q        Max 
-1.294e-04 -4.748e-05 -1.310e-06  4.645e-05  3.535e-04 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.165e+01  3.567e-02   326.5   <2e-16 ***
t           -1.591e-01  3.545e-04  -449.0   <2e-16 ***
s           -5.569e-01  2.888e-03  -192.8   <2e-16 ***
t2           1.088e-03  4.285e-06   253.9   <2e-16 ***
s2           1.016e-02  7.820e-05   129.9   <2e-16 ***
t3          -4.141e-06  4.102e-08  -101.0   <2e-16 ***
s3          -6.603e-05  7.077e-07   -93.3   <2e-16 ***
ts           5.218e-03  1.783e-05   292.7   <2e-16 ***
t2s         -1.509e-05  9.159e-08  -164.8   <2e-16 ***
ts2         -4.670e-05  2.365e-07  -197.5   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.022e-05 on 242 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 8.104e+07 on 9 and 242 DF,  p-value: < 2.2e-16


Call:
lm(formula = Z2 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)

Residuals:
       Min         1Q     Median         3Q        Max 
-4.929e-06 -2.358e-06  3.820e-08  1.774e-06  2.072e-05 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.115e-02  1.637e-03   55.69   <2e-16 ***
t           -1.806e-03  1.626e-05 -111.08   <2e-16 ***
s           -5.131e-03  1.325e-04  -38.72   <2e-16 ***
t2           1.570e-05  1.966e-07   79.86   <2e-16 ***
s2           1.005e-04  3.588e-06   28.02   <2e-16 ***
t3          -5.649e-08  1.882e-09  -30.01   <2e-16 ***
s3          -6.761e-07  3.247e-08  -20.82   <2e-16 ***
ts           6.670e-05  8.178e-07   81.56   <2e-16 ***
t2s         -2.755e-07  4.202e-09  -65.57   <2e-16 ***
ts2         -6.373e-07  1.085e-08  -58.74   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.222e-06 on 242 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 6.977e+05 on 9 and 242 DF,  p-value: < 2.2e-16


Call:
lm(formula = Z3 ~ t + s + t2 + s2 + t3 + s3 + ts + t2s + ts2)

Residuals:
       Min         1Q     Median         3Q        Max 
-5.466e-08 -1.985e-08  3.300e-10  1.677e-08  1.263e-07 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.315e-03  1.332e-05   98.72   <2e-16 ***
t           -2.357e-05  1.324e-07 -178.06   <2e-16 ***
s           -7.879e-05  1.079e-06  -73.04   <2e-16 ***
t2           1.520e-07  1.601e-09   94.97   <2e-16 ***
s2           1.613e-06  2.921e-08   55.23   <2e-16 ***
t3          -3.438e-10  1.532e-11  -22.44   <2e-16 ***
s3          -1.122e-08  2.643e-10  -42.45   <2e-16 ***
ts           9.291e-07  6.658e-09  139.54   <2e-16 ***
t2s         -2.937e-09  3.421e-11  -85.86   <2e-16 ***
ts2         -9.391e-09  8.833e-11 -106.32   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.623e-08 on 242 degrees of freedom
Multiple R-squared:      1,	Adjusted R-squared:      1 
F-statistic: 1.403e+06 on 9 and 242 DF,  p-value: < 2.2e-16

In [9]:
%%R
# 7) Retrieve coefficients from quadratic fits, and write them to output file
 
#   Z1: Coefficients, associated errors, and R^2 for Z1 = a0 + a1*t + a2*t^2
    a0 <- (fitZ1)$coeff[[1]]
    a1 <- (fitZ1)$coeff[[2]]
    a2 <- (fitZ1)$coeff[[3]]
    a3 <- (fitZ1)$coeff[[4]]
    a4 <- (fitZ1)$coeff[[5]]
    a5 <- (fitZ1)$coeff[[6]]
    a6 <- (fitZ1)$coeff[[7]]
    a7 <- (fitZ1)$coeff[[8]]
    a8 <- (fitZ1)$coeff[[9]]
    a9 <- (fitZ1)$coeff[[10]]
    ae0 <- sum1$coefficients[[1,2]]      #Std. Error of a0
    ae1 <- sum1$coefficients[[2,2]]      #Std. Error of a1
    ae2 <- sum1$coefficients[[3,2]]      #Std. Error of a2
    ae3 <- sum1$coefficients[[4,2]]      #Std. Error of a3
    ae4 <- sum1$coefficients[[5,2]]      #Std. Error of a4
    ae5 <- sum1$coefficients[[6,2]]      #Std. Error of a5
    ae6 <- sum1$coefficients[[7,2]]      #Std. Error of a6
    ae7 <- sum1$coefficients[[8,2]]      #Std. Error of a7
    ae8 <- sum1$coefficients[[9,2]]      #Std. Error of a8
    ae9 <- sum1$coefficients[[10,2]]     #Std. Error of a9
    aR2 <- sum1$adj.r.squared            #Adjusted R^2
 
#   Z2: Coefficients, associated errors, and R^2 for Z2 = a0 + a1*t + a2*t^2
    b0 <- (fitZ2)$coeff[[1]]
    b1 <- (fitZ2)$coeff[[2]]
    b2 <- (fitZ2)$coeff[[3]]
    b3 <- (fitZ2)$coeff[[4]]
    b4 <- (fitZ2)$coeff[[5]]
    b5 <- (fitZ2)$coeff[[6]]
    b6 <- (fitZ2)$coeff[[7]]
    b7 <- (fitZ2)$coeff[[8]]
    b8 <- (fitZ2)$coeff[[9]]
    b9 <- (fitZ2)$coeff[[10]]
    be0 <- sum2$coefficients[[1,2]]      #Std. Error of b0
    be1 <- sum2$coefficients[[2,2]]      #Std. Error of b1
    be2 <- sum2$coefficients[[3,2]]      #Std. Error of b2
    be3 <- sum2$coefficients[[4,2]]      #Std. Error of b3
    be4 <- sum2$coefficients[[5,2]]      #Std. Error of b4
    be5 <- sum2$coefficients[[6,2]]      #Std. Error of b5
    be6 <- sum2$coefficients[[7,2]]      #Std. Error of b6
    be7 <- sum2$coefficients[[8,2]]      #Std. Error of b7
    be8 <- sum2$coefficients[[9,2]]      #Std. Error of b8
    be9 <- sum2$coefficients[[10,2]]     #Std. Error of b9
    bR2 <- sum2$adj.r.squared            #Adjusted R^2
 
#   Z3: Coefficients, associated errors, and R^2 for Z3 = a0 + a1*t + a2*t^2
    c0 <- (fitZ3)$coeff[[1]]
    c1 <- (fitZ3)$coeff[[2]]
    c2 <- (fitZ3)$coeff[[3]]
    c3 <- (fitZ3)$coeff[[4]]
    c4 <- (fitZ3)$coeff[[5]]
    c5 <- (fitZ3)$coeff[[6]]
    c6 <- (fitZ3)$coeff[[7]]
    c7 <- (fitZ3)$coeff[[8]]
    c8 <- (fitZ3)$coeff[[9]]
    c9 <- (fitZ3)$coeff[[10]]
    ce0 <- sum3$coefficients[[1,2]]      #Std. Error of c0
    ce1 <- sum3$coefficients[[2,2]]      #Std. Error of c1
    ce2 <- sum3$coefficients[[3,2]]      #Std. Error of c2
    ce3 <- sum3$coefficients[[4,2]]      #Std. Error of c3
    ce4 <- sum3$coefficients[[5,2]]      #Std. Error of c4
    ce5 <- sum3$coefficients[[6,2]]      #Std. Error of c5
    ce6 <- sum3$coefficients[[7,2]]      #Std. Error of c6
    ce7 <- sum3$coefficients[[8,2]]      #Std. Error of c7
    ce8 <- sum3$coefficients[[9,2]]      #Std. Error of c8
    ce9 <- sum3$coefficients[[10,2]]     #Std. Error of c9
    cR2 <- sum3$adj.r.squared            #Adjusted R^2
 
 
 

The Calculation is now done. Here are the ai, bi and ci coefficients :

In [10]:
%%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", "ci")
part0  <- data.frame("0", a0, b0, c0)
part1  <- data.frame("1", a1, b1, c1)
part2  <- data.frame("2", a2, b2, c2)
part3  <- data.frame("3", a3, b3, c3)
part4  <- data.frame("4", a4, b4, c4)
part5  <- data.frame("5", a5, b5, c5)
part6  <- data.frame("6", a6, b6, c6)
part7  <- data.frame("7", a7, b7, c7)
part8  <- data.frame("8", a8, b8, c8)
part9  <- data.frame("9", a9, b9, c9)
partA0 <- data.frame( a0, a1, a2, a3, a4, a5, a6, a7, a8, a9)
partB0 <- data.frame( b0, b1, b2, b3, b4, b5, b6, b7, b8, b9)
partC0 <- data.frame( c0, c1, c2, c3, c4, c5, c6, c7, c8, c9)
partA0
 
 
        a0         a1         a2          a3         a4            a5
1 11.64714 -0.1591489 -0.5569262 0.001088194 0.01016056 -4.141395e-06
             a6       a7            a8            a9
1 -6.603189e-05 0.005218 -1.509284e-05 -4.670043e-05
In [11]:
import numpy as np
import pandas as pd
import pandas.rpy.common as com
 
In [12]:
ai=com.load_data('partA0')
bi=com.load_data('partB0')
ci=com.load_data('partC0')
 
def CreateDataSet(Number=1): 
    Output = []
    Output.extend(zip(ai, bi, ci))
    return Output
    
 
 
 

Coefficients have been calculated, and are under a python format. We first present them an easy way, and will tray to make it in a beautifull table just after.

In [13]:
ai
 
Out[13]:
a0 a1 a2 a3 a4 a5 a6 a7 a8 a9
1 11.647142 -0.159149 -0.556926 0.001088 0.010161 -0.000004 -0.000066 0.005218 -0.000015 -0.000047
In [14]:
bi
 
Out[14]:
b0 b1 b2 b3 b4 b5 b6 b7 b8 b9
1 0.091148 -0.001806 -0.005131 0.000016 0.000101 -5.648576e-08 -0.000001 0.000067 -0 -0.000001
In [15]:
ci
 
Out[15]:
c0 c1 c2 c3 c4 c5 c6 c7 c8 c9
1 0.001315 -0.000024 -0.000079 0 0.000002 -3.437711e-10 -1.121970e-08 0.000001 -2.936975e-09 -9.390577e-09

Now, we verify how much the δpCO2,o we calculate with our formula fits well the "true" one, calculated with the seacarb carbonate chemistry package. The following figure show the error (in %) between estimated and "true" ΔδpCO2 as a function of temperature, for different salinity, in the Mediterranean sea's range (from 13 to 30oC).

In [16]:
%%R
# 8) Verify on Fitted vs. True results
    x1 <- a0 + a1*t + a2*s + a3*t2 + a4*s2 + a5*t3 + a6*s3 + a7*ts + a8*t2s + a9*ts2
    x2 <- b0 + b1*t + b2*s + b3*t2 + b4*s2 + b5*t3 + b6*s3 + b7*ts + b8*t2s + b9*ts2
    x3 <- c0 + c1*t + c2*s + c3*t2 + c4*s2 + c5*t3 + c6*s3 + c7*ts + c8*t2s + c9*ts2
 
 #  ddic = 150
 #  carbd <- carb(flag=15, var1=TA, var2=dic0+ddic*1e-6, S=Salt, T=Temp, P=0, Pt=0, Sit=0, k1k2="l", kf="dr", ks="d", pHscale="T")
 #  dpCO2true <- carbd$pCO2 - pCO2o
 
    dpco2a <- 280
    carbd <- carb(flag=24, var1=pCO2o+dpco2a, var2=TA, S=Salt, T=Temp, P=0, Pt=0, Sit=0, k1k2="l", kf="dr", ks="d", pHscale="T")
    dpCO2true <- carbd$pCO2 - pCO2o
    dDICtrue <- carbd$DIC - dic0
    ddic <- dDICtrue * 1e+6
    ddic2 <- ddic*ddic
    ddic3 <- ddic*ddic2
#   dpCO2fit <- (x0 * ddic) / (1 - x1*ddic)    #This solution only for linear relationship in Sarmiento et al. (1992)
    dpCO2fit <- x1*ddic + x2*ddic2 + x3*ddic3  #This study
#   plot(x=Temp, y=dpCO2true)
#   lines(x=Temp, y=dpCO2fit, col="blue")
 
#   Rough conclusion: relative errors always less than 0.5% when delpCO2 < 280 (i.e., atm pCO2 ~ 560 ppm)
    plot(x=Temp, y=100*(dpCO2fit - dpCO2true)/dpCO2true, col="blue")
 
 
 
In [ ]: