#' @title Calculate carbonate system from alkalinity and DIC
#'
#' @description This function calculates the full carbonate system including
#' pCO2, [CO2], [HCO3], [CO3], and pH based on supplied total alkalinity (TAlk)
#' and DIC (TCO2).
#'
#' @usage
#'
#' @param TAlk Total alkalinity in units of umol per liter
#' @param TCO2 Total CO2 (DIC) in units of umol per liter
#' @param temperature Water temperature in degrees Celsius
#' @param salinity Salinity in g per kg
#' @param pressure Surface pressure in atmospheres
#'@return A list with elements for pCO2, [CO2], [HCO3], [CO3], [TCO2], TAlk, and pH
#'@note Units of concentration are micromoles per liter (umol/L)
#'@references
#'@author
#'@examples
#'@export
get.CO2Sys.TAlk_TCO2 <- function(TAlk, TCO2, temperature, salinity = 0, pressure = 1){
# Error handling
if(is.null(temperature)) stop("'temperature' missing with no defualt.")
if(is.null(TCO2)) stop("'TCO2' missing with no defualt.")
if(is.null(TAlk)) stop("'TAlk' missing with no defualt.")
if(!is.numeric(TAlk)) stop("'TAlk' must be numeric.")
if(!is.numeric(TCO2)) stop("'TCO2' must be numeric.")
vars <- c("pCO2", "CO2_", "HCO3", "CO3", "TCO2", "TAlk", "pH")
CO2Sys <- vector("list", length(vars))
names(CO2Sys) <- vars
H <- vector("numeric", length(TAlk))
CO2Sys[["TCO2"]] <- TCO2 # Units umol/L
TCO2_molL <- TCO2 * 10^-6 # Units mol/L
CO2Sys[["TAlk"]] <- TAlk # Units umol/L
TAlk_molL <- TAlk * 10^-6 # Units mol/L
K0 <- get.K0(temperature = temperature, salinity = salinity, pressure = pressure)
K1_ <- get.K1(temperature = temperature, salinity = salinity)
K2 <- get.K2(temperature = temperature, salinity = salinity)
f <- function(x, TAlk, DIC, K1_, K2, T_, S) {
# TAlk = 0.001
# print(TAlk); print(CO2_); print(K1_); print(K2); print(T_); print(S)
HCO3 <- DIC * x * K1_ / (x * x + K1_ * x + K1_ * K2)
CO3 <- DIC * K1_ * K2 / (x * x + K1_ * x + K1_ * K2)
OH <- get.OH(temperature = T_, salinity = S, H = x)
OUT <- HCO3 + 2 * CO3 + OH - TAlk
return(OUT)
}
for(i in 1:length(TAlk_molL)){
# print(f(10^(-9.5),
# TAlk = TAlk_molL[i],
# DIC = TCO2_molL[i],
# K1_ = K1_[i], K2 = K2[i],
# T_ = temperature[i],
# S = salinity[i]))
# print(f(10^(-3.5),
# TAlk = TAlk_molL[i],
# DIC = TCO2_molL[i],
# K1_ = K1_[i], K2 = K2[i],
# T_ = temperature[i],
# S = salinity[i]))
H[i] <- tryCatch(uniroot(f, c(1e-10, 10^(-3.5)), trace=1,
# H[i] <- tryCatch(uniroot(f, c(10^(-9.5), 10^(-3.5)),
# TAlk = 0.001,
TAlk = TAlk_molL[i],
DIC = TCO2_molL[i],
K1_ = K1_[i], K2 = K2[i],
T_ = temperature[i],
S = salinity[i],
tol = 1e-30)$root,
error=function(e){
# print(TAlk); print(CO2_); print(K1_); print(K2); print(T_); print(S)
# message('fail T')
return(NA)
})
}
temp <- H * H + K1_ * H + K1_ * K2
HCO3_molL <- TCO2_molL * K1_ * H / temp
CO3_molL <- TCO2_molL * K1_ * K2 / temp
CO2_molL <- H * HCO3_molL / K1_
pCO2_atm <- CO2_molL / K0
CO2Sys[["pH"]] <- -log10(H)
CO2Sys[["CO2_"]] <- CO2_molL * 10^6 #Convert to umol/L
CO2Sys[["HCO3"]] <- HCO3_molL * 10^6
CO2Sys[["CO3"]] <- CO3_molL * 10^6
CO2Sys[["pCO2"]] <- pCO2_atm * 10^6
return(CO2Sys)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.