Cell Toolbar:
The
perturbation approach developped in Sarmiento et al.,(1992) identify a
linear relation between the Ocean anthropogenic change in pCO
Below we demonstrate how the coefficients needed to calculate
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. After installing, comment out corresponding line as below.
%%R
# remove.packages('seacarb')
# install.packages('seacarb')
# install.packages('abind')
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 #).
%%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
%%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
Compute corresponding (pre-industrial) DIC field with DIC0 = f(T)
%%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
Increasing DIC incrementally and compute pCO
%%R
DIC <- numeric(0) ; pCO2 <- numeric(0)
delDIC <- numeric(0) ; delpCO2 <- numeric(0)
%%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)
}
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
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
%%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)
%%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
%%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
Show results in a nice table (using Pandas in python).
import numpy as np
import pandas as pd
import pandas.rpy.common as com
ai=com.load_data('partA0')
bi=com.load_data('partB0')
def CreateDataSet(Number=1):
Output = []
Output.extend(zip(ai, bi))
return Output
#print ai.shape
ai
bi
For a doubling of atmospheric xCO2 (280 to 560 ppm), compute relative error in perturbation approach.
%%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)