development/get.CO2Sys_pCO2_TCO2.R

In gholtgrieve/gassyPants: Functions Useful for Working With Gasses Dissolved in Water

#' @title Calculate carbonate system from pCO2 and TCO2 (DIC)
#'
#' @description This function calculates the full carbonate system including
#' pCO2, [CO2], [HCO3], [CO3], and pH based on supplied pCO2 and DIC (TCO2).
#'
#' @usage
#'
#' @param pCO2 Partial pressure of CO2 in microatmospheres
#' @param TCO2 Total CO2 (DIC) in units of umol/L
#' @param temperature Water temperature in degrees Celsius
#' @param salinity Salinity in g/kg
#' @param pressure Surface pressure in atmospheres
#' @return  A list with elements for pCO2, [CO2], [HCO3], [CO3], [TCO2], TAlk, and pH.
#'
#' @note Concentrations are in units of umol/L.
#' @references
#' @author Gordon W. Holtgrieve
#' @examples

#'@export

get.CO2Sys.pCO2_TCO2 <- function(pCO2, TCO2, temperature, salinity=0, pressure=1){

  vars <- c("pCO2", "CO2", "HCO3", "CO3", "TCO2", "TAlk", "pH")
  CO2Sys <- vector("list", length(vars))
  names(CO2Sys) <- vars

  CO2Sys[["pCO2"]] <- pCO2  # uatm
  pCO2_atm <- pCO2 / 10^6   # atm
  CO2Sys[["TCO2"]] <- TCO2  # umol/L
  TCO2_molL <- TCO2 / 10^6  # mol/L

  K0 <- get.K0(temperature = temperature, salinity = salinity, pressure = pressure)  # mol/L/atm
  K1 <- get.K1(temperature = temperature, salinity = salinity)  # mol/L
  K2 <- get.K2(temperature = temperature, salinity = salinity)  # mol/L

  CO2_molL <- pCO2_atm * K0
  K <- K1 / K2
  b <- K * K0 * pCO2_atm
  c <- (K * K0 * pCO2_atm) * (K0 * pCO2_atm - TCO2_molL)
  D <- b * b - 4 * c
  HCO3_molL <- (1/2) * (-b + sqrt(D))
  CO3_molL <- TCO2_molL - CO2_molL - HCO3_molL
  H <- K1 * CO2_molL / HCO3_molL
  OH <- get.OH(temperature, salinity, H)
  TAlk <- HCO3_molL + 2 * CO3_molL + OH - H
  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[["TAlk"]] <- TAlk * 10^6

  return(CO2Sys)
}