Markus Raitzsch 
Geoscientist
About me
I am a geologist and paleontologist, working in the field of Marine Geosciences. My expertise lies in the chemical analysis of biogenic carbonates, using inductively coupled plasma mass spectrometers (ICPMS). To facilitate the processing of ICPMS raw data, the application of statistical tests and creation of plots, I write R scripts that perform these steps virtually at the same time. Furthermore, I am passionate about wrapping the scripts in Shiny applications for non-R users and to make routine workflows more efficient.
Why R
For my work, R is the most complete language that meets all my demands for data wrangling and analysis, statistical tests, as well as production of high-quality figures. Apart from its availability across all current operating systems, R is open source, making it a great tool to write own scripts and packages, or to customize scripts from published packages. In combination with Shiny and R Markdown, R provides an all-in-one solution for most of my daily work, including data handling, writing reports, and preparing presentations.
Apps & scripts
Pulse
Pulse is a data reduction program for laser ablation time-resolved analysis of carbonates on a Nu Instruments Attom high-resolution ICP-MS. It processes raw files (*.csv), performs an outlier test, and displays the time-resolved data in an interactive plot, where the limits for blank and sample signals can be set via input fields or sliders. The blank-corrected data can then be either treated as a sample, or assigned to one of the carbonate standards JCp-1, JCt-1, or to the glass standards NIST 610, 612, or 614. At the end, element-to-Ca ratios are calculated and drift-corrected using the standard-sample-standard bracketing technique. The app can be tested here, but be aware that the current code is for offline modus, so it is not possible to transfer files between server and local machine.
ScarFace
ScarFace has been developed to facilitate the usage of the R-package seacarb that is used to calculate the carbonate chemistry of seawater requiring a command-line usage. For potential users not familiar with bare code, ScarFace enables to use seacarb via an user interface (ui) without the need for digging into R. For usage on a local computer, it can be downloaded from either Zenodo or from Github. In addition, there is a web version available, which can be used from anywhere online.
MgCaRB
MgCaRB is a tool published in Gray, W. R., & Evans, D. (2019). Nonthermal influences on Mg/Ca in planktonic foraminifera: A review of culture studies and application to the last glacial maximum. Paleoceanography and Paleoclimatology, 34, 306–315. https://doi.org/10.1029/2018PA003517, which calculates temperature from planktic foraminiferal Mg/Ca, accounting for the influences of carbonate chemistry (pH) and salinity. The key component of MgCaRB is that it iteratively solves Mg/Ca and either boron isotopes or atmospheric CO2 (using ice core data) for temperature and pH. Salinity is treated as a known, accounting for Glacial-Interglacial sea level changes. Using Monte-Carlo simulation, MgCaRB propagates uncertainties from Mg/Ca measurement, calibration, salinity, and carbonate chemistry. Age uncertainity is accounted for when sampling the sea level curve or atmospheric CO2. The Shiny web version can be found here.
Bison | sneak preview
Bison is designed to calculate the boron isotope system, perform boron isotope calibrations, as well as calculate pH and pCO2 from the isotopic composition of carbonate. On demand, this app utilizes seacarbx (see below) to account for the effect of Mg and Ca concentrations of seawater on the dissociation constants of carbonic and boric acids, which is important for more accurate deep-time pH and pCO2 reconstructions. Bison is currently under development in cooperation with some colleagues, but a beta version can be tested here. Be aware that some features (e.g., upload own files) might not work properly since it is under development for offline usage first.
seacarbx
This bundle is a seacarb extension, containing a suite of modified seacarb functions to calculate the seawater carbonate system at [Ca2+]sw and [Mg2+]sw different than today. This is particularly useful for deep-time carbonate system calculations, e.g. in Bison, or for culturing studies carried out in artificial or modified seawater. The RDA file also contains MyAMI tabulated parameters defining the temperature and salinity dependencies of the conditional equilibrium constants for [Mg2+] and [Ca2+] in the range 1–60 mM (from Hain et al., 2015, and revised by Hain et al., 2018), as well as a function for bilinear interpolation.
The seacarbx functions in detail:
carbx() is the equivalent to the function carb(), and returns parameters of the seawater carbonate system, based on a given pair of known carbonate system variables.
"carbx" <-
function(flag, var1, var2, S=35, T=25, ca=10.2821e-03, mg=52.8171e-03, Patm=1, P=0, Pt=0, Sit=0, kf='x', ks="d", pHscale="T", b="u74", gas="potential", warn="y", eos="eos80", long=1.e20, lat=1.e20){
n <- max(length(var1), length(var2), length(S), length(T), length(ca), length(mg), length(P), length(Pt), length(Sit), length(kf), length(pHscale), length(ks), length(b))
if(length(flag)!=n){flag <- rep(flag[1],n)}
if(length(var1)!=n){var1 <- rep(var1[1],n)}
if(length(var2)!=n){var2 <- rep(var2[1],n)}
if(length(S)!=n){S <- rep(S[1],n)}
if(length(T)!=n){T <- rep(T[1],n)}
if(length(ca)!=n){ca <- rep(ca[1],n)}
if(length(mg)!=n){mg <- rep(mg[1],n)}
if(length(Patm)!=n){Patm <- rep(Patm[1],n)}
if(length(P)!=n){P <- rep(P[1],n)}
if(length(Pt)!=n){Pt <- rep(Pt[1],n)}
if(length(Sit)!=n){Sit <- rep(Sit[1],n)}
# if(length(k1k2)!=n){k1k2 <- rep(k1k2[1],n)}
if(length(kf)!=n){kf <- rep(kf[1],n)}
if(length(ks)!=n){ks <- rep(ks[1],n)}
if(length(pHscale)!=n){pHscale <- rep(pHscale[1],n)}
if(length(b)!=n){b <- rep(b[1],n)}
# if(length(gas)!=n){gas <- rep(b[1],n)}
# if the concentrations of total silicate and total phosphate are NA
# they are set at 0
Sit[is.na(Sit)] <- 0
Pt[is.na(Pt)] <- 0
# Only two options for eos
if (eos != "teos10" && eos != "eos80")
stop ("invalid parameter eos: ", eos)
# if use of EOS-10 standard
if (eos == "teos10")
{
# Must convert temperature and salinity from TEOS-10 to EOS-80
# convert temperature: from Conservative (CT) to in-situ temperature
# and salinity from Absolute to Practical (SP)
STeos <- teos2eos_geo (S, T, P, long, lat)
InsT <- STeos$T
SP <- STeos$SP
}
else
{
InsT <- T
SP <- S
}
#-------Constants----------------
tk = 273.15; # [K] (for conversion [deg C] <-> [K])
TK = InsT + tk; # TK [K]; T[C]
#---- issues de equic----
Cl = SP / 1.80655; # Cl = chlorinity; S = salinity (per mille)
ST = 0.14 * Cl/96.062 # (mol/kg) total sulfate (Dickson et al., 2007, Table 2)
FLUO = 6.7e-5 * Cl/18.9984 # (mol/kg) total fluoride (Dickson et al., 2007, Table 2)
BOR = bor(S=SP , b=b); # (mol/kg) total boron
#---------------------------------------------------------------------
#--------------------- compute K's ----------------------------------
#---------------------------------------------------------------------
# Ks (free pH scale) at zero pressure and given pressure
Ks_P0 <- Ksx(S=SP, T=InsT, P=0, ca=ca, mg=mg, ks=ks, warn=warn)
Ks <- Ksx(S=SP, T=InsT, P=P, ca=ca, mg=mg, ks=ks, warn=warn)
# Kf on free pH scale
Kff_P0 <- Kf(S=SP, T=InsT, P=0, pHscale="F", kf=kf, Ks_P0, Ks, warn=warn)
Kff <- Kf(S=SP, T=InsT, P=P, pHscale="F", kf=kf, Ks_P0, Ks, warn=warn)
# Kf on given pH scale
Kf <- Kf(S=SP, T=InsT, P=P, pHscale=pHscale, kf=kf, Ks_P0, Ks, warn=warn)
# Conversion factor from total to SWS pH scale at zero pressure
ktotal2SWS_P0 <- kconv(S=SP,T=InsT,P=P,kf=kf,Ks=Ks_P0,Kff=Kff_P0)$ktotal2SWS
# Conversion factor from SWS to chosen pH scale
conv <- kconv(S=SP,T=InsT,P=P,kf=kf,Ks=Ks,Kff=Kff, warn=warn)
kSWS2chosen <- rep(1.,n)
kSWS2chosen [pHscale == "T"] <- conv$kSWS2total [pHscale == "T"]
kSWS2chosen [pHscale == "F"] <- conv$kSWS2free [pHscale == "F"]
K1 <- K1x(S=SP, T=InsT, ca=ca, mg=mg, P=P, pHscale=pHscale, kSWS2chosen, ktotal2SWS_P0, warn=warn)
K2 <- K2x(S=SP, T=InsT, ca=ca, mg=mg, P=P, pHscale=pHscale, kSWS2chosen, ktotal2SWS_P0, warn=warn)
Kw <- Kwx(S=SP, T=InsT, ca=ca, mg=mg, P=P, pHscale=pHscale, kSWS2chosen, warn=warn)
Kb <- Kbx(S=SP, T=InsT, ca=ca, mg=mg, P=P, pHscale=pHscale, kSWS2chosen, ktotal2SWS_P0, warn=warn)
K1p <- K1p(S=SP, T=InsT, P=P, pHscale=pHscale, kSWS2chosen, warn=warn)
K2p <- K2p(S=SP, T=InsT, P=P, pHscale=pHscale, kSWS2chosen, warn=warn)
K3p <- K3p(S=SP, T=InsT, P=P, pHscale=pHscale, kSWS2chosen, warn=warn)
Ksi <- Ksi(S=SP, T=InsT, P=P, pHscale=pHscale, kSWS2chosen, warn=warn)
Kspa <- Kspax(S=SP, T=InsT, ca=ca, mg=mg, P=P, warn=warn)
Kspc <- Kspcx(S=SP, T=InsT, ca=ca, mg=mg, P=P, warn=warn)
rho <- rho(S=SP,T=InsT,P=P)
# Compute "standard" K0 with S, in situ T, and atmospheric pressure
K0 <- K0x(S=SP, T=InsT, Patm=Patm, P=0, ca=ca, mg=mg, warn=warn)
# Compute potential K0 with S, potential temperature, and atmospheric pressure (usually 1 atm)
K0pot <- K0x(S=SP, T=theta(S=SP, T=InsT, P=P, Pref=0), Patm=Patm, P=0, ca=ca, mg=mg, warn=warn)
# Compute in situ K0 with S, in situ temperature, and total pressure pressure (atmospheric + hydrostatic)
K0insitu <- K0x(S=SP, T=InsT, Patm=Patm, P=P, ca=ca, mg=mg, warn=warn)
#------------------------------------------------------------------#
#------------------------------------------------------------------#
# VARIABLES #
#------------------------------------------------------------------#
#------------------------------------------------------------------#
# flag = 1 pH-CO2 given
# flag = 2 CO2-HCO3 given
# flag = 3 CO2-CO3 given
# flag = 4 CO2-ALK given
# flag = 5 CO2-DIC given
# flag = 6 pH and HCO3 given
# flag = 7 pH and CO3 given
# flag = 8 pH and ALK given
# flag = 9 pH and DIC given
# flag = 10 HCO3 and CO3 given
# flag = 11 HCO3 and ALK given
# flag = 12 HCO3 and DIC given
# flag = 13 CO3 and ALK given
# flag = 14 CO3 and DIC given
# flag = 15 ALK and DIC given
# flag = 21 pCO2-pH given
# flag = 22 pCO2-HCO3 given
# flag = 23 pCO2-CO3 given
# flag = 24 pCO2-ALK given
# flag = 25 pCO2-DIC given
# Initialise output vectors
H <- rep(NA, n)
PH <- rep(NA, n)
CO2 <- rep(NA, n)
pCO2 <- rep(NA, n)
fCO2 <- rep(NA, n)
HCO3 <- rep(NA, n)
CO3 <- rep(NA, n)
DIC <- rep(NA, n)
ALK <- rep(NA, n)
# ------------ case 1.) PH and CO2 given ----
# Indices of flag elements where flag = 1
i_flag_1 <- which (flag == 1)
PH[i_flag_1] <- var1[i_flag_1]
CO2[i_flag_1] <- var2[i_flag_1]
h <- 10^(-PH[i_flag_1])
fCO2[i_flag_1] <- CO2[i_flag_1] / K0[i_flag_1]
HCO3[i_flag_1] <- (K1[i_flag_1] * CO2[i_flag_1]) / h
CO3[i_flag_1] <- (K2[i_flag_1] * HCO3[i_flag_1]) / h
DIC[i_flag_1] <- CO2[i_flag_1] + HCO3[i_flag_1] + CO3[i_flag_1]
H[i_flag_1] <- h
# ------------ case 2.) CO2 and HCO3 given ----
# Indices of flag elements where flag = 2
i_flag_2 <- which (flag == 2)
CO2[i_flag_2] <- var1[i_flag_2]
HCO3[i_flag_2] <- var2[i_flag_2]
fCO2[i_flag_2] <- CO2[i_flag_2] / K0[i_flag_2]
h <- K0[i_flag_2] * K1[i_flag_2] * fCO2[i_flag_2] / HCO3[i_flag_2]
CO3[i_flag_2] <- K0[i_flag_2] * K1[i_flag_2] * K2[i_flag_2] * fCO2[i_flag_2] / (h^2)
DIC[i_flag_2] <- CO2[i_flag_2] + HCO3[i_flag_2] + CO3[i_flag_2]
PH[i_flag_2] <- -log10(h)
H[i_flag_2] <- h
# ------------ case 3.) CO2 and CO3 given ----
# Indices of flag elements where flag = 3
i_flag_3 <- which (flag == 3)
CO2[i_flag_3] <- var1[i_flag_3]
CO3[i_flag_3] <- var2[i_flag_3]
fCO2[i_flag_3] <- CO2[i_flag_3] / K0[i_flag_3]
k1co2 <- K1[i_flag_3] * CO2[i_flag_3]
h <- sqrt((K2[i_flag_3]*k1co2) / CO3[i_flag_3])
HCO3[i_flag_3] <- k1co2 / h
DIC[i_flag_3] <- CO2[i_flag_3] + HCO3[i_flag_3] + CO3[i_flag_3]
PH[i_flag_3] <- -log10(h)
H[i_flag_3] <- h
# ------------ case 4.) CO2 and ALK given ----
# Indices of flag elements where flag = 4
i_flag_4 <- which (flag == 4)
CO2[i_flag_4] <- var1[i_flag_4]
ALK[i_flag_4] <- var2[i_flag_4]
fALK <- function(x)# K1=K1, K2=K2, CO2=CO2, BOR=BOR, Kb=Kb, Kw=Kw, Pt=Pt, K1p=K1p, K2p=K2p, K3p=K3p, Sit=Sit, Ksi=Ksi, ST=ST, Ks=Ks, FLUO=FLUO, Kf=Kf, ALK=ALK)
{
# components for ALK
dic <- co2*(1+K1_i/x+K1_i*K2_i/(x*x))
hco3 <- dic*x*K1_i/(x*x + K1_i*x + K1_i*K2_i)
co3 <- dic*K1_i*K2_i/(x*x + K1_i*x + K1_i*K2_i)
boh4 <- bor/(1+x/Kb_i)
oh <- Kw_i/x
h3po4 <- Pt_i*x^3/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
hpo4 <- Pt_i*K1p_i*K2p_i*x/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
po4 <- Pt_i*K1p_i*K2p_i*K3p_i/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
siooh3 <- Sit_i/(1+x/Ksi_i)
## calculate Hfree and Htot
if(phs=="F")
{
hfree <- x ## if pHscale = free scale
htot <- x / (1+st/Ks_i)
}
else if(phs=="T")
{
hfree <- x * (1+st/Ks_i)
htot <- x
}
else if(phs=="SWS")
{
hfree <- x * (1 + st/Ks_i + fluo/Kff_i)
htot <- hfree / (1+st/Ks_i)
}
hso4 <- st/(1+Ks_i/hfree)
hf <- fluo/(1+Kf_i/htot)
############
OUT <- hco3+2*co3+boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4-alk
OUT
}
# Calculate [H+] from [CO2] and total alk
h <- rep(NA, length(i_flag_4))
j <- 1
for(i in (i_flag_4))
{
# Parameters used by function fAlk defined above and called below through uniroot()
K1_i <- K1[i]
K2_i <- K2[i]
K1p_i <- K1p[i]
K2p_i <- K2p[i]
K3p_i <- K3p[i]
Kb_i <- Kb[i]
Kw_i <- Kw[i]
Ksi_i <- Ksi[i]
Ks_i <- Ks[i]
Kf_i <- Kf[i]
Kff_i <- Kff[i]
Sit_i <- Sit[i]
Pt_i <- Pt[i]
co2 <- CO2[i]
alk <- ALK[i]
st <- ST[i]
bor <- BOR[i]
fluo <- FLUO[i]
phs <- pHscale[i]
# Calculate [H+] from total alk
h[j] <- uniroot(fALK,c(10^(-9.5),10^(-3.5)), tol=1e-20)$root
j <- j + 1
}
DIC[i_flag_4] <- CO2[i_flag_4] * (1+K1[i_flag_4]/h + K1[i_flag_4]*K2[i_flag_4]/(h*h))
temp <- h*h + K1[i_flag_4]*h + K1[i_flag_4]*K2[i_flag_4]
HCO3[i_flag_4] <- (DIC[i_flag_4]*K1[i_flag_4]*h) / temp
CO3[i_flag_4] <- (DIC[i_flag_4]*K1[i_flag_4]*K2[i_flag_4]) / temp
fCO2[i_flag_4] <- CO2[i_flag_4] / K0[i_flag_4]
PH[i_flag_4] <- -log10(h)
H[i_flag_4] <- h
# ------------ case 5.) CO2 and DIC given ----
# Indices of flag elements where flag = 5
i_flag_5 <- which (flag == 5)
CO2[i_flag_5] <- var1[i_flag_5]
DIC[i_flag_5] <- var2[i_flag_5]
fCO2[i_flag_5] <- CO2[i_flag_5] / K0[i_flag_5]
a <- K1[i_flag_5] * K2[i_flag_5] * CO2[i_flag_5]
b <- K1[i_flag_5] * CO2[i_flag_5]
c <- CO2[i_flag_5] - DIC[i_flag_5]
D <- b*b - 4*a*c
h <- (2*a) / (sqrt(D)-b)
HCO3[i_flag_5] <- K0[i_flag_5] * K1[i_flag_5] * fCO2[i_flag_5] / h
CO3[i_flag_5] <- DIC[i_flag_5] - CO2[i_flag_5] - HCO3[i_flag_5]
PH[i_flag_5] <- -log10(h)
H[i_flag_5] <- h
# ------------ case 6.) PH and HCO3 given ----
# Indices of flag elements where flag = 6
i_flag_6 <- which (flag == 6)
PH[i_flag_6] <- var1[i_flag_6]
HCO3[i_flag_6] <- var2[i_flag_6]
h <- 10^(-PH[i_flag_6])
CO2[i_flag_6] <- (HCO3[i_flag_6] * h)/K1[i_flag_6]
CO3[i_flag_6] <- K2[i_flag_6] * HCO3[i_flag_6] /h
DIC[i_flag_6] <- CO2[i_flag_6] + HCO3[i_flag_6] + CO3[i_flag_6]
fCO2[i_flag_6] <- CO2[i_flag_6]/K0[i_flag_6]
H[i_flag_6] <- h
# ------------ case 7.) PH and CO3 given ----
# Indices of flag elements where flag = 7
i_flag_7 <- which (flag == 7)
PH[i_flag_7] <- var1[i_flag_7]
CO3[i_flag_7] <- var2[i_flag_7]
h <- 10^(-PH[i_flag_7])
HCO3[i_flag_7] <- CO3[i_flag_7] * h/K2[i_flag_7]
CO2[i_flag_7] <- HCO3[i_flag_7] * h/K1[i_flag_7]
fCO2[i_flag_7] <- CO2[i_flag_7]/K0[i_flag_7]
DIC[i_flag_7] <- CO2[i_flag_7] + HCO3[i_flag_7] + CO3[i_flag_7]
H[i_flag_7] <- h
# ------------ case 8.) PH and ALK given ----
# Indices of flag elements where flag = 8
i_flag_8 <- which(flag == 8)
PH[i_flag_8] <- var1[i_flag_8]
ALK[i_flag_8] <- var2[i_flag_8]
h <- 10^(-PH[i_flag_8])
H[i_flag_8] <- h
## calculate Hfree anf Htot
hfree <- rep(NA, length(i_flag_8))
htot <- rep(NA, length(i_flag_8))
sc <- pHscale[i_flag_8]
st <- ST[i_flag_8]
ks <- Ks[i_flag_8]
fluo <- FLUO[i_flag_8]
kff <- Kff[i_flag_8]
# Where pHscale=="F", pHscale = free scale
i_sc_F <- which (sc == "F")
hfree[i_sc_F] <- h[i_sc_F]
htot[i_sc_F] <- h[i_sc_F] / (1+st[i_sc_F]/ks[i_sc_F])
# Where pHscale=="T", pHscale = total scale
i_sc_T <- which (sc == "T")
hfree[i_sc_T] <- h[i_sc_T] * (1+st[i_sc_T]/ks[i_sc_T])
htot[i_sc_T] <- h[i_sc_T]
# Where pHscale=="SWS", pHscale = SW scale
i_sc_S <- which (sc == "SWS")
hfree[i_sc_S] <- h[i_sc_S] * (1 + st[i_sc_S]/ks[i_sc_S] + fluo[i_sc_S]/kff[i_sc_S])
htot[i_sc_S] <- hfree[i_sc_S] / (1+st[i_sc_S]/ks[i_sc_S])
# Calculate some invariable components of total alkalinity
boh4 <- BOR[i_flag_8] / (1+h/Kb[i_flag_8])
oh <- Kw[i_flag_8]/h
temp <- (h^3 + K1p[i_flag_8]*h^2 + K1p[i_flag_8]*K2p[i_flag_8]*h +
K1p[i_flag_8]*K2p[i_flag_8]*K3p[i_flag_8])
h3po4 <- Pt[i_flag_8]*h^3 / temp
hpo4 <- Pt[i_flag_8]*K1p[i_flag_8]*K2p[i_flag_8]*h / temp
po4 <- Pt[i_flag_8]*K1p[i_flag_8]*K2p[i_flag_8]*K3p[i_flag_8] / temp
siooh3 <- Sit[i_flag_8] / (1+h/Ksi[i_flag_8])
hso4 <- st/(1+ks/hfree)
hf <- fluo / (1+Kf[i_flag_8]/htot)
# Sum of these components (partial alkalinity)
alk_p <- boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4
# As [HCO3-] and [CO3--] are both linearly dependent on DIC
# hco3 = co2*K1/h = DIC/(1+K1/h+K1*K2/(h^2)) * K1/h
# 2*co3 = 2*hco3*K2/h = DIC/(1+K1/h+K1*K2/(h^2)) * 2*K1*K2/h^2
# carbonate alk ([HCO3-] + 2*[CO3--]) is also linearly dependent on DIC
# Calculate DIC from carbonate alk (total alk - partial alk)
temp <- h*h + K1[i_flag_8]*h + K1[i_flag_8]*K2[i_flag_8]
DIC[i_flag_8] <- (ALK[i_flag_8]-alk_p) * temp /
(h*K1[i_flag_8] + 2*K1[i_flag_8]*K2[i_flag_8])
CO2[i_flag_8] <- DIC[i_flag_8] /
(1 + K1[i_flag_8]/h + K1[i_flag_8]*K2[i_flag_8]/(h^2))
HCO3[i_flag_8] <- CO2[i_flag_8]*K1[i_flag_8]/h
CO3[i_flag_8] <- HCO3[i_flag_8]*K2[i_flag_8]/h
fCO2[i_flag_8] <- CO2[i_flag_8]/K0[i_flag_8]
# ------------ case 9.) PH and DIC given ----
# Indices of flag elements where flag = 9
i_flag_9 <- which (flag == 9)
PH[i_flag_9] <- var1[i_flag_9]
DIC[i_flag_9] <- var2[i_flag_9]
h <- 10^(-PH[i_flag_9])
H[i_flag_9] <- h
temp <- h*h + K1[i_flag_9]*h + K1[i_flag_9]*K2[i_flag_9]
HCO3[i_flag_9] <- (DIC[i_flag_9]*K1[i_flag_9]*h) / temp
CO3[i_flag_9] <- (DIC[i_flag_9]*K1[i_flag_9]*K2[i_flag_9]) / temp
CO2[i_flag_9] <- h*HCO3[i_flag_9]/K1[i_flag_9]
fCO2[i_flag_9] <- CO2[i_flag_9]/K0[i_flag_9]
# ------------ case 10.) HCO3 and CO3 given ----
# Indices of flag elements where flag = 10
i_flag_10 <- which (flag == 10)
HCO3[i_flag_10] <- var1[i_flag_10]
CO3[i_flag_10] <- var2[i_flag_10]
h <- K2[i_flag_10]*HCO3[i_flag_10]/CO3[i_flag_10]
CO2[i_flag_10] <- h*HCO3[i_flag_10]/K1[i_flag_10]
DIC[i_flag_10] <- CO2[i_flag_10] + HCO3[i_flag_10] + CO3[i_flag_10]
fCO2[i_flag_10] <- CO2[i_flag_10]/K0[i_flag_10]
PH[i_flag_10] <- -log10(h)
H[i_flag_10] <- h
# ------------ case 11.) HCO3 and ALK given ----
# Indices of flag elements where flag = 11
i_flag_11 <- which (flag == 11)
HCO3[i_flag_11] <- var1[i_flag_11]
ALK[i_flag_11] <- var2[i_flag_11]
fALK <- function(x)# K1=K1, K2=K2, HCO3=HCO3, BOR=BOR, Kb=Kb, Kw=Kw, Pt=Pt, K1p=K1p, K2p=K2p, K3p=K3p, Sit=Sit, Ksi=Ksi, ST=ST, Ks=Ks, FLUO=FLUO, Kf=Kf, ALK=ALK) {
{
# components for ALK
dic <- hco3*(x^2 + K1_i*x + K1_i*K2_i) / (K1_i*x)
co3 <- dic*K1_i*K2_i / (x*x + K1_i*x + K1_i*K2_i)
boh4 <- bor / (1+x/Kb_i)
oh <- Kw_i/x
h3po4 <- Pt_i*x^3 / (x^3 + K1p_i*x^2 + K1p_i*K2p_i*x + K1p_i*K2p_i*K3p_i)
hpo4 <- Pt_i*K1p_i*K2p_i*x/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
po4 <- Pt_i*K1p_i*K2p_i*K3p_i/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
siooh3 <- Sit_i/(1+x/Ksi_i)
## calculate Hfree and Htot
if(phs=="F")
{
hfree <- x ## if pHscale = free scale
htot <- x / (1+st/Ks_i)
}
else if(phs=="T")
{
hfree <- x * (1+st/Ks_i)
htot <- x
}
else if(phs=="SWS")
{
hfree <- x * (1 + st/Ks_i + fluo/Kff_i)
htot <- hfree / (1+st/Ks_i)
}
hso4 <- st/(1+Ks_i/hfree)
hf <- fluo/(1+Kf_i/htot)
############
OUT <- hco3+2*co3+boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4-alk
OUT
}
# Calculate [H+] from [HCO3] and total alk
h <- rep(NA, length(i_flag_11))
j <- 1
for(i in (i_flag_11))
{
# Parameters used by function fAlk defined above and called below through uniroot()
K1_i <- K1[i]
K2_i <- K2[i]
K1p_i <- K1p[i]
K2p_i <- K2p[i]
K3p_i <- K3p[i]
Kb_i <- Kb[i]
Kw_i <- Kw[i]
Ksi_i <- Ksi[i]
Ks_i <- Ks[i]
Kf_i <- Kf[i]
Kff_i <- Kff[i]
Sit_i <- Sit[i]
Pt_i <- Pt[i]
hco3 <- HCO3[i]
alk <- ALK[i]
st <- ST[i]
bor <- BOR[i]
fluo <- FLUO[i]
phs <- pHscale[i]
# Calculate [H+] from total alk
h[j] <- uniroot(fALK,c(10^(-9.5),10^(-3)), tol=1e-20)$root
j <- j + 1
}
CO2[i_flag_11] <- h*HCO3[i_flag_11]/K1[i_flag_11]
CO3[i_flag_11] <- K2[i_flag_11]*HCO3[i_flag_11]/h
DIC[i_flag_11] <- CO2[i_flag_11] + HCO3[i_flag_11] + CO3[i_flag_11]
PH[i_flag_11] <- -log10(h)
H[i_flag_11] <- h
fCO2[i_flag_11] <- CO2[i_flag_11]/K0[i_flag_11]
# ------------ case 12.) HCO3 and DIC given
# Indices of flag elements where flag = 12
i_flag_12 <- which (flag == 12)
HCO3[i_flag_12] <- var1[i_flag_12]
DIC[i_flag_12] <- var2[i_flag_12]
a <- HCO3[i_flag_12]
b <- K1[i_flag_12]*(HCO3[i_flag_12]-DIC[i_flag_12])
c <- K1[i_flag_12]*K2[i_flag_12]*HCO3[i_flag_12]
D <- b*b - 4*a*c
h <- (-b-sqrt(D))/(2*a)
CO2[i_flag_12] <- h*HCO3[i_flag_12]/K1[i_flag_12]
CO3[i_flag_12] <- K2[i_flag_12]*HCO3[i_flag_12]/h
fCO2[i_flag_12] <- CO2[i_flag_12]/K0[i_flag_12]
PH[i_flag_12] <- -log10(h)
H[i_flag_12] <- h
# ------------ case 13.) CO3 and ALK given ----
# Indices of flag elements where flag = 13
i_flag_13 <- which (flag == 13)
CO3[i_flag_13] <- var1[i_flag_13]
ALK[i_flag_13] <- var2[i_flag_13]
fALK <- function(x)# K1=K1, K2=K2, HCO3=HCO3, BOR=BOR, Kb=Kb, Kw=Kw, Pt=Pt, K1p=K1p, K2p=K2p, K3p=K3p, Sit=Sit, Ksi=Ksi, ST=ST, Ks=Ks, FLUO=FLUO, Kf=Kf, ALK=ALK) {
{
# composants for ALK
dic <- co3*(x^2 + K1_i*x + K1_i*K2_i)/(K1_i*K2_i)
hco3 <- dic*K1_i*x/(x*x + K1_i*x + K1_i*K2_i)
boh4 <- bor/(1+x/Kb_i)
oh <- Kw_i/x
h3po4 <- Pt_i*x^3 / (x^3 + K1p_i*x^2 + K1p_i*K2p_i*x + K1p_i*K2p_i*K3p_i)
hpo4 <- Pt_i*K1p_i*K2p_i*x/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
po4 <- Pt_i*K1p_i*K2p_i*K3p_i/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
siooh3 <- Sit_i/(1+x/Ksi_i)
## calculate Hfree and Htot
if(phs=="F")
{
hfree <- x ## if pHscale = free scale
htot <- x / (1+st/Ks_i)
}
else if(phs=="T")
{
hfree <- x * (1+st/Ks_i)
htot <- x
}
else if(phs=="SWS")
{
hfree <- x * (1 + st/Ks_i + fluo/Kff_i)
htot <- hfree / (1+st/Ks_i)
}
hso4 <- st/(1+Ks_i/hfree)
hf <- fluo/(1+Kf_i/htot)
############
OUT <- hco3+2*co3+boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4-alk
OUT
}
# Calculate [H+] from [CO3] and total alk
h <- rep(NA, length(i_flag_13))
j <- 1
for(i in (i_flag_13))
{
# Parameters used by function fAlk defined above and called below through uniroot()
K1_i <- K1[i]
K2_i <- K2[i]
K1p_i <- K1p[i]
K2p_i <- K2p[i]
K3p_i <- K3p[i]
Kb_i <- Kb[i]
Kw_i <- Kw[i]
Ksi_i <- Ksi[i]
Ks_i <- Ks[i]
Kf_i <- Kf[i]
Kff_i <- Kff[i]
Sit_i <- Sit[i]
Pt_i <- Pt[i]
co3 <- CO3[i]
alk <- ALK[i]
st <- ST[i]
bor <- BOR[i]
fluo <- FLUO[i]
phs <- pHscale[i]
# Calculate [H+] from total alk
#h[j] <- uniroot(fALK,c(10^(-9.5),10^(-3.5)), tol=1e-20)$root
h[j] <- tryCatch(uniroot(fALK, c(1e-10, 10^(-3.5)), tol = 1e-30)$root, error = function(e) NA)
j <- j + 1
}
HCO3[i_flag_13] <- h*CO3[i_flag_13]/K2[i_flag_13]
CO2[i_flag_13] <- h*HCO3[i_flag_13]/K1[i_flag_13]
fCO2[i_flag_13] <- CO2[i_flag_13]/K0[i_flag_13]
DIC[i_flag_13] <- HCO3[i_flag_13] + CO2[i_flag_13] + CO3[i_flag_13]
PH[i_flag_13] <- -log10(h)
H[i_flag_13] <- h
# ------------ case 14.) CO3 and DIC given ----
# Indices of flag elements where flag = 14
i_flag_14 <- which (flag == 14)
CO3[i_flag_14] <- var1[i_flag_14]
DIC[i_flag_14] <- var2[i_flag_14]
a <- CO3[i_flag_14]
b <- K1[i_flag_14]*CO3[i_flag_14]
c <- K1[i_flag_14] * K2[i_flag_14] * (CO3[i_flag_14]-DIC[i_flag_14])
D <- b*b - 4*a*c
h <- (-b+sqrt(D))/(2*a)
HCO3[i_flag_14] <- h*CO3[i_flag_14]/K2[i_flag_14]
CO2[i_flag_14] <- h*HCO3[i_flag_14]/K1[i_flag_14]
fCO2[i_flag_14] <- CO2[i_flag_14]/K0[i_flag_14]
PH[i_flag_14] <- -log10(h)
H[i_flag_14] <- h
# ------------ case 15.) ALK and DIC given ----
# Indices of flag elements where flag = 15
i_flag_15 <- which(flag == 15)
ALK[i_flag_15] <- var1[i_flag_15]
DIC[i_flag_15] <- var2[i_flag_15]
fALK <- function(x) # K1=K1, K2=K2, DIC=DIC, BOR=BOR, Kb=Kb, Kw=Kw, Pt=Pt, K1p=K1p, K2p=K2p, K3p=K3p, Sit=Sit, Ksi=Ksi, ST=ST, Ks=Ks, FLUO=FLUO, Kf=Kf, ALK=ALK) {
{
# composants for ALK
hco3 <- dic*x*K1_i/(x*x + K1_i*x + K1_i*K2_i)
co3 <- dic*K1_i*K2_i/(x*x + K1_i*x + K1_i*K2_i)
boh4 <- bor/(1+x/Kb_i)
oh <- Kw_i/x
h3po4 <- Pt_i*x^3 / (x^3 + K1p_i*x^2 + K1p_i*K2p_i*x + K1p_i*K2p_i*K3p_i)
hpo4 <- Pt_i*K1p_i*K2p_i*x/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
po4 <- Pt_i*K1p_i*K2p_i*K3p_i/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
siooh3 <- Sit_i/(1+x/Ksi_i)
## calculate Hfree and Htot
if(phs=="F")
{
hfree <- x ## if pHscale = free scale
htot <- x / (1+st/Ks_i)
}
else if(phs=="T")
{
hfree <- x * (1+st/Ks_i)
htot <- x
}
else if(phs=="SWS")
{
hfree <- x * (1 + st/Ks_i + fluo/Kff_i)
htot <- hfree / (1+st/Ks_i)
}
hso4 <- st/(1+Ks_i/hfree)
hf <- fluo/(1+Kf_i/htot)
############
OUT <- hco3+2*co3+boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4-alk
OUT
}
# Calculate [H+] from DIC and total alk
h <- rep(NA, length(i_flag_15))
j <- 1
for(i in (i_flag_15))
{
# Parameters used by function fAlk defined above and called below through uniroot()
K1_i <- K1[i]
K2_i <- K2[i]
K1p_i <- K1p[i]
K2p_i <- K2p[i]
K3p_i <- K3p[i]
Kb_i <- Kb[i]
Kw_i <- Kw[i]
Ksi_i <- Ksi[i]
Ks_i <- Ks[i]
Kf_i <- Kf[i]
Kff_i <- Kff[i]
Sit_i <- Sit[i]
Pt_i <- Pt[i]
dic <- DIC[i]
alk <- ALK[i]
st <- ST[i]
bor <- BOR[i]
fluo <- FLUO[i]
phs <- pHscale[i]
# Calculate [H+] from total alk
#h[j] <- uniroot(fALK,c(1e-10,10^(-3.5)), tol=1e-30)$root
h[j] <- tryCatch(uniroot(fALK, c(1e-10, 10^(-3.5)), tol = 1e-30)$root, error = function(e) NA)
j <- j + 1
}
temp <- h*h + K1[i_flag_15]*h + K1[i_flag_15]*K2[i_flag_15]
HCO3[i_flag_15] <- (DIC[i_flag_15]*K1[i_flag_15]*h) / temp
CO3[i_flag_15] <- (DIC[i_flag_15]*K1[i_flag_15]*K2[i_flag_15]) / temp
CO2[i_flag_15] <- h*HCO3[i_flag_15]/K1[i_flag_15]
fCO2[i_flag_15] <- CO2[i_flag_15]/K0[i_flag_15]
PH[i_flag_15] <- -log10(h)
H[i_flag_15] <- h
#Initialize new 'potential' & "insitu' arrays (this are overwritten below with correct values)
fCO2pot <- fCO2
pCO2pot <- pCO2
fCO2insitu <- fCO2
pCO2insitu <- pCO2
# ------------ Compute pCO2 for cases 1 to 15 (pCO2, pCO2insitu, & pCO2pot are NOT input variables)
# Indices of flag elements where 1 <= flag <= 15
i_flag <- which (flag >= 1 & flag <= 15)
# 1) Classic pCO2: compute from in situ T and surface P (Patm)
tk <- TK[i_flag]
B <- -1636.75 + 12.0408*tk - 0.0327957*(tk*tk) + 0.0000316528*(tk*tk*tk);
fugfac <- exp((Patm[i_flag])*(B + 2*((1-fCO2[i_flag])^2)*(57.7-0.118*tk))/(82.05736*tk))
pCO2[i_flag] <- fCO2[i_flag] / fugfac
# 2) In-situ pCO2: compute from in situ T and total P (Patm + Phydrostatic)
fCO2insitu[i_flag] <- fCO2[i_flag] * K0[i_flag] / K0insitu[i_flag]
xCO2approx <- pCO2[i_flag] #Results virtually insensitive to this (do not use fCO2insitu)
fugfacinsitu <- exp((Patm[i_flag] + P[i_flag]/1.01325)*(B + 2*((1-xCO2approx)^2)*(57.7-0.118*tk))/(82.05736*tk))
pCO2insitu[i_flag] <- fCO2insitu[i_flag] / fugfacinsitu
# 3) Potential pCO2: compute from potential T & surface P (Patm only)
tkp <- theta(S=SP[i_flag], T=InsT[i_flag], P=P[i_flag], Pref=0) + 273.15 #Potential temperature in Kelvin
Bpot <- -1636.75 + 12.0408*tkp - 0.0327957*(tkp*tkp) + 0.0000316528*(tkp*tkp*tkp);
fCO2pot[i_flag] <- fCO2[i_flag] * K0[i_flag] / K0pot[i_flag]
fugfacpot <- exp((Patm[i_flag])*(Bpot + 2*((1-fCO2pot[i_flag])^2)*(57.7-0.118*tkp))/(82.05736*tkp))
pCO2pot[i_flag] <- fCO2pot[i_flag] / fugfacpot
# ------------ Compute fCO2 for cases 21 to 25 (pCO2, pCO2pot or pCO2insitu is an input variable)
# Indices of flag elements where 21 <= flag <= 25
i_flag <- which (flag >= 21 & flag <= 25)
tk <- TK[i_flag]
tkp <- theta(S=SP[i_flag], T=InsT[i_flag], P=P[i_flag], Pref=0) + 273.15 #Potential temperature in Kelvin
B <- -1636.75 + 12.0408*tk - 0.0327957*(tk*tk) + 0.0000316528*(tk*tk*tk);
Bpot <- -1636.75 + 12.0408*tkp - 0.0327957*(tkp*tkp) + 0.0000316528*(tkp*tkp*tkp);
if(gas=="insitu")
{
# From in situ pCO2 (in situ T, atm + hydro P), compute potential and standard fCO2 and pCO2
# a) define input as pCO2insitu & compute fCO2insitu
pCO2insitu[i_flag] <- var1[i_flag] * 1e-6
# xCO2approx <- pCO2insitu[i_flag] (this would be a very gross overestimate at say 5000 m)
xCO2approx <- 0. #This poor approximation has virtually no effect on computed result
fugfacinsitu <- exp((Patm[i_flag]+P[i_flag]/1.01325)*(B + 2*((1-xCO2approx)^2)*(57.7-0.118*tk)) /(82.05736*tk))
fCO2insitu[i_flag] <- pCO2insitu[i_flag] * fugfacinsitu
# b) compute fCO2pot & fCO2
fCO2pot[i_flag] <- fCO2insitu[i_flag] * K0insitu[i_flag] / K0pot[i_flag]
fCO2[i_flag] <- fCO2insitu[i_flag] * K0insitu[i_flag] / K0[i_flag]
# c) compute pCO2pot & pCO2
fugfac <- exp((Patm[i_flag]) *(B + 2*((1-fCO2[i_flag])^2) *(57.7-0.118*tk)) /(82.05736*tk))
fugfacpot <- exp((Patm[i_flag]) *(Bpot + 2*((1-fCO2pot[i_flag])^2) *(57.7-0.118*tkp))/(82.05736*tkp))
pCO2pot[i_flag] <- fCO2pot[i_flag] / fugfacpot
pCO2[i_flag] <- fCO2[i_flag] / fugfac
}
else if(gas=="potential")
{
# From potential pCO2 (potential T, atm P), compute standard and in situ fCO2 and pCO2
# a) define input as pCO2pot & compute fCO2pot
pCO2pot[i_flag] <- var1[i_flag] * 1e-6
fugfacpot <- exp((Patm[i_flag] )*(Bpot + 2*((1-pCO2pot[i_flag])^2)*(57.7-0.118*tkp))/(82.05736*tkp))
fCO2pot[i_flag] <- pCO2pot[i_flag] * fugfacpot
# b) compute fCO2 & fCO2insitu
fCO2[i_flag] <- fCO2pot[i_flag] * K0pot[i_flag] / K0[i_flag]
fCO2insitu[i_flag] <- fCO2pot[i_flag] * K0pot[i_flag] / K0insitu[i_flag]
# c) Compute pCO2 & pCO2insitu
fugfac <- exp((Patm[i_flag] )*(B + 2*((1-fCO2[i_flag])^2)*(57.7-0.118*tk))/(82.05736*tk))
fugfacinsitu <- exp((Patm[i_flag]+P[i_flag]/1.01325)*(B + 2*((1-fCO2[i_flag])^2)*(57.7-0.118*tk))/(82.05736*tk))
pCO2[i_flag] <- fCO2[i_flag] / fugfac
pCO2insitu[i_flag] <- fCO2insitu[i_flag] / fugfacinsitu
}
else if(gas=="standard")
{
# From standard pCO2 (in situ T, atm P), compute potential and in situ fCO2 and pCO2
# a) define input as pCO2 & compute fCO2
pCO2[i_flag] <- var1[i_flag] * 1e-6
fugfac <- exp((Patm[i_flag] )*(B + 2*((1-pCO2[i_flag])^2)*(57.7-0.118*tk))/(82.05736*tk))
fCO2[i_flag] <- pCO2[i_flag] * fugfac
# b) compute fCO2pot & fCO2insitu
fCO2pot[i_flag] <- fCO2[i_flag] * K0[i_flag] / K0pot[i_flag]
fCO2insitu[i_flag] <- fCO2[i_flag] * K0[i_flag] / K0insitu[i_flag]
# c) Compute pCO2pot & pCO2insitu
fugfacpot <- exp((Patm[i_flag] )*(Bpot + 2*((1-pCO2[i_flag])^2)*(57.7-0.118*tkp))/(82.05736*tkp))
fugfacinsitu <- exp((Patm[i_flag]+P[i_flag]/1.01325)*(B + 2*((1-pCO2[i_flag])^2)*(57.7-0.118*tk))/(82.05736*tk))
pCO2pot[i_flag] <- fCO2pot[i_flag] / fugfacpot
pCO2insitu[i_flag] <- fCO2insitu[i_flag] / fugfacinsitu
}
# ------------ case 21.) PH and pCO2 given ----
# Indices of flag elements where flag = 21
i_flag_21 <- which (flag == 21)
PH[i_flag_21] <- var2[i_flag_21]
h <- 10^(-PH[i_flag_21])
H[i_flag_21] <- h
CO2[i_flag_21] <- K0[i_flag_21]*fCO2[i_flag_21]
HCO3[i_flag_21] <- K1[i_flag_21]*CO2[i_flag_21]/h
CO3[i_flag_21] <- K2[i_flag_21]*HCO3[i_flag_21]/h
DIC[i_flag_21] <- CO2[i_flag_21] + HCO3[i_flag_21] + CO3[i_flag_21]
# ------------ case 22.) HCO3 and pCO2 given ----
# Indices of flag elements where flag = 22
i_flag_22 <- which (flag == 22)
HCO3[i_flag_22] <- var2[i_flag_22]
CO2[i_flag_22] <- fCO2[i_flag_22]*K0[i_flag_22]
h <- CO2[i_flag_22]*K1[i_flag_22]/HCO3[i_flag_22]
CO3[i_flag_22] <- HCO3[i_flag_22]*K2[i_flag_22]/h
DIC[i_flag_22] <- CO2[i_flag_22] + HCO3[i_flag_22] + CO3[i_flag_22]
PH[i_flag_22] <- -log10(h)
H[i_flag_22] <- h
# ------------ case 23.) CO3 and pCO2 given ----
# Indices of flag elements where flag = 23
i_flag_23 <- which (flag == 23)
CO3[i_flag_23] <- var2[i_flag_23]
h <- sqrt(K0[i_flag_23]*K1[i_flag_23]*K2[i_flag_23]*fCO2[i_flag_23]/CO3[i_flag_23])
HCO3[i_flag_23] <- h*CO3[i_flag_23]/K2[i_flag_23]
CO2[i_flag_23] <- h*HCO3[i_flag_23]/K1[i_flag_23]
DIC[i_flag_23] <- CO2[i_flag_23] + HCO3[i_flag_23] + CO3[i_flag_23]
PH[i_flag_23] <- -log10(h)
H[i_flag_23] <- h
# ------------ case 24.) ALK and pCO2 given ----
# Indices of flag elements where flag = 24
i_flag_24 <- which (flag == 24)
ALK[i_flag_24] <- var2[i_flag_24]
CO2[i_flag_24] <- fCO2[i_flag_24]*K0[i_flag_24]
# From this line on, this case is similar to case 4
fALK <- function(x)# K1=K1, K2=K2, CO2=CO2, BOR=BOR, Kb=Kb, Kw=Kw, Pt=Pt, K1p=K1p, K2p=K2p, K3p=K3p, Sit=Sit, Ksi=Ksi, ST=ST, Ks=Ks, FLUO=FLUO, Kf=Kf, ALK=ALK)
{
# components for ALK
dic <- co2*(1+K1_i/x+K1_i*K2_i/(x*x))
hco3 <- dic*x*K1_i/(x*x + K1_i*x + K1_i*K2_i)
co3 <- dic*K1_i*K2_i/(x*x + K1_i*x + K1_i*K2_i)
boh4 <- bor/(1+x/Kb_i)
oh <- Kw_i/x
h3po4 <- Pt_i*x^3/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
hpo4 <- Pt_i*K1p_i*K2p_i*x/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
po4 <- Pt_i*K1p_i*K2p_i*K3p_i/(x^3+K1p_i*x^2+K1p_i*K2p_i*x+K1p_i*K2p_i*K3p_i)
siooh3 <- Sit_i/(1+x/Ksi_i)
## calculate Hfree and Htot
if(phs=="F")
{
hfree <- x ## if pHscale = free scale
htot <- x / (1+st/Ks_i)
}
else if(phs=="T")
{
hfree <- x * (1+st/Ks_i)
htot <- x
}
else if(phs=="SWS")
{
hfree <- x * (1 + st/Ks_i + fluo/Kff_i)
htot <- hfree / (1+st/Ks_i)
}
hso4 <- st/(1+Ks_i/hfree)
hf <- fluo/(1+Kf_i/htot)
############
OUT <- hco3+2*co3+boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4-alk
OUT
}
# Calculate [H+] from [CO2] and total alk
h <- rep(NA, length(i_flag_24))
j <- 1
for(i in (i_flag_24))
{
# Parameters used by function fAlk defined above and called below through uniroot()
K1_i <- K1[i]
K2_i <- K2[i]
K1p_i <- K1p[i]
K2p_i <- K2p[i]
K3p_i <- K3p[i]
Kb_i <- Kb[i]
Kw_i <- Kw[i]
Ksi_i <- Ksi[i]
Ks_i <- Ks[i]
Kf_i <- Kf[i]
Kff_i <- Kff[i]
Sit_i <- Sit[i]
Pt_i <- Pt[i]
co2 <- CO2[i]
alk <- ALK[i]
st <- ST[i]
bor <- BOR[i]
fluo <- FLUO[i]
phs <- pHscale[i]
# Calculate [H+] from total alk
h[j] <- uniroot(fALK,c(1e-10,10^(-3.5)), tol=1e-20)$root
j <- j + 1
}
HCO3[i_flag_24] <- K1[i_flag_24]*CO2[i_flag_24]/h
CO3[i_flag_24] <- K2[i_flag_24]*HCO3[i_flag_24]/h
PH[i_flag_24] <- -log10(h)
H[i_flag_24] <- h
DIC[i_flag_24] <- CO2[i_flag_24] + HCO3[i_flag_24] + CO3[i_flag_24]
# ------------ case 25.) DIC and pCO2 given ----
# Indices of flag elements where flag = 25
i_flag_25 <- which (flag == 25)
DIC[i_flag_25] <- var2[i_flag_25]
CO2[i_flag_25] <- K0[i_flag_25]*fCO2[i_flag_25]
# Though case 25 is the same as case 5, computations are made in a different way
K <- K1[i_flag_25]/K2[i_flag_25]
b <- K*K0[i_flag_25]*fCO2[i_flag_25]
c <- (K*K0[i_flag_25]*fCO2[i_flag_25]) *
(K0[i_flag_25]*fCO2[i_flag_25]-DIC[i_flag_25])
D <- b*b - 4*c
HCO3[i_flag_25] <- (1/2)*(-b + sqrt(D))
CO3[i_flag_25] <- DIC[i_flag_25] - CO2[i_flag_25] - HCO3[i_flag_25]
h <- K1[i_flag_25]*CO2[i_flag_25]/HCO3[i_flag_25]
PH[i_flag_25] <- -log10(h)
H[i_flag_25] <- h
# ------------ CALCULATION OF ALK in cases ----
cases <- c(1, 2, 3, 5, 6, 7, 9, 10, 12, 14, 21, 22, 23, 25)
# Indices of flag elements in these cases
i_flag <- which (flag %in% cases)
h <- H[i_flag]
# HCO3[i_flag] <- DIC[i_flag]*h*K1[i_flag]/(h*h + K1[i_flag]*h + K1[i_flag]*K2[i_flag])
# CO3[i_flag] <- HCO3[i_flag] * K2[i_flag] / h
boh4 <- BOR[i_flag]/(1+h/Kb[i_flag])
oh <- Kw[i_flag]/h
temp <- h^3 + K1p[i_flag]*h^2 + K1p[i_flag]*K2p[i_flag]*h + K1p[i_flag]*K2p[i_flag]*K3p[i_flag]
h3po4 <- Pt[i_flag]*(h^3) / temp
hpo4 <- Pt[i_flag]*K1p[i_flag]*K2p[i_flag]*h / temp
po4 <- Pt[i_flag]*K1p[i_flag]*K2p[i_flag]*K3p[i_flag] / temp
siooh3 <- Sit[i_flag]/(1+h/Ksi[i_flag])
## calculate Hfree anf Htot
hfree <- rep(NA, length(i_flag))
htot <- rep(NA, length(i_flag))
sc <- pHscale[i_flag]
st <- ST[i_flag]
ks <- Ks[i_flag]
fluo <- FLUO[i_flag]
kff <- Kff[i_flag]
# Where pHscale=="F", pHscale = free scale
i_sc_F <- which (sc == "F")
hfree[i_sc_F] <- h[i_sc_F]
htot[i_sc_F] <- h[i_sc_F] / (1+st[i_sc_F]/ks[i_sc_F])
# Where pHscale=="T", pHscale = total scale
i_sc_T <- which (sc == "T")
hfree[i_sc_T] <- h[i_sc_T] * (1+st[i_sc_T]/ks[i_sc_T])
htot[i_sc_T] <- h[i_sc_T]
# Where pHscale=="SWS", pHscale = SW scale
i_sc_S <- which (sc == "SWS")
hfree[i_sc_S] <- h[i_sc_S] * (1 + st[i_sc_S]/ks[i_sc_S] + fluo[i_sc_S]/kff[i_sc_S])
htot[i_sc_S] <- hfree[i_sc_S] / (1+st[i_sc_S]/ks[i_sc_S])
hso4 <- st/(1+ks/hfree)
hf <- fluo/(1+Kf[i_flag]/htot)
ALK[i_flag] <- HCO3[i_flag] + 2*CO3[i_flag] + boh4+oh+hpo4+2*po4+siooh3-hfree-hso4-hf-h3po4
##########################################################
# CALCULATION OF ARAGONITE AND CALCITE SATURATION STATE #
##########################################################
Ca = (0.02128/40.078) * SP/1.80655 #Improved formula from Dickson et al. (2007) as discussed in Orr & Epitalon (2014, revised)
Oa <- (Ca*CO3)/Kspa
Oc <- (Ca*CO3)/Kspc
#PCO2 and fCO2 - convert from atm to microatm
pCO2 <- pCO2*1e6
fCO2 <- fCO2*1e6
pCO2pot <- pCO2pot*1e6
fCO2pot <- fCO2pot*1e6
pCO2insitu <- pCO2insitu*1e6
fCO2insitu <- fCO2insitu*1e6
RES <- data.frame(flag, S, T, Patm, P, PH, CO2, fCO2, pCO2, fCO2pot, pCO2pot, fCO2insitu, pCO2insitu, HCO3, CO3, DIC, ALK, Oa, Oc)
names(RES) <- c("flag", "S", "T", "Patm", "P", "pH", "CO2", "fCO2", "pCO2", "fCO2pot", "pCO2pot", "fCO2insitu", "pCO2insitu", "HCO3", "CO3", "DIC", "ALK", "OmegaAragonite", "OmegaCalcite")
return(RES)
}K0x() is the equivalent to the function K0(), and returns Henry’s constant mol/(kg/atm).
"K0x" <-
function(S=35,T=25,P=0,Patm=1,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(Patm), length(ca), length(mg))
##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
if(length(Patm)!=nK){Patm <- rep(Patm[1], nK)}
#------- Constants ----------------
#---- issues of equic----
tk <- 273.15; # [K] (for conversion [deg C] <-> [K])
TK <- T + tk; # TC [C]; T[K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
#---------------------------------------------------------------------
#---------------------- K0 (K Henry) ---------------------------------
#
# CO2(g) <-> CO2(aq.)
# K0 = [CO2]/ p CO2
#
# Weiss (1974) [mol/kg/atm]
#
#
tmp = constants$K0p0 + constants$K0p1*100/TK + constants$K0p2*log(TK/100);
nK0we74x = tmp + S*(constants$K0p3 + constants$K0p4*TK/100 + constants$K0p5*(TK/100)^2);
##------------Warnings
is_w <- warn == "y"
if (any (is_w & (T>45 | S>45 | T<(-1) | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) ) {warning("mg, ca, S and/or T is outside the range of validity of the formulation available for K0x in seacarbx.")}
K0x = exp(nK0we74x);
#---------------------------------------------------------------------
#---------------------- Pressure correction ----------------------------
Phydro_atm = P / 1.01325 # convert hydrostatic pressure from bar to atm (1.01325 bar / atm)
Ptot = Patm + Phydro_atm # total pressure (in atm) = atmospheric pressure + hydrostatic pressure
R = 82.05736 # (cm3 * atm) / (mol * K) CODATA (2006)
vbarCO2 = 32.3 # partial molal volume (cm3 / mol) from Weiss (1974, Appendix, paragraph 3)
K0x = K0x * exp( ((1-Ptot)*vbarCO2)/(R*TK) ) # Weiss (1974, equation 5), TK is absolute temperature (K)
attr(K0x,"unit") = "mol/kg"
return(K0x)
}K1x() is the equivalent to the function K1(), and returns the first dissociation constant of carbonic acid (mol/kg).
"K1x" <-
function(S=35,T=25,P=0,pHscale="T",kSWS2scale=0,ktotal2SWS_P0=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(pHscale), length(kSWS2scale) ,length(ktotal2SWS_P0), length(ca), length(mg))
##-------- Create vectors for all input (if vectorial)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
if(length(pHscale)!=nK){pHscale <- rep(pHscale[1], nK)}
##---------- pHsc : this vector is not the actual pHscale because it can change during processing
pHsc <- rep(NA,nK)
pHlabel <- rep(NA,nK)
K1x <- rep(NA, nK)
#------- Constants ----------------
#---- issues of equic----
tk <- 273.15; # [K] (for conversion [deg C] <-> [K])
TK <- T + tk; # TC [C]; T[K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
# --------------------- K1x ---------------------------------------
# first acidity constant:
# [H^+] [HCO_3^-] / [CO2] = K_1
#
# Mehrbach et al (1973) refit by Lueker et al. (2000).
#
# (Lueker et al., 2000 in Guide to the Best Practices for Ocean CO2 Measurements
# Dickson, Sabine and Christian , 2007, Chapter 5, p. 13)
#
# pH-scale: 'total'. mol/kg-soln
# Effect of [Ca] and [Mg] is taken into account in the following formula (see Hain et al., 2014)
logK1xlue <- (constants$K1p0 + constants$K1p1/TK + constants$K1p2*log(TK) + constants$K1p3*S + constants$K1p4*(S^2))
K1x <- 10^logK1xlue
pHsc <- "T"
#-------------------- pH scale -------------------------------
##----------------- Conversion from total to SWS scale
## if pressure correction needed
## or pH scale conversion required anyway
## (in which case SWS may be an intermediate stage of conversion)
convert <- (pHsc == "T") & ((P > 0) | (pHscale != pHsc))
if (any (convert))
{
##------------- Convert from total to SWS scale
# if correction factor (from Total scale to seawater at P=0) not given
if (missing(ktotal2SWS_P0))
{
# Compute it
ktotal2SWS_P0 <- rep(1.0,nK)
ktotal2SWS_P0 <- kconv(S=S[convert], T=T[convert], P=0)$ktotal2SWS
}
else
{
# Check its length
if(length(ktotal2SWS_P0)!=nK) ktotal2SWS_P0 <- rep(ktotal2SWS_P0[1], nK)
# Filter
ktotal2SWS_P0 <- ktotal2SWS_P0[convert]
}
K1x[convert] <- K1x[convert] * ktotal2SWS_P0
pHsc[convert] <- "SWS"
# Just computed constant is on Free scale only if required scale is free scale
# and if no pressure correction needed
# --> No need to convert from free to other scale
# --> No need to determine conversion factor from free to SWS scale
}
# ------------------- Pressure effect --------------------------------
i_press <- which (P > 0)
if (length(i_press) > 0)
{
# Call Pcorrect() on SWS scale
K1x[i_press] <- Pcorrect(Kvalue=K1x[i_press], Ktype="K1", T=T[i_press],
S=S[i_press], P=P[i_press], pHscale=pHsc[i_press], 1., 1., warn=warn)
}
## --------------- Last conversion in the required pH scale ----------------
convert <- (pHscale != pHsc)
# At this stage, if any conversion is needed, it is from SWS scale
# Which pH scales are required ?
is_total <- convert & (pHscale=="T")
is_free <- convert & (pHscale=="F")
# if any pH scale correction required (from SWS scale)
if (any(convert))
{
# if pH scale correction factor not given
if (missing(kSWS2scale))
{
# Compute it
kSWS2scale <- rep(1.0,nK)
if (any(is_total)){ kSWS2scale[is_total] <- kconv(S=S[is_total], T=T[is_total], P=P[is_total])$kSWS2total }
if (any(is_free)){ kSWS2scale[is_free] <- kconv(S=S[is_free], T=T[is_free], P=P[is_free])$kSWS2free }
}
else
{
# Check its length
if (length(kSWS2scale)!=nK) kSWS2scale <- rep(kSWS2scale[1], nK)
}
# Apply pH scale correction
K1x[convert] <- K1x[convert] * kSWS2scale[convert]
}
# Return full name of pH scale
is_total <- (pHscale=="T")
is_free <- (pHscale=="F")
pHlabel[is_total] <- "total scale"
pHlabel[is_free] <- "free scale"
pHlabel[!is_total & !is_free] <- "seawater scale"
##------------ Warnings --------------
if (any(T>35 | T<2 | S<19 | S>43 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) {
warning("mg, ca, S and/or T is outside the range of validity of the formulation chosen for K1x.")}
##--------------- Attributes ------------
method <- rep(NA, nK)
method <- "Luecker et al. (2000)"
attr(K1x,"unit") <- "mol/kg-soln"
attr(K1x,"pH scale") <- pHlabel
attr(K1x, "method") <- method
return(K1x)
}K2x() is the equivalent to the function K2(), and returns the second dissociation constant of carbonic acid (mol/kg).
"K2x" <-
function(S=35,T=25,P=0,pHscale="T",kSWS2scale=0,ktotal2SWS_P0=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(pHscale), length(kSWS2scale) ,length(ktotal2SWS_P0), length(ca), length(mg))
##-------- Create vectors for all input (if vectorial)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
if(length(pHscale)!=nK){pHscale <- rep(pHscale[1], nK)}
##---------- Initialisation of vectors
pHsc <- rep(NA,nK) # pHsc : this vector note the actual pHscale because it can change during processing
pHlabel <- rep(NA,nK)
K2x <- rep(NA, nK)
#-------Constants----------------
#---- issues of equic----
tk <- 273.15; # [K] (for conversion [deg C] <-> [K])
TK <- T + tk; # T [C]; T [K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
# --------------------- K2x ------------------------------
# second acidity constant:
# [H^+] [CO_3^--] / [HCO_3^-] = K_2
#
# Mehrbach et al. (1973) refit by Lueker et al. (2000).
#
# (Lueker et al., 2000 in Guide to the Best Practices for Ocean CO2 Measurements
# Dickson, Sabin and Christian , 2007, Chapter 5, p. 14)
#
# pH-scale: 'total'. mol/kg-soln
# Effect of [Ca] and [Mg] is taken into account in the following formula (see Hain et al., 2014)
logK2xlue <- (constants$K2p0 + constants$K2p1/TK + constants$K2p2*log(TK) + constants$K2p3*S + constants$K2p4*(S^2))
K2x <- 10^logK2xlue
pHsc <- "T"
#-------------------- pH scale -------------------------------
##----------------- Conversion from total to SWS scale
## if pressure correction needed
## or pH scale conversion required anyway
## (in which case SWS may be an intermediate stage of conversion)
convert <- (pHsc == "T") & ((P > 0) | (pHscale != pHsc))
if (any (convert))
{
##------------- Convert from total to SWS scale
# if correction factor (from Total scale to seawater at P=0) not given
if (missing(ktotal2SWS_P0))
{
# Compute it
ktotal2SWS_P0 <- rep(1.0,nK)
ktotal2SWS_P0 <- kconv(S=S[convert], T=T[convert], P=0)$ktotal2SWS
}
else
{
# Check its length
if(length(ktotal2SWS_P0)!=nK) ktotal2SWS_P0 <- rep(ktotal2SWS_P0[1], nK)
# Filter
ktotal2SWS_P0 <- ktotal2SWS_P0[convert]
}
K2x[convert] <- K2x[convert] * ktotal2SWS_P0
pHsc[convert] <- "SWS"
# Just computed constant is on Free scale only if required scale is free scale
# and if no pressure correction needed
# --> No need to convert from free to other scale
# --> No need to determine conversion factor from free to SWS scale
}
# ------------------- Pressure effect --------------------------------
i_press <- which (P > 0)
if (length(i_press) > 0)
{
# Call Pcorrect() on SWS scale
K2x[i_press] <- Pcorrect(Kvalue=K2x[i_press], Ktype="K2", T=T[i_press],
S=S[i_press], P=P[i_press], pHscale=pHsc[i_press], 1., 1., warn=warn)
}
## ------------------- Last conversion in the require pHscale ------------------
convert <- (pHscale != pHsc)
# At this stage, if any conversion is needed, it is from SWS scale
# Which pH scales are required ?
is_total <- convert & (pHscale=="T")
is_free <- convert & (pHscale=="F")
# if any pH scale correction required (from SWS scale)
if (any(convert))
{
# if pH scale correction factor not given
if (missing(kSWS2scale))
{
# Compute it
kSWS2scale <- rep(1.0,nK)
if (any(is_total)){ kSWS2scale[is_total] <- kconv(S=S[is_total], T=T[is_total], P=P[is_total])$kSWS2total }
if (any(is_free)){ kSWS2scale[is_free] <- kconv(S=S[is_free], T=T[is_free], P=P[is_free])$kSWS2free }
}
else
# Check its length
if (length(kSWS2scale)!=nK) kSWS2scale <- rep(kSWS2scale[1], nK)
# Apply pH scale correction
K2x[convert] <- K2x[convert] * kSWS2scale[convert]
}
# Return full name of pH scale
is_total <- (pHscale=="T")
is_free <- (pHscale=="F")
pHlabel[is_total] <- "total scale"
pHlabel[is_free] <- "free scale"
pHlabel[!is_total & !is_free] <- "seawater scale"
##------------ Warnings --------------
if (any(T>35 | T<2 | S<19 | S>43 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) {
warning("mg, ca, S and/or T is outside the range of validity of the formulation chosen for K2x.")}
##---------------Attributes
method <- rep(NA, nK)
method <- "Luecker et al. (2000)"
attr(K2x,"unit") = "mol/kg-soln"
attr(K2x,"pH scale") = pHlabel
attr(K2x,"method") = method
return(K2x)
}Ksx() is the equivalent to the function Ks(), and returns the stability constant of hydrogen sulfate (mol/kg).
"Ksx" <-
function(S=35,T=25,P=0, ks="d",warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(ks), length(ca), length(mg))
##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ks)!=nK){ks <- rep(ks[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
#-------Constants----------------
#---- issues de equic----
tk = 273.15; # [K] (for conversion [deg C] <-> [K])
TK = T + tk; # T [C]; TK [K]
Cl = S / 1.80655; # Cl = chlorinity; S = salinity (per mille)
iom0 = 19.924*S/(1000-1.005*S);
ST = 0.14 * Cl/96.062 # (mol/kg) total sulfate (Dickson et al., 2007, Table 2)
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
#--------------------------------------------------------------
#------------------ Ks ----------------------------------------
# Dickson (1990)
#
# Equilibrium constant for HSO4- = H+ + SO4--
#
# K_S = [H+]free [SO4--] / [HSO4-]
# pH-scale: free scale !!!
#
# the term log(1-0.001005*S) converts from
# mol/kg H2O to mol/kg soln
tmp1 = constants$KSO4p1 / TK + constants$KSO4p0 + constants$KSO4p2*log(TK);
tmp2 = +(constants$KSO4p3 / TK + constants$KSO4p4 + constants$KSO4p5 * log(TK))*sqrt(iom0);
tmp3 = +(constants$KSO4p6 / TK + constants$KSO4p7 + constants$KSO4p8 * log(TK))*iom0;
tmp4 = constants$KSO4p9 / TK *iom0^1.5 + constants$KSO4p10 / TK *iom0^2;
lnKsx = tmp1 + tmp2 + tmp3 + tmp4 + log(1-0.001005*S);
Ks_d = exp(lnKsx);
#--------------------------------------------------------------
#------------------ Ks ----------------------------------------
# Khoo et al. 1977
#
# Equilibrium constant for HSO4- = H+ + SO4--
#
# K_S = [H+]free [SO4--] / [HSO4-]
# pH-scale: free scale !!!
#
# correct for T range : 5 - 40°C
# correct for S range : 20 - 45
I35 <- 0.7227
I <- I35*(27.570*S)/(1000-1.0016*S) #formal ionic strength
logbeta <- 647.59/TK - 6.3451 + 0.019085*TK - 0.5208*I^(0.5)
beta <- 10^(logbeta)
Ks_kx <- 1/beta
#-------------------------------------------------------------------
#--------------- choice between the formulations -------------------
is_k <- (ks=='k') # everything that is not "k" is "d" by default
Ksx <- Ks_d
Ksx[is_k] <- Ks_kx[is_k]
method <- rep(NA, nK)
method[!is_k] <- "Dickson (1990)"
method[is_k] <- "Khoo et al. (1977)"
# ------------------- Pressure effect --------------------------------
if (any(P != 0))
Ksx <- Pcorrect(Kvalue=Ksx, Ktype="Ks", T=T, S=S, P=P, pHscale="F", warn=warn)
##------------Warnings
is_w <- warn == "y"
if (any (is_w & is_k & (T<5 | T>40 | S<20 | S>45 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)))
{warning("S and/or T is outside the range of validity of the formulation chosen for Ks.")}
if (any (is_w & (T>45 | S>45 | T<0 | S<5 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)))
{warning("mg, ca, S and/or T is outside the range of validity of the formulations available for Ksx in seacarbx.")}
##------------Attributes
attr(Ksx,"unit") = "mol/kg-soln"
attr(Ksx,"pH scale") = "free scale"
attr(Ksx, "method") = method
return(Ksx)
}Kwx() is the equivalent to the function Kw(), and returns the ion product of water (mol2/kg2)
"Kwx" <-
function(S=35,T=25,P=0,pHscale="T",kSWS2scale=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(pHscale), length(kSWS2scale), length(ca), length(mg))
##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
if(length(pHscale)!=nK){pHscale <- rep(pHscale[1], nK)}
#-------Constantes----------------
#---- issues de equic----
tk = 273.15; # [K] (for conversion [deg C] <-> [K])
TK = T + tk; # T [C]; TK[K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
#-------------------------------------------------------------------
# --------------------- Kwater -------------------------------------
#
# Millero (1995)(in Guide to Best Practices in Ocean CO2 Measurements (2007, Chapter 5, p.16))
# $K_w$ in mol/kg-soln.
# pH-scale: pH$_{total}$ ('total` scale).
# *** J. Orr (15 Jan 2013): Formulation changed to be on the SWS scale (without later conversion)
#tmp1 = -13847.26/TK + 148.9652 - 23.6521 * log(TK); ## second term is 148.9652 as given in Dickson (2007) instead of 148.96502 as recommanded in DOE (1994).
# From J. C. Orr on 15 Jan 2013:
# The formulation above was a modified version of Millero (1995) where Dickson et al. (2007) subtracted 0.015
# from Millero's original constant (148.9802) to give 148.9652 (the 2nd term above). BUT Dickson's reason for that
# operation was to "convert--approximately--from theSWS pH scale (including HF) used by Millero (1995) to the 'total'
# scale ...".
# This subtraction of 0.015 to switch from the SWS to Total scale is not good for 2 reasons:
# (1) The 0.015 value is inexact (not constant), e.g., it is 0.022 at T=25, S=35, P=0;
# (2) It makes no sense to switch to the Total scale when just below you switch back to the SWS scale.
# The best solution is to reestablish the original equation (SWS scale) and delete the subsequent scale conversion.
tmp1 = constants$Kwp1/TK + constants$Kwp0 + constants$Kwp2*log(TK); # now the original formulation: Millero (1995)
tmp2 = + (constants$Kwp3/TK + constants$Kwp4 + constants$Kwp5*log(TK))*sqrt(S) + constants$Kwp6*S;
lnKwx = tmp1 + tmp2;
Kwx = exp(lnKwx);
# ---- Conversion from Total scale to seawater scale before pressure corrections
# *** JCO: This is no longer necessary: with original formulation (Millero, 1995), Kw is on "seawater scale"!
#factor <- kconv(S=S, T=T, P=rep(0,nK))$ktotal2SWS
#Kw <- Kw * factor
# If needed, pressure correction
if (any (P != 0))
Kwx <- Pcorrect(Kvalue=Kwx, Ktype="Kw", T=T, S=S, P=P, pHscale="SWS", 1., 1., warn=warn)
###----------------pH scale corrections
# Which pH scales are required ?
is_total <- pHscale=="T"
is_free <- pHscale=="F"
# if any pH scale correction required (from SWS scale)
if (any(is_total) || any(is_free))
{
# if pH scale correction factor not given
if (missing(kSWS2scale))
{
# Compute it
kSWS2scale <- rep(1.0,nK)
if (any(is_total))
kSWS2scale[is_total] <- kconv(S=S[is_total], T=T[is_total], P=P[is_total])$kSWS2total
if (any(is_free))
kSWS2scale[is_free] <- kconv(S=S[is_free], T=T[is_free], P=P[is_free])$kSWS2free
}
else
# Check its length
if(length(kSWS2scale)!=nK){kSWS2scale <- rep(kSWS2scale[1], nK)}
# Apply pH scale correction
Kwx <- Kwx*kSWS2scale
}
# Return full name of pH scale
pHsc <- rep(NA,nK)
pHsc[is_total] <- "total scale"
pHsc[is_free] <- "free scale"
pHsc[!is_total & !is_free] <- "seawater scale"
##------------Warnings
is_w <- warn == "y"
if (any (is_w & (T>45 | S>45 | T<0 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) ) {warning("mg, ca, S and/or T is outside the range of validity of the formulation available for Kwx in seacarbx.")}
attr(Kwx,"unit") = "mol/kg-soln"
attr(Kwx,"pH scale") = pHsc
return(Kwx)
}Kbx() is the equivalent to the function Kb(), and returns the dissociation constant of boric acid (mol/kg)
"Kbx" <-
function(S=35,T=25,P=0,pHscale="T",kSWS2scale=0,ktotal2SWS_P0=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(pHscale), length(kSWS2scale) ,length(ktotal2SWS_P0), length(ca), length(mg))
##-------- Create vectors for all input (if vectorial)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
if(length(pHscale)!=nK){pHscale <- rep(pHscale[1], nK)}
#-------Constants----------------
#---- issues of equic----
tk <- 273.15; # [K] (for conversion [deg C] <-> [K])
TK <- T + tk; # T [C]; TK [K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
#---------------------------------------------------------------------
# --------------------- Kb --------------------------------------------
# Kbor = [H+][B(OH)4-]/[B(OH)3]
#
# (Dickson, 1990 in Guide to Best Practices in Ocean CO2 Measurements 2007)
# pH-scale: 'total'. mol/kg-soln
tmp1 = (constants$KBp3 + constants$KBp4*sqrt(S) + constants$KBp5*S + constants$KBp6*S^(3/2) + constants$KBp7*S^2);
tmp2 = constants$KBp0 + constants$KBp1*sqrt(S) + constants$KBp2*S;
tmp3 = +(constants$KBp8 + constants$KBp9*sqrt(S) + constants$KBp10*S)*log(TK);
lnKbx = tmp1 / TK + tmp2 + tmp3 + constants$KBp11*sqrt(S)*TK;
Kbx <- exp(lnKbx)
## ---- Conversion from Total scale to seawater scale before pressure corrections
# if correction factor (from Total scale to seawater at P=0) not given
if (missing(ktotal2SWS_P0))
{
# Compute it
ktotal2SWS_P0 <- kconv(S=S, T=T, P=0)$ktotal2SWS
}
else
# Check its length
if(length(ktotal2SWS_P0)!=nK) ktotal2SWS_P0 <- rep(ktotal2SWS_P0[1], nK)
Kbx <- Kbx * ktotal2SWS_P0
# -------------------Correct for Pressure effect, if needed ---------------------------
if (any (P != 0))
Kbx <- Pcorrect(Kvalue=Kbx, Ktype="Kb", T=T, S=S, P=P, pHscale="SWS", 1., 1., warn=warn)
###----------------pH scale corrections
# Which pH scales are required ?
is_total <- pHscale=="T"
is_free <- pHscale=="F"
# if any pH scale correction required (from SWS scale)
if (any(is_total) || any(is_free))
{
# if pH scale correction factor not given
if (missing(kSWS2scale))
{
# Compute it
kSWS2scale <- rep(1.0,nK)
if (any(is_total)) {kSWS2scale[is_total] <- kconv(S=S[is_total], T=T[is_total], P=P[is_total])$kSWS2total}
if (any(is_free)) {kSWS2scale[is_free] <- kconv(S=S[is_free], T=T[is_free], P=P[is_free])$kSWS2free}
}
else
# Check its length
if(length(kSWS2scale)!=nK){kSWS2scale <- rep(kSWS2scale[1], nK)}
# Apply pH scale correction
Kbx <- Kbx*kSWS2scale
}
# Return full name of pH scale
pHsc <- rep(NA,nK)
pHsc[is_total] <- "total scale"
pHsc[is_free] <- "free scale"
pHsc[!is_total & !is_free] <- "seawater scale"
##------------Warnings
is_w <- warn == "y"
if (any (is_w & (T>45 | S>45 | T<0 | S<5 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) ) {warning("mg, ca, S and/or T is outside the range of validity of the formulation available for Kbx in seacarbx.")}
attr(Kbx,"unit") = "mol/kg-soln"
attr(Kbx,"pH scale") = pHsc
return(Kbx)
}Kspax() is the equivalent to the function Kspa(), and returns the solubility product of aragonite (mol/kg)
"Kspax" <-
function(S=35,T=25,P=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(ca), length(mg))
##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
#---- issues de equic----
tk = 273.15; # [K] (for conversion [deg C] <-> [K])
TC = T + tk; # TC [C]; T[K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
# --------------------- Kspa (aragonite) ----------------------------
#
# apparent solubility product of aragonite
#
# Kspa = [Ca2+]T [CO32-]T
#
# where $[]_T$ refers to the equilibrium total
# (free + complexed) ion concentration.
#
# Mucci 1983 mol/kg-soln
tmp1 = constants$KspAp0 + constants$KspAp1*TC + constants$KspAp2/TC + constants$KspAp3*log10(TC);
tmp2 = +(constants$KspAp4 + constants$KspAp5*TC + constants$KspAp6/TC)*sqrt(S);
tmp3 = constants$KspAp7*S + constants$KspAp8*S^1.5;
log10Kspax = tmp1 + tmp2 + tmp3;
Kspax = 10^(log10Kspax);
# ----------------- Pressure Correction ------------------
Kspax <- Pcorrect(Kvalue=Kspax, Ktype="Kspa", T=T, S=S, P=P, warn=warn)
##------------Warnings
is_w <- warn == "y"
if (any (is_w & (T>40 | S>44 | T<5 | S<5 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) ){
warning("mg, ca, S and/or T is outside the range of validity of the formulation available for Kspax in seacarbx.")}
attr(Kspax,"unit") = "mol/kg"
return(Kspax)
}Kspcx() is the equivalent to the function Kspc(), and returns the solubility product of calcite (mol/kg)
"Kspcx" <-
function(S=35,T=25,P=0,warn="y",ca=10.2821e-03,mg=52.8171e-03) {
nK <- max(length(S), length(T), length(P), length(ca), length(mg))
##-------- Creation de vecteur pour toutes les entrees (si vectorielles)
if(length(S)!=nK){S <- rep(S[1], nK)}
if(length(T)!=nK){T <- rep(T[1], nK)}
if(length(P)!=nK){P <- rep(P[1], nK)}
if(length(ca)!=nK){ca <- rep(ca[1], nK)}
if(length(mg)!=nK){mg <- rep(mg[1], nK)}
#---- issues de equic----
tk = 273.15; # [K] (for conversion [deg C] <-> [K])
TC = T + tk; # TC [C]; T[K]
#------- MyAMI constants ------------
constants <- if (!exists("constants")) {
lookup_fn(ca=ca, mg=mg)
} else if (ca < 0.001 || mg < 0.001 || ca > 0.06 || mg > 0.06) {
constants[0]
} else if (ca == round(constants$ca, digits=8) && mg == round(constants$mg, digits=8)) {
constants
} else {
lookup_fn(ca=ca, mg=mg)
}
# --------------------- Kspc (calcite) ----------------------------
#
# apparent solubility product of calcite
#
# Kspc = [Ca2+]T [CO32-]T
#
# where $[]_T$ refers to the equilibrium total
# (free + complexed) ion concentration.
#
# Mucci 1983 mol/kg-soln
tmp1 = constants$KspCp0 + constants$KspCp1*TC + constants$KspCp2/TC + constants$KspCp3*log10(TC);
tmp2 = +(constants$KspCp4 + constants$KspCp5*TC + constants$KspCp6/TC)*sqrt(S);
tmp3 = constants$KspCp7*S + constants$KspCp8*S^1.5;
log10Kspcx = tmp1 + tmp2 + tmp3;
Kspcx = 10^(log10Kspcx);
# ----------------- Pressure Correction ------------------
Kspcx <- Pcorrect(Kvalue=Kspcx, Ktype="Kspc", T=T, S=S, P=P, warn=warn)
##------------Warnings
is_w <- warn == "y"
if (any (is_w & (T>40 | S>44 | T<5 | S<5 | mg<1e-03 | mg>60e-03 | ca<1e-03 | ca>60e-03)) ){
warning("mg, ca, S and/or T is outside the range of validity of the formulation available for Kspcx in seacarbx.")}
attr(Kspcx,"unit") = "mol/kg"
return(Kspcx)
}The usage of seacarbx is very simple, with a code written in R base and hence independent of further packages. Download seacarbx from here (file size 1.8 MB), from GitHub or from Zenodo, unzip it to your working directory, and load all functions and the tabulated data from within R with:
load("~/path-to-your-directory/seacarbx.rda")NOTE that seacarbx only works when the original seacarb package is loaded as well:
library(seacarb)The syntax is basically the same as for the original functions, except the additional input parameters [Mg2+] and [Ca2+] (in mol/L), which are, as an example, called like this:
carbx(flag=1, var1=8.1, var2=12e-6, S=35, T=25, ca=10.2821e-03, mg=52.8171e-03,...)Every time that the input variables [Ca2+] and/or [Mg2+] are changed in the seacarbx functions, the bilinear interpolation function lookup_fn is called, which computes accurate parameters from the MyAMI lookup table (see also figure).
lookup_fn <- function(ca, mg, print=FALSE) {
cas <- 1e-12 + ca
mgs <- 1e-12 + mg
# First look for the closest ca match in lookup table
df <- myami[which(abs(myami$ca-cas) == min(abs(myami$ca-cas))),]
# Then check whether search ca value is smaller or larger than next table value. Depending on this, look for the next lower or next larger ca value and row-bind with first table
df2 <- if (df$ca[1] > cas || cas == 0.06 + 1e-12) {
rbind(df, myami[which(abs(myami$ca-cas+0.001) == min(abs(myami$ca-cas+0.001))),])
} else {
rbind(df, myami[which(abs(myami$ca-cas-0.001) == min(abs(myami$ca-cas-0.001))),])
}
# Then the same with mg
df3 <- df2[which(abs(df2$mg-mgs) == min(abs(df2$mg-mgs))),]
df4 <- if (df3$mg[1] > mgs || mgs == 0.06 + 1e-12) {
rbind(df3, df2[which(abs(df2$mg-mgs+0.001) == min(abs(df2$mg-mgs+0.001))),])
} else {
rbind(df3, df2[which(abs(df2$mg-mgs-0.001) == min(abs(df2$mg-mgs-0.001))),])
}
# Find parameters for 4 corner points surrounding target point P(ca,mg)
Q11 <- df4[which(df4$ca == min(df4$ca) & df4$mg == min(df4$mg)),] # lower ca, lower mg
Q21 <- df4[which(df4$ca == max(df4$ca) & df4$mg == min(df4$mg)),] # upper ca, lower mg
Q12 <- df4[which(df4$ca == min(df4$ca) & df4$mg == max(df4$mg)),] # lower ca, upper mg
Q22 <- df4[which(df4$ca == max(df4$ca) & df4$mg == max(df4$mg)),] # upper ca, upper mg
# Boundaries
ca1 <- Q11$ca # lower ca boundary
ca2 <- Q21$ca # upper ca boundary
mg1 <- Q11$mg # lower mg boundary
mg2 <- Q12$mg # upper mg boundary
# Parameters at auxiliary points fR
fR1 <- (ca2-ca)/(ca2-ca1)*Q11 + (ca-ca1)/(ca2-ca1)*Q21 # interpolated parameter value between ca1 and ca2 at mg1 (i.e. at low mg)
fR1
fR2 <- (ca2-ca)/(ca2-ca1)*Q12 + (ca-ca1)/(ca2-ca1)*Q22 # interpolated parameter value between ca1 and ca2 at mg2 (i.e. at high mg)
fR2
# Interpolated parameters at target point P(ca,mg)
fP <- (mg2-mg)/(mg2-mg1)*fR1 + (mg-mg1)/(mg2-mg1)*fR2 # interpolated parameter value between R1 and R2 (i.e. along y(=mg) axis)
fP
# write table into txt file, if required
if (print==TRUE) {
write.table(fP, file = "MyAMI_param_output.txt", sep = "\t",
row.names = FALSE)
}
# save dataframe into working directory
assign("constants", fP, envir = .GlobalEnv)
return(fP)
}lookup_fn can also be used as a separate lookup tool to calculate the MyAMI parameters for a given set of [Mg2+] and [Ca2+], which can optionally be saved as a text file to your working directory by using the argument print:
lookup_fn(ca=10.2821e-03, mg=52.8171e-03, print=TRUE)
MyAMI model results in e.g. ~2.2 µmol/kg steps for [CO32-] if tabulated parameters are used (left), while seacarbx applies a bilinear interpolation function to compute parameters that continuously vary with [Ca2+,Mg2+] (right).
Lab & software skills
Laboratory
Inductively Coupled Plasma Mass Spectrometry (ICP-MS), both Single- and Multicollector
Nano- and femtosecond laser ablation
Thermal Ionization Mass Spectrometry (TIMS)
Scanning Electron Microscopy (SEM)
Confocal laser scanning microscopy (CLSM)
X-ray diffractometry (XRD) analyses of clay and bulk samples
Culturing benthic and planktic foraminifera under controlled laboratory conditions
Software & technologies
Programming: R, Python, SQL
Big Data Analytics: MySQL workbench, Redis, ML, Matplotlib, Cloudera, MapReduce, Hadoop, Spark
Data engineering: BI, SQL, NoSQL, MS Access, SQLite, Data Warehouse modeling, ETL
Cloud computing; Microsoft Azure
Mapping: GMT, ODV, leaflet, sf
Graphics and Layout: Adobe CS, Gimp, Inkscape, Scribus, drawio
Office: MS Office, OpenOffice, LibreOffice
Others: VBA, Analyseries, OriginPro
OS: Linux, Windows, MAC
Résumé
| When | Where | With whom |
|---|---|---|
| since 2022 | Dettmer Group KG, Data Engineer/BI Developer | |
| 2021 – 2022 | Advanced training as Big Data Engineer | alfatraining |
| 2019 – 2021 | MARUM, University of Bremen | Prof. Heiko Pälike |
| 2019 | AWI: Isotope lab work (50%) and Program Management (50%) | Prof. Jelle Bijma, Martina Wilde |
| 2017 – 2019 | LU Hannover, Institut für Mineralogie | PD Ingo Horn |
| 2016 – 2017 | AWI | Prof. Jelle Bijma |
| 2014 – 2016 | MARUM, University of Bremen | Prof. Michal Kučera |
| 2011 – 2014 | AWI | Prof. Jelle Bijma |
| 2011 | University of Bremen | Prof. Michael Schulz |
| 2009 – 2010 | LDEO, Columbia University (USA) | Prof. Bärbel Hönisch |
| 2008 – 2009 | MARUM, University of Bremen | Dr. Torsten Bickert |
| 2007 – 2008 | Utrecht University (NL) | Prof. Gert-Jan Reichart |
| 2005 – 2008 | MARUM, University of Bremen: PhD studies | Prof. Gerold Wefer |
| 1995 – 2004 | FU Berlin: Studies of Geology and Paleontology (Diplom) | Prof. Christoph Heubeck |
Publications
Boer, W., Nordstad, S., Weber, M., Mertz-Kraus, R., Hönisch, B., Bijma, J., Raitzsch, M., Wilhelms-Dick, D., Foster, G. L., Goring-Harford, H., Nürnberg, D., Hauff, F., Kuhnert, H., Lugli, F., Spero, H., Rosner, M., van Gaever, P., de Nooijer, L. J., and Reichart, G.-J. (2022): A New Calcium Carbonate Nano-Particulate Pressed Powder Pellet (NFHS-2-NP) for LA-ICP-OES, LA-(MC)-ICP-MS and µXRF. – Geostandards and Geoanalytical Research, https://doi.org/10.1111/ggr.12425.
Raitzsch, M., Bijma, J., Bickert, T., Schulz, M., Holbourn, A., and Kučera, M. (2021): Atmospheric carbon dioxide variations across the middle Miocene climate transition. – Climate of the Past, 17(2): 703-719, https://doi.org/10.5194/cp-17-703-2021.
Raitzsch, M., Hain, M., Henehan, M., and Gattuso, J.-P. (2021): seacarbx - seacarb extension for deep-time carbonate system calculations. – Zenodo, http://doi.org/10.5281/zenodo.4432146.
van Dijk, I., Raitzsch, M., Brummer, G.-J., and Bijma, J. (2020): Novel method to image and quantify cogwheel structures in foraminiferal shells. – Frontiers in Ecology and Evolution, 8: https://doi.org/10.3389/fevo.2020.567231.
Raitzsch, M., Rollion-Bard, C., Horn, I., Steinhoefel, G., Benthien, A., Richter, K.-U., Buisson, M., Louvat, P., and Bijma, J. (2020): Technical Note: Single-shell δ11B analysis of Cibicidoides wuellerstorfi using femtosecond laser ablation MC-ICPMS and secondary ion mass spectrometry. – Biogeosciences, 17(21): 5365-5375, https://doi.org/10.5194/bg-17-5365-2020.
Zamelczyk, K., Rasmussen, T. L., Raitzsch, M., and Chierici, M. (2020): The last two millennia: climate, ocean circulation and paleoproductivity inferred from planktic foraminifera, southwestern Svalbard margin. – Polar Research, 39: https://doi.org/10.33265/polar.v39.3715.
Raitzsch, M. and Gattuso, J.-P. (2020): ScarFace - seacarb calculations with R Shiny user interface. – Zenodo, https://doi.org/10.5281/zenodo.3662139.
Tyszka, J., Bickmeyer, U., Raitzsch, M., Bijma, J., Höher, N., Topa, P., Kaczmarek, K., and Mewes, A. (2019): Form and function of F-actin during biomineralization revealed from live experiments on foraminifera. – PNAS, 116(10): 4111-4116, https://doi.org/10.1073/pnas.1810394116.
Raitzsch, M; Bijma, J., Benthien, A., Richter, K.-U., Steinhoefel, G., and Kučera, M. (2018): Boron isotope-based seasonal paleo-pH reconstruction for the Southeast Atlantic – A Multispecies approach using habitat preference of planktonic foraminifera. – Earth and Planetary Science Letters, 487: 138-150, https://doi.org/10.1016/j.epsl.2018.02.002.
Nakajima, K., Suzuki, M., Nagai, Y., Oaki, Y., Toyofuku, T., Bijma, J., Nehrke, G., Raitzsch, M., Tani, K., and Imai, H. (2017): Hierarchical textures on aragonitic shells of the hyaline radial foraminifer Hoeglundina elegans. – CrystEngComm, 19(47): 7191-7196, https://doi.org/10.1039/C7CE01870C.
Howes, E.L., Kaczmarek, K., Raitzsch, M., Mewes, A., Bijma, N., Horn, I., Misra, S., Gattuso, G.-P., and Bijma, J. (2017): Decoupled carbonate chemistry controls on the incorporation of boron into Orbulina universa. – Biogeosciences, 14: 415–430, https://doi.org/10.5194/bg-14-415-2017.
Wollenburg, J., Raitzsch, M., and Tiedemann, R. (2015): Novel high-pressure culture experiments on deep-sea benthic foraminifera – Evidence for methane seepage-related δ13C of Cibicides wuellerstorfi. – Marine Micropaleontology, 117: 47-64, https://doi.org/10.1016/j.marmicro.2015.04.003.
Milker, Y., Rachmayani, R., Weinkauf, M.F.G., Prange, M., Raitzsch, M., Schulz, M., and Kučera, M. (2015): Global Synthesis of Sea-Surface Temperature Trends During Marine Isotope Stage 11. – In M. Schulz and A. Paul (eds.): Integrated Analysis of Interglacial Climate Dynamics (INTERDYNAMIC), SpringerBriefs in Earth System Sciences, 139 pp, ISBN 978-3-319-00692-5.
Raitzsch, M. and Hönisch, B. (2013): Cenozoic boron isotope variations in benthic foraminifers. – Geology, 41(5): 591-594, https://doi.org/10.1130/G34031.1.
Milker, Y., Rachmayani, R., Weinkauf, M., Prange, M., Raitzsch, M., Schulz,M., and Kučera, M. (2013): Global and regional sea surface temperature trends during Marine Isotope Stage 11. – Climate of the Past, 9(5): 2231-2252, https://doi.org/10.5194/cp-9-2231-2013.
Raitzsch, M., Hathorne, E.C., Kuhnert, H., Groeneveld, J., and Bickert, T. (2011): Modern and Late Pleistocene B/Ca ratios of the benthic foraminifer Planulina wuellerstorfi determined with laser ablation ICP-MS. – Geology, 39(11): 1039-1042, https://doi.org/10.1130/G32009.1.
Dueñas-Bohórquez, A., Raitzsch, M., De Nooijer, L.J., and Reichart, G.-J. (2011): Independent impacts of calcium and carbonate ion concentration on Mg and Sr incorporation in cultured benthic foraminifera. – Marine Micropalentology, 81(3-4): 122-130, https://doi.org/10.1016/j.marmicro.2011.08.002.
Raitzsch, M., Kuhnert, H., Hathorne, E.C., Groeneveld, J., and Bickert, T. (2011): U/Ca in benthic foraminifers: A proxy for the deep-sea carbonate saturation. – Geochemistry Geophysics Geosystems, 12(6), Q06019, https://doi.org/10.1029/2010GC003344.
van Raden, U.J., Groeneveld, J., Raitzsch, M., and Kučera, M. (2011): Mg/Ca in the planktonic foraminifera Globorotalia inflata and Globigerinoides bulloides from Western Mediterranean plankton tow and core top samples. – Marine Micropaleontology, 78(3-4): 101-112, https://doi.org/10.1016/j.marmicro.2010.11.002.
Heubeck, C., Raitzsch, M., Völker, D., and Wiedicke-Hombach, M. (2010): Sedimentación. – In J. Díaz-Naveas and J. Frutos (eds.): Geología Marina de Chile. – Comité Oceanográfico Nacional de Chile - Pontificia Universidad Católica de Valparaíso - Servicio Nacional de Geología y Minería de Chile, 117 pp, ISBN 978-956-235-026-6.
Raitzsch, M., Dueñas-Bohórquez, A., Reichart, G.-J., De Nooijer, L.J., and Bickert, T. (2010): Incorporation of Mg and Sr in calcite of cultured benthic foraminifera: impact of calcium concentration and associated calcite saturation state. – Biogeosciences, 7: 869-881, https://doi.org/10.5194/bg-7-869-2010.
Raitzsch, M., Kuhnert, H., Groeneveld, J., and Bickert, T. (2008): Benthic foraminifer Mg/Ca anomalies in South Atlantic core top sediments and their implications for paleothermometry. – Geochemistry Geophysics Geosystems, 9(5), Q05010, https://doi.org/10.1029/2007GC001788.
Raitzsch, M., Völker, D., and Heubeck, C. (2007): Neogene sedimentary and mass-wasting processes on the continental margin off south-central Chile inferred from dredge samples. – Marine Geology, 244(1-4): 166-183, https://doi.org/10.1016/j.margeo.2007.06.007.
Theses
Raitzsch, M. (2008): Mg/Ca, B/Ca and U/Ca in shells of benthic foraminifers: Tracers of temperature and seawater carbonate chemistry – PhD thesis. Universität Bremen, FB Geowissenschaften, available from http://nbn-resolving.de/urn:nbn:de:gbv:46-diss000112160.
Raitzsch, M. (2004): Untersuchung subrezenter Sedimentations- und Akkretionsprozesse am aktiven Kontinentalrand des südlichen Zentralchile (36°–40°S) anhand von Dredgeproben (SO161-5). – Diploma thesis. Freie Universität Berlin, FR Geologie, 114 pp.
Raitzsch, M. (2004): Diplomkartierung in den Corbières Orientales (Dept. Aude/Frankreich). – Diploma mapping. Freie Universität Berlin, FR Geologie, 63 pp.
Privacy policy & Legal notice
Read more
Privacy policy
Personal data (usually referred to just as “data” below) will only be processed by me to the extent necessary and for the purpose of providing a functional and user-friendly website, including its contents, and the services offered there.
Per Art. 4 No. 1 of Regulation (EU) 2016/679, i.e. the General Data Protection Regulation (hereinafter referred to as the “GDPR”), “processing” refers to any operation or set of operations such as collection, recording, organization, structuring, storage, adaptation, alteration, retrieval, consultation, use, disclosure by transmission, dissemination, or otherwise making available, alignment, or combination, restriction, erasure, or destruction performed on personal data, whether by automated means or not.
The following privacy policy is intended to inform you in particular about the type, scope, purpose, duration, and legal basis for the processing of such data either under my own control or in conjunction with others. I also inform you below about the third-party components I use to optimize my website and improve the user experience which may result in said third parties also processing data they collect and control.
My privacy policy is structured as follows:
I. Information about me as controller of your data
II. The rights of users and data subjects
III. Information about the data processing
I. Information about me as controller of your data
The party responsible for this website (the “controller”) for purposes of data protection law is:
Markus Raitzsch
Papenbütteler Weg 8
27729 Hambergen
Telephone: 0176 32349183
E-Mail: 
II. The rights of users and data subjects
With regard to the data processing to be described in more detail below, users and data subjects have the right
- to confirmation of whether data concerning them is being processed, information about the data being processed, further information about the nature of the data processing, and copies of the data (cf. also Art. 15 GDPR);
- to correct or complete incorrect or incomplete data (cf. also Art. 16 GDPR);
- to the immediate deletion of data concerning them (cf. also Art. 17 DSGVO), or, alternatively, if further processing is necessary as stipulated in Art. 17 Para. 3 GDPR, to restrict said processing per Art. 18 GDPR;
- to receive copies of the data concerning them and/or provided by them and to have the same transmitted to other providers/controllers (cf. also Art. 20 GDPR);
- to file complaints with the supervisory authority if they believe that data concerning them is being processed by the controller in breach of data protection provisions (see also Art. 77 GDPR).
In addition, the controller is obliged to inform all recipients to whom it discloses data of any such corrections, deletions, or restrictions placed on processing the same per Art. 16, 17 Para. 1, 18 GDPR. However, this obligation does not apply if such notification is impossible or involves a disproportionate effort. Nevertheless, users have a right to information about these recipients.
Likewise, under Art. 21 GDPR, users and data subjects have the right to object to the controller’s future processing of their data pursuant to Art. 6 Para. 1 lit. f) GDPR. In particular, an objection to data processing for the purpose of direct advertising is permissible.
III. Information about the data processing
Your data processed when using my website will be deleted or blocked as soon as the purpose for its storage ceases to apply, provided the deletion of the same is not in breach of any statutory storage obligations or unless otherwise stipulated below.
Server data
For technical reasons, the following data sent by your internet browser to me or to my server provider will be collected, especially to ensure a secure and stable website: These server log files record the type and version of your browser, operating system, the website from which you came (referrer URL), the webpages on my site visited, the date and time of your visit, as well as the IP address from which you visited my site.
The data thus collected will be temporarily stored, but not in association with any other of your data.
The basis for this storage is Art. 6 Para. 1 lit. f) GDPR. My legitimate interest lies in the improvement, stability, functionality, and security of my website.
The data will be deleted within no more than seven days, unless continued storage is required for evidentiary purposes. In which case, all or part of the data will be excluded from deletion until the investigation of the relevant incident is finally resolved.
I maintain an online presence on Twitter to present my company and my services and to communicate with customers/prospects. Twitter is a service provided by Twitter Inc., 1355 Market Street, Suite 900, San Francisco, CA 94103, USA.
I would like to point out that this might cause user data to be processed outside the European Union, particularly in the United States. This may increase risks for users that, for example, may make subsequent access to the user data more difficult. I also do not have access to this user data. Access is only available to Twitter.
The privacy policy of Twitter can be found at:
https://twitter.com/privacy
I maintain an online presence on LinkedIn to present my company and my services and to communicate with customers/prospects. LinkedIn is a service of LinkedIn Ireland Unlimited Company, Wilton Plaza, Wilton Place, Dublin 2, Irland, a subsidiary of LinkedIn Corporation, 1000 W. Maude Avenue, Sunnyvale, CA 94085, USA.
I would like to point out that this might cause user data to be processed outside the European Union, particularly in the United States. This may increase risks for users that, for example, may make subsequent access to the user data more difficult. I also do not have access to this user data. Access is only available to LinkedIn.
The LinkedIn privacy policy can be found here:
https://www.linkedin.com/legal/privacy-policy
Legacy notice
Information in accordance with Section 5 TMG
Markus Raitzsch
Papenbütteler Weg 8
27729 Hambergen
Contact information
Telephone: 0176 32349183
E-Mail:
Internet address: https://markus-r.com
Disclaimer
Accountability for content
The contents of my pages have been created with the utmost care. However, I cannot guarantee the contents’ accuracy, completeness or topicality. According to statutory provisions, I am furthermore responsible for my own content on these web pages. In this matter, please note that I am not obliged to monitor the transmitted or saved information of third parties, or investigate circumstances pointing to illegal activity. My obligations to remove or block the use of information under generally applicable laws remain unaffected by this as per §§ 8 to 10 of the Telemedia Act (TMG).
Accountability for links
Responsibility for the content of external links (to web pages of third parties) lies solely with the operators of the linked pages. No violations were evident to me at the time of linking. Should any legal infringement become known to me, I will remove the respective link immediately.
Copyright
Our web pages and their contents are subject to German copyright law. Unless expressly permitted by law, every form of utilizing, reproducing or processing works subject to copyright protection on our web pages requires the prior consent of the respective owner of the rights. Individual reproductions of a work are only allowed for private use. The materials from these pages are copyrighted and any unauthorized use may violate copyright laws.
Social media links via graphics
I also integrate the following social media sites into my website. The integration takes place via a linked graphic of the respective site. The use of these graphics stored on my own servers prevents the automatic connection to the servers of these networks for their display. Only by clicking on the corresponding graphic will you be forwarded to the service of the respective social network.
Once you click, that network may record information about you and your visit to my site. It cannot be ruled out that such data will be processed in the United States.
Initially, this data includes such things as your IP address, the date and time of your visit, and the page visited. If you are logged into your user account on that network, however, the network operator might assign the information collected about your visit to my site to your personal account. If you interact by clicking Like, Share, etc., this information can be stored your personal user account and possibly posted on the respective network. To prevent this, you need to log out of your social media account before clicking on the graphic. The various social media networks also offer settings that you can configure accordingly.
The following social networks are integrated into my site by linked graphics:
Twitter
Twitter Inc., 795 Folsom St., Suite 600, San Francisco, CA 94107, USA.
Privacy Policy: https://twitter.com/privacy
LinkedIn
LinkedIn Ireland Unlimited Company, Wilton Plaza, Wilton Place, Dublin 2, Irland, a subsidiary of LinkedIn Corporation, 1000 W. Maude Avenue, Sunnyvale, CA 94085 USA.
Privacy Policy: https://www.linkedin.com/legal/privacy-policy