R/get.CO2Sys_TAlk_TCO2.R

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

#' @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)
}