Nothing
#' Spectral response for a sealed well
#'
#' This is the primary function to calculate the response
#' for a sealed well.
#'
#' @details
#' The response depends strongly on the physical properties
#' given. Default values are assumed where reasonable--for instance,
#' the pore-fluid is assumed to be water--but considerable care
#' should be invested in the choice of
#' parameters, unless the function is used in an optimization scheme.
#'
#' Assumed values are:
#' \tabular{rlrl}{
#' \code{Avs} \tab 1 \tab \tab amplification factor for volumetric strain\cr
#' \code{Aw} \tab 1 \tab \tab amplification factor for water well\cr
#' }
#'
#' The responses returned here are,
#' effectively, the amplification of water levels in a well, relative to
#' the aquifer strain; or
#' \deqn{Z = \frac{z}{\epsilon} \equiv \frac{p}{\rho g \epsilon}}
#' If \code{as.pressure=TRUE}, then the responses are scaled by
#' \code{rho*grav} so that they represent water pressure relative to
#' aquifer strain:
#' \deqn{Z = \frac{p}{\epsilon}}
#'
#' Not all parameters need to be given, but should be.
#' For example, if
#' either \code{rho} or \code{grav} are not specified, they
#' are taken from \code{\link{constants}}.
#' \emph{Parameters which do not end in \code{.} do
#' not need to be specified (they may be excluded); if
#' they are missing, warnings will be thrown.}
#'
#' @name well_response
#' @export
#'
#' @param omega frequency, (see \code{freq.units})
#' @param T. effective aquifer transmissivity \eqn{[m^2/s]}
#' @param S. well storativity, \eqn{[unitless]}
#' @param Vw. well volume, \eqn{[m^3]}
#' @param Rs. radius of screened portion, \eqn{[m]}
#' @param Ku. undrained bulk modulus, \eqn{[Pa]}
#' @param B. Skempton's coefficient, \eqn{[unitless, bounded]}
#' @param Avs amplification factor for volumetric strain \eqn{E_{kk,obs}/E_{kk}}, \eqn{[]}
#' @param Aw amplification factor of well volume change for \eqn{E_{kk}}, \eqn{[]}
#' @param rho fluid density \eqn{[kg/m^3]}
#' @param Kf bulk modulus of fluid, \eqn{[Pa]}
#' @param grav local gravitational acceleration \eqn{[m/s^2]}
#' @param freq.units set the units of \code{omega}
#' @param as.pressure logical; should the response for water pressure? (default is water height)
#' @param ... additional arguments
#'
#' @return An object with class \code{'wrsp'}
#'
#' @author A. J. Barbour
#'
#' @references See \code{\link{kitagawa-package}} for references and more background.
#' @seealso \code{\link{open_well_response}} for the open-well equivalents
#' \code{\link{wrsp-methods}} for a description of the class \code{'wrsp'} and its methods,
#' \code{\link{sensing_volume}} to easily estimate the volume \code{Vw.}, and
#' \code{\link{kitagawa-package}} for references and more background.
#' @family WellResponseFunctions
#'
#' @examples
#' #### dummy example
#' well_response(1:10, T.=1, S.=1, Vw.=1, Rs.=1, Ku.=1, B.=1)
#'
#' #### a more physically realistic calculation:
#' # Physical params applicable for B084 borehole
#' # (see: http://pbo.unavco.org/station/overview/B084/ for details)
#' #
#' Rc <- 0.0508 # m, radius of water-sensing (2in)
#' Lc <- 146.9 # m, length of grouted region (482ft)
#' Rs <- 3*Rc # m, radius of screened region (6in)
#' Ls <- 9.14 # m, length of screened region (30ft)
#' #
#' # calculate the sensing volume for the given well parameters
#' Volw <- sensing_volume(Rc, Lc, Rs, Ls) # m**3, ~= 1.8
#' #
#' Frqs <- 10**seq.int(from=-4,to=0,by=0.1) # log10-space
#' head(Rsp <- well_response(omega=Frqs, T.=1e-6, S.=1e-5,
#' Vw.=Volw, Rs.=Rs, Ku.=40e9, B.=0.2, freq.units="Hz"))
#'
#' # Access plot.wrsp:
#' plot(Rsp)
#'
well_response <- function(omega, T., S., Vw., Rs., Ku., B., ...) UseMethod("well_response")
#' @rdname well_response
#' @export
well_response.default <- function(omega, T., S., Vw., Rs., Ku., B.,
Avs, Aw, rho, Kf, grav,
freq.units=c("rad_per_sec","Hz"),
as.pressure=TRUE, ...){
# Enforce units of omega to be radians/sec
freq.units <- match.arg(freq.units)
fc <- switch(freq.units, rad_per_sec=1, Hz=2*pi)
omega <- fc*omega
#
# Setup constants
defA <- 1
if (missing(Avs)){
Avs <- defA
}
if (missing(Aw)){
Aw <- defA
}
const <- kitagawa::constants(FALSE)
if (missing(rho)){
rho <- const$water$density
}
if (missing(Kf)){
Kf <- const$water$bulkmod
}
if (missing(grav)){
grav <- const$gravity
}
rhog <- rho*grav
#
# Alpha function
Alpha. <- omega_constants(omega, c.type="alpha", S.=S., T.=T., Rs.=Rs.)
# A1, and A2 (functions of Phi and Psi, calculated internally)
Amat <- alpha_constants(Alpha., c.type="A")
stopifnot(ncol(Amat)==7)
# A1,2 are in Mod(A.[,6:7])
A12 <- matrix(Mod(Amat[,6:7]),ncol=2) # is complex, but imag is zero, so == abs
stopifnot(ncol(A12)==2)
A1 <- A12[,1]
A2 <- A12[,2]
#
TVFRG <- 2 * pi * T. / omega / Vw. / rhog
#
#check Skemptons coefficient: it should be bound to [0,1]
.in0to1(B.)
#
XX. <- Ku. * B. / Aw * TVFRG - A2
rNum. <- XX. * XX. + A1 * A1
#
YY. <- Kf * TVFRG - A2
rDen. <- YY. * YY. + A1 * A1
##
## complex response EQ 17
cNum <- complex(real=(Ku. * B. / Aw * TVFRG - A2), imaginary=A1)
cDen <- complex(real=(Kf * TVFRG - A2), imaginary=A1)
#
# z / epsilon == p / epsilon * rho * g
# but is there a mistake?
# / rhog
wellresp <- -1 * Kf * Aw / Avs * cNum / cDen
##
# optionally scale to p/epsilon
rhog <- ifelse(as.pressure, rhog, 1)
Zunits <- ifelse(as.pressure, "P/E", "Z/E")
wellresp <- wellresp * rhog
##
omega <- omega/fc
toret <- list(Aquifer=list(Transmiss=T., Storativ=S., Diffusiv=T./S.),
Well=list(Volume=Vw., ScreenRad=Rs.),
Solid=list(BulkModU=Ku.,Skemp=B.),
Fluid=list(BulkMod=Kf, Density=rho),
Amplification=list(Ekk=Avs, Well=Aw),
Omega=list(Units=freq.units),
Gravity=grav,
Model=list(Model="kitagwawa", Pressure=as.pressure),
Response=cbind(omega=omega, wellresp=wellresp),
Response.units=Zunits
)
class(toret) <- "wrsp"
return(toret)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.