# 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.