R/get.CO2Sys_TAlk_pCO2.R

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

#' @title Calculate carbonate system from alkalinity and pCO2
#'
#' @description This function calculates the full carbonate system including
#' pCO2, [CO2], [HCO3], [CO3], and pH based on supplied total alkalinity (TAlk)
#' and partial pressure of dissolved CO2 (pCO2).
#'
#' @usage
#'
#' @param TAlk Total alkalinity in units of umol per liter
#' @param pCO2 partial pressure of CO2 in units of uatm
#' @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 Gordon W. Holtgrieve
#' @examples

#' @export

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

  # Error handling
  if(is.null(temperature)) stop("'temperature' missing with no defualt.")
  if(is.null(pCO2)) stop("'pCO2' 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(pCO2)) stop("'pCO2' 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[["pCO2"]] <- pCO2  # uatm
  pCO2_atm <- pCO2 / 10^6   # atm
  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)

  CO2Sys[["CO2_"]] <- pCO2 * K0         # Units umol/L
  CO2_molL <- CO2Sys[["CO2_"]] / 10^6  # Units mol/L

  f <- function(x, TAlk, CO2_, K1_, K2,  T_, S) {
      # TAlk = 0.001
      # print(TAlk); print(CO2_); print(K1_); print(K2); print(T_); print(S)
    DIC <- CO2_ * (1 + K1_/x + K1_ * K2/(x * x))
    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],
      #       CO2_ = CO2_molL[i],
      #       K1_ = K1_[i], K2 = K2[i],
      #       T_ = temperature[i],
      #       S = salinity[i]))
      # print(f(10^(-3.5),
      #       TAlk  = TAlk_molL[i],
      #       CO2_ = CO2_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] <- uniroot(f, c(10^(-9.5), 10^(-3.5)),
                    # TAlk  = 0.001,
                    TAlk  = TAlk_molL[i],
                    CO2_ = CO2_molL[i],
                    K1_ = K1_[i], K2 = K2[i],
                    T_ = temperature[i],
                    S = salinity[i],
                    tol = 1e-30)$root,
        error=function(e){
            # message('fail p')
            return(NA)
        })
  }

  HCO3_molL <- K1_ * CO2_molL / H
  CO3_molL <- K2 * HCO3_molL / H
  DIC_molL <- CO2_molL + HCO3_molL + CO3_molL

  CO2Sys[["pH"]] <- -log10(H)
  CO2Sys[["HCO3"]] <- HCO3_molL * 10^6  #Convert to umol/L
  CO2Sys[["CO3"]] <- CO3_molL * 10^6
  CO2Sys[["TCO2"]] <- DIC_molL * 10^6

  return(CO2Sys)
}