supplementary_material_CO2_anthr_GLO

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 pCO $_2$ and its ratio with the corresponding ocean change in total carbon :

δ p C O 2 o δ C T = z 0 + z 1 δ p C O 2 o .

$\begin{eqnarray} \frac{\delta p CO_{2_o}}{\delta C_T} = z_0 + z_1 \delta p CO_{2_o} \hbox{.} \end{eqnarray}$ where

z0 $z_0$ and

z1 $z_1$ are quadratic functions of temperature, i.e :

z 0 = a 0 + a 1 T + a 2 T 2,

$\begin{eqnarray} z_0 = a_0 + a_1 T + a_2 T^2 \hbox{,} \end{eqnarray}$ and

z 1 = b 0 + b 1 T + b 2 T 2 .

$\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.

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 pCO $_{2,o}$ , which varies with temperature (due to humidity correction)

In [4]:

%%Rtk     <- Temp + 273.15   #Absolute temperature (K)LNespa <- 20.1050 - 0.0097982*tk - 6163.10/tkespa   <- 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]:

%%Rcarb0 <- 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 pCO $_2$ for each increment until we get $\delta$ pCO $_2$ from 0 to 280 ppm and the corresponding change in surface DIC.

In [6]:

%%RDIC  <- numeric(0) ; pCO2 <- numeric(0)delDIC  <- numeric(0) ; delpCO2 <- numeric(0)
 
 

In [7]:

%%RDIC  <- 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 $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.

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 $a_i$ and $b_i$ 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 $a_i$ and $b_i$ coefficients, and their respective standard error ( $ae_i$ and $be_i$ ):

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 npimport pandas as pdimport 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.shapeai
 
 

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 $\delta pCO_{2,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 $\delta$ pCO $_2$ 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 [ ]:

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

Assess the error in the δpCO2,o\delta pCO_{2,o} calculated from the perturbation approach relative to "true" values from the seacarb carbonate chemistry package.¶

Figure 1: Relative error as a function of temperature for the δ\deltapCO2_2 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.¶

Assess the error in the $\delta pCO_{2,o}$ calculated from the perturbation approach relative to "true" values from the seacarb carbonate chemistry package.¶

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. Relative errors remain less than 0.5% across the observed temperature range of the Mediterranean Sea.¶