#' Error function
#'
#' @param x x
#' @return Error function result
#' @references \url{https://stat.ethz.ch/R-manual/R-devel/library/stats/html/Normal.html}
#' @importFrom stats pnorm
erf <- function(x) {
2 * stats::pnorm(x * sqrt(2)) - 1
}
#' Helper function hantushS
#'
#' @param x distance between 0 and half length of recharge basin
#' @param alpha alpha
#' @param beta beta
#' @return Hantush star
#' @references p.22, \url{https://pubs.usgs.gov/sir/2010/5102/support/sir2010-5102.pdf}
hantushS <- function(x, alpha, beta) {
xRoot <- sqrt(x)
erf(alpha / xRoot) * erf(beta / xRoot)
}
#' Hantush function Sstar
#' @param alpha alpha
#' @param beta beta
#' @param dbg If True additional messages on integration quality of function
#' hantushSstar are printed on screen
#' @return Hantush Sstar result
#' @references p.22, \url{https://pubs.usgs.gov/sir/2010/5102/support/sir2010-5102.pdf}
#' @importFrom stats integrate pnorm
hantushSstar <- function(alpha, beta, dbg) {
x <- stats::integrate(
f = hantushS, lower = 0, upper = 1,
alpha = alpha, beta = beta # arguments passed to hansuhS
)
if (dbg) {
msg <- sprintf("Sstar: %3.5f with an absolute error < %f", x$value, x$abs.error)
message(msg)
}
return(x$value)
}
#' Hantush equation base properties
#'
#' @param time time elapsed since recharge began (T), (Default: 1.5 days)
#' @param basinWidth half width of the recharge basin (L), (Default: 10 m)
#' @param basinLength half length of the recharge basin (L), (Default: 10 m)
#' @param infiltrationRate recharge (infiltration) rate (L/T), (Default: 0.5 m/d)
#' @param horizConductivity horizontal hydraulic conductivity (L/T), (Default: 10 m/d)
#' @param iniHead initial head (height of the water table above the base of the
#' aquifer);(L), (Default: 10)
#' @param specificYield specific yield (Default: 0.2)
#' @param numberTimeSteps number of time steps to be used for average aquifer
#' thickness calculation (Default: 150)
#' @return Base properties for Hantush equation
#' @export
#' @references p.22, \url{https://pubs.usgs.gov/sir/2010/5102/support/sir2010-5102.pdf}
baseProperties <- function(time = 10,
basinWidth = 10, ### m
basinLength = 10, ### m
infiltrationRate = 0.5, ### 0.5 m/day
horizConductivity = 10, ### m/d
iniHead = 10, ### meter
specificYield = 0.2,
numberTimeSteps = 150) {
x <- list(
time = time,
basinWidth = basinWidth,
basinLength = basinLength,
infiltrationRate = infiltrationRate,
horizConductivity = horizConductivity,
iniHead = iniHead,
specificYield = specificYield,
numberTimeSteps = numberTimeSteps
)
return(x)
}
#' Hantush equation
#'
#' @param x distance from the center of the recharge basin in the x direction (L)
#' @param y distance from the center of the recharge basin in the y direction (L)
#' @param baseProps basic model properties as retrieved by baseProperties()
#' @param dbg If True additional messages on integration quality of function
#' hantushSstar are printed on screen
#' @return Head at a given time after recharge begins
#' @export
#' @references p.22, \url{https://pubs.usgs.gov/sir/2010/5102/support/sir2010-5102.pdf}
#' @seealso \code{\link{baseProperties}} for basic model properties
hantush <- function(x = 0,
y = 0,
baseProps = baseProperties(),
dbg = TRUE) {
indices <- seq_len(baseProps$numberTimeSteps)
timeSteps <- baseProps$time / baseProps$numberTimeSteps * indices
h <- vector(length = baseProps$numberTimeSteps)
L <- baseProps$basinLength
sy <- baseProps$specificYield
W <- baseProps$basinWidth
for (i in indices)
{
if (i == 1) {
avgHead <- baseProps$iniHead
} else {
avgHead <- 0.5 * (baseProps$iniHead + h[i - 1])
}
if (dbg) {
cat(sprintf(
"Elapsed time: %4.5f ; Havg: %4.5f Hi-1: %4.5f\n\n",
timeSteps[i], avgHead, h[i - 1]
))
}
diffusivity <- baseProps$horizConductivity * avgHead / sy
factor1 <- (baseProps$infiltrationRate * avgHead * timeSteps[i]) / (2 * sy)
divisor <- sqrt(4 * timeSteps[i] * baseProps$horizConductivity * avgHead / baseProps$specificYield)
alphaPlus <- (L + x) / divisor
alphaMinus <- (L - x) / divisor
betaPlus <- (W + y) / divisor
betaMinus <- (W - y) / divisor
factor2 <-
hantushSstar(alphaPlus, betaPlus, dbg) +
hantushSstar(alphaPlus, betaMinus, dbg) +
hantushSstar(alphaMinus, betaPlus, dbg) +
hantushSstar(alphaMinus, betaMinus, dbg)
h[i] <- sqrt(factor1 * factor2 + baseProps$iniHead^2)
}
res <- data.frame(
timeSteps = timeSteps,
x = rep(x, baseProps$numberTimeSteps),
y = rep(y, baseProps$numberTimeSteps),
head = h
)
return(list(
res = res,
baseProps = baseProps
))
}
#' Hantush distance: for multiple coordinates
#'
#' @param x vector with distances from the center of the recharge basin in the x
#' direction (L), (Default: each meter between 0-100m)
#' @param y vector with distances from the center of the recharge basin in the y
#' direction (L), (Default: 0 times length of x)
#' @param baseProps basic model properties as retrieved by baseProperties()
#' @param dbg If True additional debug messages are printed on screen
#' @return Head at a given time after recharge begins
#' @export
#' @seealso \code{\link{hantush}} for parameterizing the Hantush equation
hantushDistances <- function(x = 0:10,
y = rep(0, length(x)),
baseProps = baseProperties(),
dbg = TRUE) {
out <- NULL
for (i in 1:length(x))
{
tmp <- hantush(x = x[i], y = y[i], baseProps = baseProps, dbg = dbg)$res
out <- rbind(out, tmp)
}
if (dbg) {
cat(sprintf("Head of %4.5f m at position %5.2f m (x), %5.2f (y)\n", out$head[i], x[i], y[i]))
}
maxSimTime <- max(out$timeSteps)
simTime <- out[out$timeSteps == maxSimTime, c("x", "y", "head") ]
out <- list(
timeSteps = out,
simTime = simTime,
baseProps = baseProps
)
return(out)
}
#' Hantush distances & base properties: allows input of vector of x,y coordinates and
#' also a vector for one of the base properties
#'
#' @param x vector with distances from the center of the recharge basin in the x
#' direction (L), (Default: every 5 meter between 0-200m)
#' @param y vector with distances from the center of the recharge basin in the y
#' direction (L), (Default: 0 times length of x)
#' @param baseProps as retrieved by baseProperties(), but one property is allowed
#' to be a vector (e.g. infiltrationRate = c(1,2,4))
#' @param dbg If True additional debug messages are printed on screen
#' @return List with sublists "dat" (x,y,head & WLincrease), "changedBaseProp.Name"
#' (name of base property with multiple values) and "baseProps" (complete base
#' parameterisation)
#' @seealso \code{\link{baseProperties}} for basic model properties
#' @export
#' @examples
#' baseProps <- baseProperties(
#' time = 2^(0:6),
#' infiltrationRate = 1,
#' basinWidth = 10,
#' basinLength = 50,
#' horizConductivity = 10,
#' iniHead = 10,
#' specificYield = 0.2,
#' numberTimeSteps = 10
#' )
#' res <- hantushDistancesBaseProps(baseProps = baseProps)
#' cols <- length(unique(res$dat[[res$changedBaseProp.Name]]))
#' mainTxt <- sprintf("Changed baseProperty: %s", res$changedBaseProp.Name)
#' xyplot(WLincrease ~ x,
#' groups = res$dat[[res$changedBaseProp.Name]],
#' data = res$dat,
#' type = "b",
#' auto.key = list(columns = cols),
#' main = mainTxt
#' )
hantushDistancesBaseProps <- function(x = seq(0, 200, 5),
y = rep(0, length(x)),
baseProps = baseProperties(
time = 2^(0:6),
infiltrationRate = 1,
basinWidth = 10,
basinLength = 50,
horizConductivity = 10,
iniHead = 10,
specificYield = 0.2
),
dbg = FALSE) {
newRes <- NULL
multiplePars <- unlist(lapply(X = baseProps, FUN = length))
multipleParVals <- NULL
changedBaseProp.Name <- ""
if (any(multiplePars > 1)) {
selMultiplePars <- which(multiplePars > 1)
if (length(selMultiplePars) == 1) {
multipleParVals <- baseProps[selMultiplePars]
changedBaseProp.Name <- names(multipleParVals)
} else {
variedPars <- paste(names(selMultiplePars), collapse = ", ")
msg <- sprintf("Only one base property can be a vector!
Now you defined multiple values the following base properties:\n%s", variedPars)
stop(msg)
}
} else {
multipleParVals <- baseProps$time ### select a random variable from baseProp
}
baseProps.sel <- baseProps
for (parVal in as.vector(multipleParVals[[1]]))
{
if (any(multiplePars > 1)) {
cat(sprintf(
"Calculating results for changed baseProperty '%s': %3.5f ... ",
changedBaseProp.Name,
parVal
))
}
baseProps.sel[[ changedBaseProp.Name]] <- parVal
res <- hantushDistances(
x = x,
y = y,
baseProps = baseProps.sel,
dbg = dbg
)
res$simTime$WLincrease <- res$simTime$head - res$baseProps$iniHead
if (any(multiplePars > 1)) {
changedBaseProp <- data.frame(rep(parVal, nrow(res$simTime)))
names(changedBaseProp) <- names(multipleParVals)
tmpRes <- cbind(changedBaseProp, res$simTime)
} else {
tmpRes <- res$simTime
}
newRes <- rbind(newRes, tmpRes)
if (any(multiplePars > 1)) {
cat("Done!\n")
}
}
return(list(
dat = newRes,
changedBaseProp.Name = changedBaseProp.Name,
baseProps = baseProps
))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.