    x

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

    x

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 :

    x

δ p C O 2 o δ C T=z 0+z 1 δ p C O 2 o.
where z 0
and z 1
are quadratic functions of temperature, i.e :
z 0=a 0+a 1 T+a 2 T 2,
and
z 1=b 0+b 1 T+b 2 T 2.

Below we demonstrate how the coefficients needed to calculate z 0
and z 1
(i.e : a i
and b i
) have been calculated.

    x

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

In [1]:

    %reload_ext rpy2.ipython

. . .

    x

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

. . .

    x

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.

. . .

    x

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]:

    %%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

. . .

    x

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

. . .

    x

Increasing DIC incrementally and compute pCO 2
for each increment until we get δ
pCO 2
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)

       }

. . .

    x

Compute z 0
and z 1
, i.e., the intercept and slope relating the ratio δ p C O 2 δ C T
and δ 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)

      }

. . .

    x

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

. . .

    x

The calculation is now completed. Below are the a i
and b i
coefficients, and their respective standard error (a e i
and b e 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

. . .

    x

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

. . .

    x

#### Assess the error in the δ p C O 2,o
calculated from the perturbation approach relative to "true" values from the seacarb carbonate chemistry package.¶

    x

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)

. . .

    x

##### Figure 1: Relative error as a function of temperature for the δ
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 [ ]:

. . .
Close Expand Open in Pager Close