    x

## Routine to compute coefficients for anthropogenic CO2 perturbation approach (VAR 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=0+z 1 δ C T+z 2 δ C 2 T+z 3 δ C 3 T.
where z 0
and z 1
are quadratic functions of temperature, i.e :
z 1=a 0+a 1 T+a 2 S+a 3 T 2+a 4 S 2+a 5 T 3+a 6 S 3+a 7 T S+a 8 T 2 S+a 9 T S 2,
z 2=b 0+b 1 T+b 2 S+b 3 T 2+b 4 S 2+b 5 T 3+b 6 S 3+b 7 T S+b 8 T 2 S+b 9 T S 2,
and
z 3=c 0+c 1 T+c 2 S+c 3 T 2+c 4 S 2+c 5 T 3+c 6 S 3+c 7 T S+c 8 T 2 S+c 9 T S 2.

In the following, we will show how the coefficients needed to calculate z 1
, z 2
and z 3
(i.e : a i
, b i
and c i
) have been calculated.

    x

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

. . .

    x

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

. . .

    x

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

. . .

    x

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

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

. . .

    x

We continu, increasing DIC incrementally, and computing pCO 2
for each increment, until we get δ
pCO 2
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"

. . .

    x

Now we compute z 1
, z 2
and z 3
, respectively the intercept, the slope, and the second term coefficient of the cubic relation between Mediterranean δ p C O 2
and δ C T
, 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)

        }

      }

. . .

    x

Then we fit each coefficient to a cubic equation in Temperature and Salinity, so we can deduce the a i
, b i
and c i
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

. . .

    x

The Calculation is now done. Here are the a i
, b i
and c i
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

. . .

    x

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

. . .

    x

Now, we verify how much the δ p C O 2,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" Δ δ
pCO 2
as a function of temperature, for different salinity, in the Mediterranean sea's range (from 13 to 30 o
C).

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

. . .
Close Expand Open in Pager Close