R/KSSL_VG_model.R

In soilDB: Soil Database Interface

# define van Genuchten model as a function
# this is tailored to the parameters stored in our KSSL data
# https://en.wikipedia.org/wiki/Water_retention_curve
.vg <- function(phi, theta_r, theta_s, alpha, n) {
  theta_r + ((theta_s - theta_r) / ((1 + (alpha * phi) ^ n) ^ (1 - 1 / n)))
}




## notes: Rosetta units for alpha and npar are log10(1/cm) and log10([-])
# VG_params: table of VG parameters from KSSL / Rosetta:  alpha and npar are in log10 form
# phi_min: lower limit for water potential in kPa
# phi_max: upper limit for water potential in kPa
# pts: number of points to include in curve


#' @title Develop a Water Retention Curve from KSSL Data
#' 
#' @description Water retention curve modeling via van Genuchten model and KSSL data.
#'
#' @param VG_params \code{data.frame} or \code{list} object with the parameters of the van Genuchten model, see details
#' 
#' @param phi_min lower limit for water potential in kPa
#' 
#' @param phi_max upper limit for water potential in kPa
#' 
#' @param pts number of points to include in estimated water retention curve
#' 
#' @details 
#' 
#' This function was developed to work with measured or estimated parameters of the \href{https://en.wikipedia.org/wiki/Water_retention_curve}{van Genuchten model}, as generated by the \href{https://www.ars.usda.gov/pacific-west-area/riverside-ca/agricultural-water-efficiency-and-salinity-research-unit/docs/model/rosetta-model/}{Rosetta model}. As such, \code{VG_params} should have the following format and conventions:
#' \describe{
#'   \item{theta_r}{saturated water content, values should be in the range of \{0, 1\}}
#'   \item{theta_s}{residual water content, values should be in the range of \{0, 1\}}
#'   \item{alpha}{related to the inverse of the air entry suction, function expects log10-transformed values with units of 1/cm}
#'   \item{npar}{index of pore size distribution, function expects log10-transformed values (dimensionless)}
#' }
#'
#' @return 
#' A list with the following components:
#' \describe{
#'   \item{VG_curve}{estimated water retention curve: paired estimates of water potential (phi) and water content (theta)}
#'   \item{VG_function}{spline function for converting water potential (phi, units of kPa) to estimated volumetric water content (theta, units of percent, range: \{0, 1\})}
#'   \item{VG_inverse_function}{spline function for converting volumetric water content (theta, units of percent, range: \{0, 1\}) to estimated water potential (phi, units of kPa)}
#' }
#' 
#' @author D.E. Beaudette
#' 
#' @note A practical example is given in the \href{http://ncss-tech.github.io/AQP/soilDB/fetchSCAN-demo.html}{fetchSCAN tutorial}.
#' 
#' @references 
#' 
#' \href{https://en.wikipedia.org/wiki/Water_retention_curve}{water retention curve estimation}
#' 
# 'van Genuchten, M.Th. (1980). "A closed-form equation for predicting the hydraulic conductivity of unsaturated soils". Soil Science Society of America Journal. 44 (5): 892-898. 
#' 
#' 
#' 
#' @export
#'
#' @examples
#' 
#' # basic example
#' d <- data.frame(
#'   theta_r = 0.0337216, 
#'   theta_s = 0.4864061, 
#'   alpha = -1.581517, 
#'   npar = 0.1227247
#' )
#' 
#' vg <- KSSL_VG_model(d)
#' 
#' str(vg)
#' 
KSSL_VG_model <- function(VG_params, phi_min = 10^-6, phi_max = 10^8, pts = 100) {
  
  # sanity check, expected columns
  if( any(! c('theta_r', 'theta_s', 'alpha', 'npar') %in% names(VG_params)) ) {
    message('one or more required column is missing')
    return(list(VG_curve = NULL, VG_inverse_function = NULL))
  }
  
  # subset to to expected columns
  VG_params <- VG_params[, c('theta_r', 'theta_s', 'alpha', 'npar')]
  
  # sanity check: no NA allowed
  # return NULL if present
  if(any(is.na(VG_params))) {
    message('one or more required value is NA')
    return(list(VG_curve = NULL, VG_inverse_function = NULL))
  }
    
  
  # useful range in kPa suctions
  phi <- 10^seq(log(phi_min, base=10), log(phi_max, base=10), length.out = pts)
  m <- data.frame(phi=phi)
  
  # Rosetta units for alpha and npar are 1/cm and [-]
  # convert kPa to cm of H20
  h <- m$phi * 10.19716
  
  # compute theta given measured parameters and sequence of phi
  m$theta <- .vg(
    phi = h,
    theta_r = VG_params$theta_r,
    theta_s = VG_params$theta_s,
    alpha = 10 ^ (VG_params$alpha),
    n = 10 ^ (VG_params$npar)
  )
  
  # use splines to fit standard model: phi -> theta
  # scale of theta {0,1}
  # phi: units of kPa
  vg.fwd <- splinefun(m$phi, m$theta)
  
  # use splines to fit inverse model: theta -> phi
  # scale of theta {0,1}
  # phi: units of kPa
  vg.inv <- splinefun(m$theta, m$phi)
  
  # return curve and spline function
  return(list(VG_curve=m, VG_function=vg.fwd, VG_inverse_function=vg.inv))
}