Cell Toolbar:
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 :
In the following, we will show how the coefficients needed to calculate
So let's use R
to compute the partials. Start by loading Rmagic
, which is part of Rpy2
for iPython
.
%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.
%%R
# remove.packages('seacarb')
# install.packages('seacarb')
# remove.packages('Hmisc')
# install.packages('Hmisc')
# install.packages('abind')
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).
%%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)
now we correct for atmospheric humidity (100%) at equilibrium, see Sarmiento et al. (1992 - Eq. 9)
Compute initial field pCO
%%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)
%%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 pCO
%%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)
}
Now we compute
%%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
%%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)
%%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
%%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
import numpy as np
import pandas as pd
import pandas.rpy.common as com
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.
ai
bi
ci
Now, we verify how much the
%%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")