Nothing
R/tleaves.R
In tealeaves: Solve for Leaf Temperature Using Energy Balance
Defines functions
tleaves
Documented in
tleaves
#' \code{tleaves}: find leaf temperatures for multiple parameter sets
#'
#' @param leaf_par A list of leaf parameters. This can be generated using the \code{make_leafpar} function.
#'
#' @param enviro_par A list of environmental parameters. This can be generated using the \code{make_enviropar} function.
#'
#' @param constants A list of physical constants. This can be generated using the \code{make_constants} function.
#'
#' @param progress Logical. Should a progress bar be displayed?
#'
#' @param quiet Logical. Should messages be displayed?
#'
#' @param set_units Logical. Should \code{units} be set? The function is faster when FALSE, but input must be in correct units or else results will be incorrect without any warning.
#'
#' @param parallel Logical. Should parallel processing be used via \code{\link[furrr]{future_map}}?
#'
#' @return
#'
#' \code{tleaves}: \cr
#' \cr
#' A tibble with the following \code{units} columns \cr
#'
#' \tabular{ll}{
#' \bold{Input:} \tab \cr
#' \code{abs_l} \tab Absorbtivity of longwave radiation (unitless) \cr
#' \code{abs_s} \tab Absorbtivity of shortwave radiation (unitless) \cr
#' \code{g_sw} \tab Stomatal conductance to H2O (\eqn{\mu}mol H2O / (m^2 s Pa)) \cr
#' \code{g_uw} \tab Cuticular conductance to H2O (\eqn{\mu}mol H2O / (m^2 s Pa)) \cr
#' \code{leafsize} Leaf characteristic dimension \tab (m) \cr
#' \code{logit_sr} \tab Stomatal ratio (logit transformed; unitless) \cr
#' \code{P} \tab Atmospheric pressure (kPa) \cr
#' \code{RH} \tab Relative humidity (unitless) \cr
#' \code{S_lw} \tab incident long-wave radiation flux density (W / m^2) \cr
#' \code{S_sw} \tab incident short-wave (solar) radiation flux density (W / m^2) \cr
#' \code{T_air} \tab Air temperature (K) \cr
#' \code{wind} \tab Wind speed (m / s) \cr
#' \bold{Output:} \tab \cr
#' \code{T_leaf} \tab Equilibrium leaf tempearture (K) \cr
#' \code{value} \tab Leaf energy balance (W / m^2) at \code{tleaf} \cr
#' \code{convergence} \tab Convergence code (0 = converged) \cr
#' \code{R_abs} \tab Total absorbed radiation (W / m^2; see \code{\link{.get_Rabs}}) \cr
#' \code{S_r} \tab Thermal infrared radiation loss (W / m^2; see \code{\link{.get_Sr}}) \cr
#' \code{H} \tab Sensible heat flux density (W / m^2; see \code{\link{.get_H}}) \cr
#' \code{L} \tab Latent heat flux density (W / m^2; see \code{\link{.get_L}}) \cr
#' \code{E} \tab Evapotranspiration (mol H2O/ (m^2 s))
#' }
#'
#' \code{tleaf}: \cr
#' \cr
#' A data.frame with the following numeric columns: \cr
#'
#' \tabular{ll}{
#' \code{T_leaf} \tab Equilibrium leaf temperature (K) \cr
#' \code{value} \tab Leaf energy balance (W / m^2) at \code{tleaf} \cr
#' \code{convergence} \tab Convergence code (0 = converged) \cr
#' \code{R_abs} \tab Total absorbed radiation (W / m^2; see \code{\link{.get_Rabs}}) \cr
#' \code{S_r} \tab Longwave re-radiation (W / m^2; see \code{\link{.get_Sr}}) \cr
#' \code{H} \tab Sensible heat flux density (W / m^2; see \code{\link{.get_H}}) \cr
#' \code{L}\tab Latent heat flux density (W / m^2; see \code{\link{.get_L}}) \cr
#' \code{E} \tab Evapotranspiration (mol H2O/ (m^2 s))
#' }
#'
#' @examples
#' # tleaf for single parameter set:
#'
#' leaf_par <- make_leafpar()
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' tleaf(leaf_par, enviro_par, constants)
#'
#' # tleaves for multiple parameter set:
#'
#' enviro_par <- make_enviropar(
#' replace = list(
#' T_air = set_units(c(293.15, 298.15), K)
#' )
#' )
#' tleaves(leaf_par, enviro_par, constants)
#'
#' @export
#'
tleaves <- function(leaf_par, enviro_par, constants, progress = TRUE,
quiet = FALSE, set_units = TRUE, parallel = FALSE) {
if (set_units) {
leaf_par %<>% tealeaves::leaf_par()
enviro_par %<>% tealeaves::enviro_par()
constants %<>% tealeaves::constants()
} else {
if (!quiet) warning("tleaves: units have not been checked prior to solving")
}
pars <- c(leaf_par, enviro_par)
tsky_function <- NULL
if (is.function(pars$T_sky)) {
tsky_function <- pars$T_sky
pars$T_sky <- NULL
}
par_units <- purrr::map(pars, units) %>%
magrittr::set_names(names(pars))
pars %<>%
purrr::map_if(~ inherits(.x, "units"), drop_units) %>%
make_parameter_sets()
if (!is.null(tsky_function)) {
# For this to work, tsky_function has to work with arg pars
# check with function like check_parameter_function(.f)
# also, this seems pretty inefficient, but I don't have a better solution
tidyr::crossing(
par = names(pars[[1]]),
i = 1:length(pars)
) %>%
purrr::pwalk(~ {
glue::glue("units(pars[[{i}]]${par}) <<- par_units${par}",
par = ..1, i = ..2) %>%
parse(text = .) %>%
eval()
})
pars <- pars %>%
purrr::map(~ {
.x$T_sky <- tsky_function(.x)
.x
})
par_units$T_sky <- units(pars[[1]]$T_sky)
pars <- pars %>%
purrr::map(~ {purrr::map_if(.x, ~ inherits(.x, "units"), drop_units)})
}
soln <- find_tleaves(pars, constants, progress, quiet, parallel)
pars %<>% purrr::map_dfr(purrr::flatten_dfr)
colnames(pars) %>%
glue::glue("units(pars${x}) <<- par_units${x}", x = .) %>%
parse(text = .) %>%
eval()
pars %<>% dplyr::bind_cols(soln)
pars$T_leaf %<>% set_units(K)
pars$E %<>% set_units(mol / m ^ 2 / s)
pars$R_abs %<>% set_units(W / m ^ 2)
pars$S_r %<>% set_units(W / m ^ 2)
pars$H %<>% set_units(W / m ^ 2)
pars$L %<>% set_units(W / m ^ 2)
pars
}
#' \code{tleaf}: find leaf temperatures for a single parameter set
#' @rdname tleaves
#' @export
tleaf <- function(leaf_par, enviro_par, constants, quiet = FALSE,
set_units = TRUE) {
if (set_units) {
leaf_par %<>% tealeaves::leaf_par()
enviro_par %<>% tealeaves::enviro_par()
constants %<>% tealeaves::constants()
} else {
if (!quiet) warning("tleaf: units have not been checked prior to solving")
}
if (is.function(enviro_par$T_sky)) {
enviro_par$T_sky <- enviro_par$T_sky(enviro_par)
}
# Drop units for solving
ulp <- purrr::map_if(leaf_par, ~ inherits(.x, "units"), drop_units)
uep <- purrr::map_if(enviro_par, ~ inherits(.x, "units"), drop_units)
ucs <- purrr::map_if(constants, ~ inherits(.x, "units"), drop_units)
soln <- find_tleaf(ulp, uep, ucs, quiet)
# Check results -----
if (soln$convergence == 1) {
"stats::uniroot did not converge, NA returned. Inspect parameters carefully." %>%
crayon::red() %>%
message()
}
# Energy balance components -----
components <- suppressWarnings(
if (is.na(soln$T_leaf)) {
list(R_abs = NA, S_r = NA, H = NA, L = NA, E = NA)
} else {
soln$T_leaf %>%
energy_balance(leaf_par = ulp, enviro_par = uep,
constants = ucs, quiet = TRUE, components = TRUE,
set_units = FALSE) %>%
magrittr::use_series("components")
}
)
soln %<>% dplyr::bind_cols(components)
if (is.na(soln$T_leaf)) {
soln %<>% dplyr::bind_cols(data.frame(Ar = NA, Gr = NA, Re = NA))
} else {
soln %<>% dplyr::bind_cols(Ar(soln$T_leaf, c(ulp, uep, ucs),
unitless = TRUE))
}
soln$T_leaf %<>% set_units(K)
# Adding g_bw to output
if (is.na(soln$T_leaf)) {
soln$g_bw <- NA
} else {
soln$g_bw <- .get_gbw(
T_leaf = drop_units(soln$T_leaf),
surface = "lower",
pars = c(ulp, uep, ucs),
unitless = TRUE
) +
.get_gbw(
T_leaf = drop_units(soln$T_leaf),
surface = "upper",
pars = c(ulp, uep, ucs),
unitless = TRUE
)
}
# Return -----
soln
}
#' Calculate leaf energy balance
#'
#' @inheritParams tleaves
#'
#' @param tleaf Leaf temperature in Kelvin. If input is numeric, it will be automatically converted to \code{units}.
#'
#' @param quiet Logical. Should a message appear about conversion from \code{numeric} to \code{units}? Useful for finding leaf temperature that balances heat transfer using \code{\link[stats]{uniroot}}.
#'
#' @param components Logical. Should leaf energy components be returned? Transpiration (in mol / (m^2 s)) also returned.
#'
#' @return A numeric value in W / m^2. Optionally, a named list of energy balance components in W / m^2 and transpiration in mol / (m^2 s).
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#' ep$T_sky <- ep$T_sky(ep)
#'
#' T_leaf <- set_units(298.15, K)
#'
#' energy_balance(T_leaf, lp, ep, cs, FALSE, TRUE, TRUE)
#'
#' @export
#'
energy_balance <- function(tleaf, leaf_par, enviro_par, constants,
quiet = FALSE, components = FALSE, set_units = FALSE) {
# Checks -----
stopifnot(length(set_units) == 1L & is.logical(set_units))
if (set_units) {
tleaf %<>%
set_units(K) %>%
drop_units()
leaf_par %<>% leaf_par()
enviro_par %<>% enviro_par()
constants %<>% constants()
}
stopifnot(length(quiet) == 1L & is.logical(quiet))
stopifnot(length(components) == 1L & is.logical(components))
pars <- c(leaf_par, enviro_par, constants) %>%
purrr::map_if(~ inherits(.x, "units"), drop_units)
# R_abs: total absorbed radiation (W m^-2) -----
R_abs <- .get_Rabs(pars, unitless = TRUE)
# S_r: longwave re-radiation (W m^-2) -----
S_r <- .get_Sr(tleaf, pars)
# H: sensible heat flux density (W m^-2) -----
H <- .get_H(tleaf, pars, unitless = TRUE)
# L: latent heat flux density (W m^-2) -----
L <- .get_L(tleaf, pars, unitless = TRUE)
# Return -----
if (components) {
# E: transpiration (mol / (m^2 s))
E <- .get_gtw(tleaf, pars, unitless = TRUE) *
.get_dwv(tleaf, pars, unitless = TRUE)
ret <- list(
energy_balance = R_abs - (S_r + H + L),
components = list(R_abs = R_abs, S_r = S_r, H = H, L = L, E = E)
)
ret$energy_balance %<>% set_units(W / m ^ 2)
ret$components$R_abs %<>% set_units(W / m ^ 2)
ret$components$S_r %<>% set_units(W / m ^ 2)
ret$components$H %<>% set_units(W / m ^ 2)
ret$components$L %<>% set_units(W / m ^ 2)
ret$components$E %<>% set_units(mol / m ^ 2 / s)
return(ret)
}
R_abs - (S_r + H + L)
}
#' R_abs: total absorbed radiation (W / m^2)
#'
#' @param pars Concatenated parameters (\code{leaf_par}, \code{enviro_par}, and \code{constants})
#'
#' @param unitless Logical. Should function use parameters with \code{units}? The function is faster when FALSE, but input must be in correct units or else results will be incorrect without any warning.
#'
#' @inheritParams tleaves
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' The following treatment follows Okajima et al. (2012):
#'
#' \deqn{R_\mathrm{abs} = \alpha_\mathrm{s} (1 + r) S_\mathrm{sw} + \alpha_\mathrm{l} \sigma (T_\mathrm{sky} ^ 4 + T_\mathrm{air} ^ 4)}{R_abs = \alpha_s (1 + r) S_sw + \alpha_l \sigma (T_sky ^ 4 + T_air ^ 4)}
#'
#' The incident longwave (aka thermal infrared) radiation is modeled from sky and air temperature \eqn{\sigma (T_\mathrm{sky} ^ 4 + T_\mathrm{air} ^ 4)}{\sigma (T_sky ^ 4 + T_air ^ 4)} where \eqn{T_\mathrm{sky}}{T_sky} is function of the air temperature and incoming solar shortwave radiation:
#'
#' \deqn{T_\mathrm{sky} = T_\mathrm{air} - 20 S_\mathrm{sw} / 1000}{T_sky = T_air - 20 S_sw / 1000}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\alpha_\mathrm{s}}{\alpha_s} \tab \code{abs_s} \tab absorbtivity of shortwave radiation (0.3 - 4 \eqn{\mu}m) \tab none \tab 0.80\cr
#' \eqn{\alpha_\mathrm{l}}{\alpha_l} \tab \code{abs_l} \tab absorbtivity of longwave radiation (4 - 80 \eqn{\mu}m) \tab none \tab 0.97\cr
#' \eqn{r} \tab \code{r} \tab reflectance for shortwave irradiance (albedo) \tab none \tab 0.2 \cr
#' \eqn{\sigma} \tab \code{s} \tab Stefan-Boltzmann constant \tab W / (m\eqn{^2} K\eqn{^4}) \tab 5.67e-08 \cr
#' \eqn{S_\mathrm{sw}}{S_sw} \tab \code{S_sw} \tab incident short-wave (solar) radiation flux density \tab W / m\eqn{^2} \tab 1000 \cr
#' \eqn{S_\mathrm{lw}}{S_lw} \tab \code{S_lw} \tab incident long-wave radiation flux density \tab W / m\eqn{^2} \tab calculated \cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15 \cr
#' \eqn{T_\mathrm{sky}}{T_sky} \tab \code{T_sky} \tab sky temperature \tab K \tab calculated
#' }
#'
#' @references
#'
#' Okajima Y, H Taneda, K Noguchi, I Terashima. 2012. Optimum leaf size predicted by a novel leaf energy balance model incorporating dependencies of photosynthesis on light and temperature. Ecological Research 27: 333-46.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#' ep$T_sky <- ep$T_sky(ep)
#'
#' tealeaves:::.get_Rabs(c(cs, ep, lp), FALSE)
#'
.get_Rabs <- function(pars, unitless) {
if (unitless) {
R_abs <- pars$abs_s * (1 + pars$r) * pars$S_sw + pars$abs_l * pars$s * (pars$T_sky ^ 4 + pars$T_air ^ 4)
} else {
R_abs <- pars$abs_s * (set_units(1) + pars$r) * pars$S_sw + pars$abs_l * pars$s * (pars$T_sky ^ 4 + pars$T_air ^ 4)
R_abs %<>% set_units(W / m ^ 2)
}
R_abs
}
#' S_r: longwave re-radiation (W / m^2)
#'
#' @inheritParams .get_Rabs
#' @param T_leaf Leaf temperature in Kelvin
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' \deqn{S_\mathrm{r} = 2 \sigma \alpha_\mathrm{l} T_\mathrm{air} ^ 4}{S_r = 2 \sigma \alpha_l T_air ^ 4 }
#'
#' The factor of 2 accounts for re-radiation from both leaf surfaces (Foster and Smith 1986). \cr
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\alpha_\mathrm{l}}{\alpha_l} \tab \code{abs_l} \tab absorbtivity of longwave radiation (4 - 80 \eqn{\mu}m) \tab none \tab 0.97\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{\sigma} \tab \code{s} \tab Stefan-Boltzmann constant \tab W / (m\eqn{^2} K\eqn{^4}) \tab 5.67e-08
#' }
#'
#' Note that leaf absorbtivity is the same value as leaf emissivity
#'
#' @references
#'
#' Foster JR, Smith WK. 1986. Influence of stomatal distribution on transpiration in low-wind environments. Plant, Cell & Environment 9: 751-9.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_Sr(T_leaf, c(cs, ep, lp))
#'
.get_Sr <- function(T_leaf, pars) 2 * pars$s * pars$abs_l * T_leaf ^ 4
#' H: sensible heat flux density (W / m^2)
#'
#' @inheritParams tleaves
#' @inheritParams .get_Rabs
#' @inheritParams .get_Sr
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' \deqn{H = P_\mathrm{a} c_p g_\mathrm{h} (T_\mathrm{leaf} - T_\mathrm{air})}{H = P_a c_p g_h * (T_leaf - T_air)}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{c_p} \tab \code{c_p} \tab heat capacity of air \tab J / (g K) \tab 1.01\cr
#' \eqn{g_\mathrm{h}}{g_h} \tab \code{g_h} \tab boundary layer conductance to heat \tab m / s \tab \link[=.get_gh]{calculated}\cr
#' \eqn{P_\mathrm{a}}{P_a} \tab \code{P_a} \tab density of dry air \tab g / m^3 \tab \link[=.get_Pa]{calculated}\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @seealso \code{\link{.get_gh}}, \code{\link{.get_Pa}}
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_H(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_H <- function(T_leaf, pars, unitless) {
# Density of dry air
P_a <- .get_Pa(T_leaf, pars, unitless)
# Boundary layer conductance to heat
g_h <- sum(.get_gh(T_leaf, "lower", pars, unitless),
.get_gh(T_leaf, "upper", pars, unitless))
H <- P_a * pars$c_p * g_h * (T_leaf - pars$T_air)
H
}
#' P_a: density of dry air (g / m^3)
#'
#' @inheritParams .get_H
#'
#' @return Value in g / m\eqn{^3} of class \code{units}
#'
#' @details
#'
#' \deqn{P_\mathrm{a} = P / (R_\mathrm{air} (T_\mathrm{leaf} - T_\mathrm{air}) / 2)}{P_a = P / (R_air (T_leaf - T_air) / 2)}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{P} \tab \code{P} \tab atmospheric pressure \tab kPa \tab 101.3246\cr
#' \eqn{R_\mathrm{air}}{R_air} \tab \code{R_air} \tab specific gas constant for dry air \tab J / (kg K) \tab 287.058\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_Pa(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_Pa <- function(T_leaf, pars, unitless) {
P_a <- pars$P / (pars$R_air * (pars$T_air + T_leaf) / 2)
if (unitless) {
P_a %<>% magrittr::multiply_by(1e6)
} else {
P_a %<>% set_units(g / m ^ 3)
}
P_a
}
#' g_h: boundary layer conductance to heat (m / s)
#'
#' @inheritParams .get_H
#' @param surface Leaf surface (lower or upper)
#'
#' @return Value in m/s of class \code{units}
#'
#' @details
#'
#' \deqn{g_\mathrm{h} = D_\mathrm{h} Nu / d}{g_h = D_h Nu / d}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{h}}{D_h} \tab \code{D_h} \tab diffusion coefficient for heat in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{Nu} \tab \code{Nu} \tab Nusselt number \tab none \tab \link[=.get_nu]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gh(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_gh <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
# Calculate diffusion coefficient to heat
D_h <- .get_Dx(pars$D_h0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
# Calculate Nusselt numbers
Nu <- .get_nu(T_leaf, surface, pars, unitless)
D_h * Nu / pars$leafsize
}
#' D_x: Calculate diffusion coefficient for a given temperature and pressure
#'
#' @param D_0 Diffusion coefficient at 273.15 K (0 °C) and 101.3246 kPa
#' @param Temp Temperature in Kelvin
#' @param eT Exponent for temperature dependence of diffusion
#' @param P Atmospheric pressure in kPa
#' @inheritParams .get_Rabs
#' @inheritParams tleaves
#'
#' @return Value in m\eqn{^2}/s of class \code{units}
#'
#' @details
#'
#' \deqn{D = D_\mathrm{0} (T / 273.15) ^ {eT} (101.3246 / P)}{D = D_0 [(T / 273.15) ^ eT] (101.3246 / P)}
#' \cr
#' According to Montieth & Unger (2013), eT is generally between 1.5 and 2. Their data in Appendix 3 indicate \eqn{eT = 1.75} is reasonable for environmental physics.
#'
#' @references
#'
#' Monteith JL, Unsworth MH. 2013. Principles of Environmental Physics. 4th edition. Academic Press, London.
#'
#' @examples
#'
#' tealeaves:::.get_Dx(
#' D_0 = set_units(2.12e-05, m^2/s),
#' Temp = set_units(298.15, K),
#' eT = set_units(1.75),
#' P = set_units(101.3246, kPa),
#' unitless = FALSE
#' )
#'
.get_Dx <- function(D_0, Temp, eT, P, unitless) {
# See Eq. 3.10 in Monteith & Unger ed. 4
if (unitless) {
Dx <- D_0 * (Temp / 273.15) ^ eT * (101.3246 / P)
} else {
Dx <- D_0 *
drop_units((set_units(Temp, K) / set_units(273.15, K))) ^ drop_units(eT) *
drop_units((set_units(101.3246, kPa) / set_units(P, kPa)))
}
Dx
}
#' Gr: Grashof number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' \deqn{Gr = t_\mathrm{air} G d ^ 3 |T_\mathrm{v,leaf} - T_\mathrm{v,air}| / D_\mathrm{m} ^ 2}{Gr = t_air G d ^ 3 abs(Tv_leaf - Tv_air) / D_m ^ 2}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{m}}{D_m} \tab \code{D_m} \tab diffusion coefficient of momentum in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{G} \tab \code{G} \tab gravitational acceleration \tab m / s\eqn{^2} \tab 9.8\cr
#' \eqn{t_\mathrm{air}}{t_air} \tab \code{t_air} \tab coefficient of thermal expansion of air \tab 1 / K \tab 1 / Temp \cr
#' \eqn{T_\mathrm{v,air}}{Tv_air} \tab \code{Tv_air} \tab virtual air temperature \tab K \tab \link[=.get_Tv]{calculated}\cr
#' \eqn{T_\mathrm{v,leaf}}{Tv_leaf} \tab \code{Tv_leaf} \tab virtual leaf temperature \tab K \tab \link[=.get_Tv]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gr(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_gr <- function(T_leaf, pars, unitless) {
# Calculate virtual temperature
# Assumes inside of leaf is 100% RH
Tv_leaf <- .get_Tv(T_leaf, .get_ps(T_leaf, pars$P, unitless), pars$P,
pars$epsilon, unitless)
Tv_air <- .get_Tv(pars$T_air, pars$RH * .get_ps(pars$T_air, pars$P, unitless), pars$P,
pars$epsilon, unitless)
D_m <- .get_Dx(pars$D_m0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
if (unitless) {
Gr <- (1 / pars$T_air) * pars$G * pars$leafsize ^ 3 *
abs(Tv_leaf - Tv_air) / D_m ^ 2
} else {
Gr <- (set_units(1) / pars$T_air) * pars$G * pars$leafsize ^ 3 *
abs(Tv_leaf - Tv_air) / D_m ^ 2
}
Gr
}
#' Calculate virtual temperature
#'
#' @inheritParams .get_Dx
#' @param p water vapour pressure in kPa
#' @param epsilon ratio of water to air molar masses (unitless)
#'
#' @return Value in K of class \code{units}
#'
#' @details
#'
#' \deqn{T_\mathrm{v} = T / [1 - (1 - \epsilon) (p / P)]}{T_v = T / [1 - (1 - epsilon) (p / P)]}
#'
#' Eq. 2.35 in Monteith & Unsworth (2013) \cr
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\epsilon} \tab \code{epsilon} \tab ratio of water to air molar masses \tab unitless \tab 0.622 \cr
#' \eqn{p} \tab \code{p} \tab water vapour pressure \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{P} \tab \code{P} \tab atmospheric pressure \tab kPa \tab 101.3246
#' }
#'
#' @references
#'
#' Monteith JL, Unsworth MH. 2013. Principles of Environmental Physics. 4th edition. Academic Press, London.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#' p <- ep$RH * tealeaves:::.get_ps(T_leaf, ep$P, FALSE)
#' tealeaves:::.get_Tv(T_leaf, p, ep$P, cs$epsilon, FALSE)
#'
.get_Tv <- function(Temp, p, P, epsilon, unitless) {
if (unitless) {
Tv <- Temp / (1 - (1 - epsilon) * (p / P))
} else {
Tv <- set_units(Temp, K) /
(set_units(1) - (set_units(1) - epsilon) *
(set_units(p, kPa) / set_units(P, kPa)))
}
Tv
}
#' Saturation water vapour pressure (kPa)
#'
#' @inheritParams .get_Dx
#'
#' @return Value in kPa of class \code{units}
#'
#' @details
#'
#' Goff-Gratch equation (see http://cires1.colorado.edu/~voemel/vp.html) \cr
#' \cr
#' This equation assumes P = 1 atm = 101.3246 kPa, otherwise boiling temperature needs to change \cr
#'
#' @references \url{http://cires1.colorado.edu/~voemel/vp.html}
#'
#' @examples
#'
#' T_leaf <- set_units(298.15, K)
#' P <- set_units(101.3246, kPa)
#' tealeaves:::.get_ps(T_leaf, P, FALSE)
#'
.get_ps <- function(Temp, P, unitless) {
# Goff-Gratch equation (see http://cires1.colorado.edu/~voemel/vp.html)
# This assumes P = 1 atm = 101.3246 kPa, otherwise boiling temperature needs to change
# This returns p_s in hPa
if (unitless) {
P %<>% magrittr::multiply_by(10)
} else {
Temp %<>% set_units(K) %>% drop_units()
P %<>% set_units(hPa) %>% drop_units()
}
p_s <- 10 ^ (-7.90298 * (373.16 / Temp - 1) +
5.02808 * log10(373.16 / Temp) -
1.3816e-7 * (10 ^ (11.344 * (1 - Temp / 373.16) - 1)) +
8.1328e-3 * (10 ^ (-3.49149 * (373.16 / Temp - 1)) - 1) +
log10(P))
# Convert to kPa
if (unitless) {
p_s %<>% magrittr::multiply_by(0.1)
} else {
p_s %<>%
set_units(hPa) %>%
set_units(kPa)
}
p_s
}
#' Re: Reynolds number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' \deqn{Re = u d / D_\mathrm{m}}{Re = u d / D_m}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{m}}{D_m} \tab \code{D_m} \tab diffusion coefficient of momentum in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{u} \tab \code{wind} \tab windspeed \tab m / s \tab 2
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_re(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_re <- function(T_leaf, pars, unitless) {
D_m <- .get_Dx(pars$D_m0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
Re <- pars$wind * pars$leafsize / D_m
Re
}
#' Nu: Nusselt number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' The Nusselt number depends on a combination how much free or forced convection predominates. For mixed convection: \cr
#' \cr
#' \deqn{Nu = (a Re ^ b) ^ {3.5} + (c Gr ^ d) ^ {3.5}) ^ {1 / 3.5}}{Nu = (a Re ^ b) ^ 3.5 + (c Gr ^ d) ^ 3.5) ^ (1 / 3.5)}
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{a, b, c, d} \tab \code{a, b, c, d} \tab empirical coefficients \tab none \tab \link[=make_constants]{calculated}\cr
#' \eqn{Gr} \tab \code{Gr} \tab Grashof number \tab none \tab \link[=.get_gr]{calculated}\cr
#' \eqn{Re} \tab \code{Re} \tab Reynolds number \tab none \tab \link[=.get_re]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_nu(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_nu <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
Gr <- .get_gr(T_leaf, pars, unitless)
Re <- .get_re(T_leaf, pars, unitless)
# Archemides number
# Ar <- Gr / Re ^ 2
cons <- pars$nu_constant(Re, "forced", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_forced <- cons$a * Re ^ cons$b
} else {
Nu_forced <- cons$a * drop_units(Re) ^ cons$b
}
cons <- pars$nu_constant(Re, "free", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_free <- cons$a * Gr ^ cons$b
} else {
Nu_free <- cons$a * drop_units(Gr) ^ cons$b
}
Nu <- (Nu_forced ^ 3.5 + Nu_free ^ 3.5) ^ (1 / 3.5)
if (!unitless) Nu %<>% set_units()
Nu
}
#' L: Latent heat flux density (W / m^2)
#'
#' @inheritParams .get_H
#'
#' @return Value in W / m^2 of class \code{units}
#'
#' @details
#'
#' \deqn{L = h_\mathrm{vap} g_\mathrm{tw} d_\mathrm{wv}}{L = h_vap g_tw d_wv}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d_\mathrm{wv}}{d_wv} \tab \code{d_wv} \tab water vapour gradient \tab mol / m ^ 3 \tab \link[=.get_dwv]{calculated} \cr
#' \eqn{h_\mathrm{vap}}{h_vap} \tab \code{h_vap} \tab latent heat of vaporization \tab J / mol \tab \link[=.get_hvap]{calculated} \cr
#' \eqn{g_\mathrm{tw}}{g_tw} \tab \code{g_tw} \tab total conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab \link[=.get_gtw]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_L(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_L <- function(T_leaf, pars, unitless) {
h_vap <- .get_hvap(T_leaf, unitless)
g_tw <- .get_gtw(T_leaf, pars, unitless)
d_wv <- .get_dwv(T_leaf, pars, unitless)
L <- h_vap * g_tw * d_wv
if (!unitless) L %<>% set_units(W / m ^ 2)
L
}
#' d_wv: water vapour gradient (mol / m ^ 3)
#'
#' @inheritParams .get_H
#'
#' @return Value in mol / m^3 of class \code{units}
#'
#' @details
#'
#' \bold{Water vapour gradient:} The water vapour pressure differential from inside to outside of the leaf is the saturation water vapor pressure inside the leaf (p_leaf) minus the water vapor pressure of the air (p_air):
#'
#' \deqn{d_\mathrm{wv} = p_\mathrm{leaf} / (R T_\mathrm{leaf}) - RH p_\mathrm{air} / (R T_\mathrm{air})}{d_wv = p_leaf / (R T_leaf) - RH p_air / (R T_air)}
#'
#' Note that water vapor pressure is converted from kPa to mol / m^3 using ideal gas law. \cr
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{p_\mathrm{air}}{p_air} \tab \code{p_air} \tab saturation water vapour pressure of air \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{p_\mathrm{leaf}}{p_leaf} \tab \code{p_leaf} \tab saturation water vapour pressure inside the leaf \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{R} \tab \code{R} \tab ideal gas constant \tab J / (mol K) \tab 8.3144598\cr
#' \eqn{\mathrm{RH}}{RH} \tab \code{RH} \tab relative humidity \tab \% \tab 0.50\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' # Water vapour gradient:
#'
#' leaf_par <- make_leafpar()
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' pars <- c(leaf_par, enviro_par, constants)
#' T_leaf <- set_units(300, K)
#' T_air <- set_units(298.15, K)
#' p_leaf <- set_units(35.31683, kPa)
#' p_air <- set_units(31.65367, kPa)
#'
#' d_wv <- p_leaf / (pars$R * T_leaf) - pars$RH * p_air / (pars$R * T_air)
#'
.get_dwv <- function(T_leaf, pars, unitless) {
# Water vapour differential converted from kPa to mol m ^ -3 using ideal gas law
d_wv <- .get_ps(T_leaf, pars$P, unitless) / (pars$R * T_leaf) -
pars$RH * .get_ps(pars$T_air, pars$P, unitless) / (pars$R * pars$T_air)
if (unitless) {
d_wv %<>% magrittr::multiply_by(1e3)
} else {
d_wv %<>% set_units(mol / m ^ 3)
}
d_wv
}
#' g_tw: total conductance to water vapour (m/s)
#'
#' @inheritParams .get_H
#'
#' @return Value in m/s of class \code{units}
#'
#' @details
#'
#' \bold{Total conductance to water vapor:} The total conductance to water vapor (\eqn{g_\mathrm{tw}}{g_tw}) is the sum of the parallel lower (abaxial) and upper (adaxial) conductances:
#'
#' \deqn{g_\mathrm{tw} = g_\mathrm{w,lower} + g_\mathrm{w,upper}}{g_tw = gw_lower + gw_upper}
#'
#' The conductance to water vapor on each surface is a function of parallel stomatal (\eqn{g_\mathrm{sw}}{g_sw}) and cuticular (\eqn{g_\mathrm{uw}}{g_uw}) conductances in series with the boundary layer conductance (\eqn{g_\mathrm{bw}}{g_bw}). The stomatal, cuticular, and boundary layer conductance on the lower surface are:
#'
#' \deqn{g_\mathrm{sw,lower} = g_\mathrm{sw} (1 - sr) R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{gsw_lower = g_sw (1 - sr) R (T_leaf + T_air) / 2}
#' \deqn{g_\mathrm{uw,lower} = g_\mathrm{uw} / 2 R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{guw_lower = g_uw / 2 R (T_leaf + T_air) / 2}
#' \cr
#' See \code{\link{.get_gbw}} for details on calculating boundary layer conductance. The equations for the upper surface are:
#'
#' \deqn{g_\mathrm{sw,upper} = g_\mathrm{sw} sr R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{gsw_upper = g_sw sr R (T_leaf + T_air) / 2}
#' \deqn{g_\mathrm{uw,upper} = g_\mathrm{uw} / 2 R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{guw_upper = g_uw / 2 R (T_leaf + T_air) / 2}
#' \cr
#' Note that the stomatal and cuticular conductances are given in units of (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) (see \code{\link{make_leafpar}}) and converted to m/s using the ideal gas law. The total leaf stomatal (\eqn{g_\mathrm{sw}}{g_sw}) and cuticular (\eqn{g_\mathrm{uw}}{g_uw}) conductances are partitioned across lower and upper surfaces. The stomatal conductance on each surface depends on stomatal ratio (sr); the cuticular conductance is assumed identical on both surfaces.
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{g_\mathrm{sw}}{g_sw} \tab \code{g_sw} \tab stomatal conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab 5\cr
#' \eqn{g_\mathrm{uw}}{g_uw} \tab \code{g_uw} \tab cuticular conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab 0.1\cr
#' \eqn{R} \tab \code{R} \tab ideal gas constant \tab J / (mol K) \tab 8.3144598\cr
#' \eqn{\mathrm{logit}(sr)}{logit(sr)} \tab \code{logit_sr} \tab stomatal ratio (logit transformed) \tab none \tab 0 = logit(0.5)\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' # Total conductance to water vapor
#'
#' ## Hypostomatous leaf; default parameters
#' leaf_par <- make_leafpar(replace = list(logit_sr = set_units(-Inf)))
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' pars <- c(leaf_par, enviro_par, constants)
#' T_leaf <- set_units(300, K)
#'
#' ## Fixing boundary layer conductance rather than calculating
#' gbw_lower <- set_units(0.1, m / s)
#' gbw_upper <- set_units(0.1, m / s)
#'
#' # Lower surface ----
#' ## Note that pars$logit_sr is logit-transformed! Use stats::plogis() to convert to proportion.
#' gsw_lower <- set_units(pars$g_sw * (set_units(1) - stats::plogis(pars$logit_sr)) * pars$R *
#' ((T_leaf + pars$T_air) / 2), "m / s")
#' guw_lower <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
#' gtw_lower <- 1 / (1 / (gsw_lower + guw_lower) + 1 / gbw_lower)
#'
#' # Upper surface ----
#' gsw_upper <- set_units(pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
#' ((T_leaf + pars$T_air) / 2), m / s)
#' guw_upper <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
#' gtw_upper <- 1 / (1 / (gsw_upper + guw_upper) + 1 / gbw_upper)
#'
#' ## Lower and upper surface are in parallel
#' g_tw <- gtw_lower + gtw_upper
#'
.get_gtw <- function(T_leaf, pars, unitless) {
# Lower surface ----
gbw_lower <- .get_gbw(T_leaf, "lower", pars, unitless)
# Convert stomatal and cuticular conductance from molar to 'engineering' units
# See email from Tom Buckley (July 4, 2017)
if (unitless) {
gsw_lower <- (pars$g_sw * (1 - stats::plogis(pars$logit_sr)) * pars$R *
((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
guw_lower <- (pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
} else {
gsw_lower <- set_units(pars$g_sw * (set_units(1) - stats::plogis(pars$logit_sr)) *
pars$R * ((T_leaf + pars$T_air) / 2), m / s)
guw_lower <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
}
gtw_lower <- 1 / (1 / (gsw_lower + guw_lower) + 1 / gbw_lower)
# Upper surface ----
gbw_upper <- .get_gbw(T_leaf, "upper", pars, unitless)
# Convert stomatal and cuticular conductance from molar to 'engineering' units
# See email from Tom Buckley (July 4, 2017)
if (unitless) {
gsw_upper <- (pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
guw_upper <- (pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
} else {
gsw_upper <- set_units(pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
((T_leaf + pars$T_air) / 2), m / s)
guw_upper <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
}
gtw_upper <- 1 / (1 / (gsw_upper + guw_upper) + 1 / gbw_upper)
# Lower and upper surface are in parallel
g_tw <- gtw_lower + gtw_upper
g_tw
}
#' h_vap: heat of vaporization (J / mol)
#'
#' @inheritParams .get_H
#'
#' @return Value in J/mol of class \code{units}
#'
#' @details
#'
#' \bold{Heat of vaporization:} The heat of vaporization (\eqn{h_\mathrm{vap}}{h_vap}) is a function of temperature. I used data from on temperature and \eqn{h_\mathrm{vap}}{h_vap} from Nobel (2009, Appendix 1) to estimate a linear regression. See Examples.
#'
#' @examples
#'
#' # Heat of vaporization and temperature
#' ## data from Nobel (2009)
#'
#' T_K <- 273.15 + c(0, 10, 20, 25, 30, 40, 50, 60)
#' h_vap <- 1e3 * c(45.06, 44.63, 44.21, 44.00,
#' 43.78, 43.35, 42.91, 42.47) # (in J / mol)
#'
#' fit <- lm(h_vap ~ T_K)
#'
#' ## coefficients are 56847.68250 J / mol and 43.12514 J / (mol K)
#'
#' coef(fit)
#'
#' T_leaf <- 298.15
#' h_vap <- set_units(56847.68250, J / mol) -
#' set_units(43.12514, J / mol / K) * set_units(T_leaf, K)
#'
#' ## h_vap at 298.15 K is 43989.92 [J/mol]
#'
#' set_units(h_vap, J / mol)
#'
#' tealeaves:::.get_hvap(set_units(298.15, K), FALSE)
#'
#' @references Nobel PS. 2009. Physicochemical and Environmental Plant Physiology. 4th Edition. Academic Press.
#'
.get_hvap <- function(T_leaf, unitless) {
# Equation from Foster and Smith 1986 seems to be off:
# h_vap <- 4.504e4 - 41.94 * T_leaf
# Instead, using regression based on data from Nobel (2009, 4th Ed, Appendix 1)
# T_K <- 273.15 + c(0, 10, 20, 25, 30, 40, 50, 60)
# h_vap <- 1e3 * c(45.06, 44.63, 44.21, 44, 43.78, 43.35, 42.91, 42.47) # (in J / mol)
# fit <- lm(h_vap ~ T_K)
if (unitless) {
h_vap <- 56847.68250 - 43.12514 * T_leaf
} else {
h_vap <- set_units(56847.68250, J / mol) -
set_units(43.12514, J / mol / K) * set_units(T_leaf, K)
h_vap %<>% set_units(J / mol)
}
h_vap
}
#' g_bw: Boundary layer conductance to water vapour (m / s)
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return Value in m / s of class \code{units}
#'
#' @details
#'
#' \deqn{g_\mathrm{bw} = D_\mathrm{w} Sh / d}{g_bw = D_w Sh / d}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{w}}{D_w} \tab \code{D_w} \tab diffusion coefficient for water vapour \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{Sh} \tab \code{Sh} \tab Sherwood number \tab none \tab \link[=.get_sh]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gbw(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_gbw <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
D_w <- .get_Dx(pars$D_w0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
# Calculate Sherwood numbers
Sh <- .get_sh(T_leaf, surface, pars, unitless)
D_w * Sh / pars$leafsize
}
#' Sh: Sherwood number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' The Sherwood number depends on a combination how much free or forced convection predominates. For mixed convection: \cr
#' \cr
#' \deqn{Sh = (a Re ^ b) ^ {3.5} + (c Gr ^ d) ^ {3.5}) ^ {1 / 3.5}}{Sh = (a Re ^ b) ^ 3.5 + (c Gr ^ d) ^ 3.5) ^ (1 / 3.5)}
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{a, b, c, d} \tab \code{a, b, c, d} \tab empirical coefficients \tab none \tab \link[=make_constants]{calculated}\cr
#' \eqn{Gr} \tab \code{Gr} \tab Grashof number \tab none \tab \link[=.get_gr]{calculated}\cr
#' \eqn{Re} \tab \code{Re} \tab Reynolds number \tab none \tab \link[=.get_re]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_sh(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_sh <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
Gr <- .get_gr(T_leaf, pars, unitless)
Re <- .get_re(T_leaf, pars, unitless)
# Archemides number
# Ar <- Gr / Re ^ 2
D_h <- .get_Dx(pars$D_h0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
D_w <- .get_Dx(pars$D_w0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
cons <- pars$nu_constant(Re, "forced", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_forced <- cons$a * Re ^ cons$b
Sh_forced <- Nu_forced * (D_h / D_w) ^ pars$sh_constant("forced")
} else {
Nu_forced <- cons$a * drop_units(Re) ^ cons$b
Sh_forced <- Nu_forced * drop_units(D_h / D_w) ^ pars$sh_constant("forced")
}
cons <- pars$nu_constant(Re, "free", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_free <- cons$a * Gr ^ cons$b
Sh_free <- Nu_free * (D_h / D_w) ^ pars$sh_constant("free")
} else {
Nu_free <- cons$a * drop_units(Gr) ^ cons$b
Sh_free <- Nu_free * drop_units(D_h / D_w) ^ pars$sh_constant("free")
}
Sh <- (Sh_forced ^ 3.5 + Sh_free ^ 3.5) ^ (1 / 3.5)
if (!unitless) Sh %<>% set_units()
Sh
}
make_parameter_sets <- function(pars) {
pars %<>%
names() %>%
glue::glue("{x} = pars${x}", x = .) %>%
stringr::str_c(collapse = ", ") %>%
glue::glue("tidyr::crossing({x})", x = .) %>%
parse(text = .) %>%
eval() %>%
purrr::transpose()
pars
}
find_tleaves <- function(par_sets, constants, progress, quiet, parallel) {
# This function is not intended to be called directly by users.
# It assumes all parameters have correct units
if (!quiet) {
glue::glue("\nSolving for T_leaf from {n} parameter set{s}...",
n = length(par_sets),
s = dplyr::if_else(length(par_sets) > 1, "s", "")) %>%
crayon::green() %>%
message(appendLF = FALSE)
}
if (parallel) future::plan("multisession")
if (progress & !parallel) pb <- dplyr::progress_estimated(length(par_sets))
soln <- suppressWarnings(
par_sets %>%
furrr::future_map_dfr(~{
ret <- tleaf(.x, .x, constants, quiet = TRUE, set_units = FALSE)
if (progress & !parallel) pb$tick()$print()
ret
}, .progress = progress)
)
soln
}
find_tleaf <- function(leaf_par, enviro_par, constants, quiet) {
# This function is not intended to be called directly by users.
# It assumes all parameters have correct units
# Balance energy fluxes -----
if (!quiet) {
"\nSolving for T_leaf ..." %>%
crayon::green() %>%
message(appendLF = FALSE)
}
.f <- function(T_leaf, ...) {
eb <- energy_balance(T_leaf, ...)
eb
}
fit <- safely_uniroot(
f = .f,
leaf_par = leaf_par, enviro_par = enviro_par, constants = constants,
quiet = TRUE, set_units = FALSE,
lower = enviro_par$T_air - 30,
upper = enviro_par$T_air + 30
)
if (is.null(fit$result)) {
fit <- list(root = NA, f.root = NA, convergence = 1)
} else {
fit <- fit$result
}
soln <- data.frame(T_leaf = fit$root, value = fit$f.root,
convergence = dplyr::if_else(is.null(fit$convergence), 0, 1))
if (!quiet) {
" done" %>%
crayon::green() %>%
message()
}
soln
}
#' Evaporation (mol / (m^2 s))
#'
#' @inheritParams .get_Rabs
#' @inheritParams .get_Sr
#' @inheritParams tleaves
#'
#' @return
#' \code{unitless = TRUE}: A value in units of mol / (m ^ 2 / s) number of class \code{numeric}
#' \code{unitless = FALSE}: A value in units of mol / (m ^ 2 / s) of class \code{units}
#'
#' @details
#' The leaf evaporation rate is the product of the total conductance to water vapour (m / s) and the water vapour gradient (mol / m^3):
#'
#' \deqn{E = g_\mathrm{tw} D_\mathrm{wv}}{E = g_tw D_wv}
#'
#' If \code{unitless = TRUE}, \code{T_leaf} is assumed in degrees K without checking.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' E(T_leaf, c(cs, ep, lp), FALSE)
#'
#' @export
#'
E <- function(T_leaf, pars, unitless) {
if (unitless & is(T_leaf, "units")) T_leaf %<>% drop_units()
if (!unitless) T_leaf %<>% set_units(K)
E <- .get_gtw(T_leaf, pars, unitless) * .get_dwv(T_leaf, pars, unitless)
if (!unitless) E %<>% set_units(mol / m ^ 2 / s)
E
}
Try the tealeaves package in your browser
Any scripts or data that you put into this service are public.
tealeaves documentation built on July 20, 2022, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions tleaves
Documented in tleaves
#' \code{tleaves}: find leaf temperatures for multiple parameter sets
#'
#' @param leaf_par A list of leaf parameters. This can be generated using the \code{make_leafpar} function.
#'
#' @param enviro_par A list of environmental parameters. This can be generated using the \code{make_enviropar} function.
#'
#' @param constants A list of physical constants. This can be generated using the \code{make_constants} function.
#'
#' @param progress Logical. Should a progress bar be displayed?
#'
#' @param quiet Logical. Should messages be displayed?
#'
#' @param set_units Logical. Should \code{units} be set? The function is faster when FALSE, but input must be in correct units or else results will be incorrect without any warning.
#'
#' @param parallel Logical. Should parallel processing be used via \code{\link[furrr]{future_map}}?
#'
#' @return
#'
#' \code{tleaves}: \cr
#' \cr
#' A tibble with the following \code{units} columns \cr
#'
#' \tabular{ll}{
#' \bold{Input:} \tab \cr
#' \code{abs_l} \tab Absorbtivity of longwave radiation (unitless) \cr
#' \code{abs_s} \tab Absorbtivity of shortwave radiation (unitless) \cr
#' \code{g_sw} \tab Stomatal conductance to H2O (\eqn{\mu}mol H2O / (m^2 s Pa)) \cr
#' \code{g_uw} \tab Cuticular conductance to H2O (\eqn{\mu}mol H2O / (m^2 s Pa)) \cr
#' \code{leafsize} Leaf characteristic dimension \tab (m) \cr
#' \code{logit_sr} \tab Stomatal ratio (logit transformed; unitless) \cr
#' \code{P} \tab Atmospheric pressure (kPa) \cr
#' \code{RH} \tab Relative humidity (unitless) \cr
#' \code{S_lw} \tab incident long-wave radiation flux density (W / m^2) \cr
#' \code{S_sw} \tab incident short-wave (solar) radiation flux density (W / m^2) \cr
#' \code{T_air} \tab Air temperature (K) \cr
#' \code{wind} \tab Wind speed (m / s) \cr
#' \bold{Output:} \tab \cr
#' \code{T_leaf} \tab Equilibrium leaf tempearture (K) \cr
#' \code{value} \tab Leaf energy balance (W / m^2) at \code{tleaf} \cr
#' \code{convergence} \tab Convergence code (0 = converged) \cr
#' \code{R_abs} \tab Total absorbed radiation (W / m^2; see \code{\link{.get_Rabs}}) \cr
#' \code{S_r} \tab Thermal infrared radiation loss (W / m^2; see \code{\link{.get_Sr}}) \cr
#' \code{H} \tab Sensible heat flux density (W / m^2; see \code{\link{.get_H}}) \cr
#' \code{L} \tab Latent heat flux density (W / m^2; see \code{\link{.get_L}}) \cr
#' \code{E} \tab Evapotranspiration (mol H2O/ (m^2 s))
#' }
#'
#' \code{tleaf}: \cr
#' \cr
#' A data.frame with the following numeric columns: \cr
#'
#' \tabular{ll}{
#' \code{T_leaf} \tab Equilibrium leaf temperature (K) \cr
#' \code{value} \tab Leaf energy balance (W / m^2) at \code{tleaf} \cr
#' \code{convergence} \tab Convergence code (0 = converged) \cr
#' \code{R_abs} \tab Total absorbed radiation (W / m^2; see \code{\link{.get_Rabs}}) \cr
#' \code{S_r} \tab Longwave re-radiation (W / m^2; see \code{\link{.get_Sr}}) \cr
#' \code{H} \tab Sensible heat flux density (W / m^2; see \code{\link{.get_H}}) \cr
#' \code{L}\tab Latent heat flux density (W / m^2; see \code{\link{.get_L}}) \cr
#' \code{E} \tab Evapotranspiration (mol H2O/ (m^2 s))
#' }
#'
#' @examples
#' # tleaf for single parameter set:
#'
#' leaf_par <- make_leafpar()
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' tleaf(leaf_par, enviro_par, constants)
#'
#' # tleaves for multiple parameter set:
#'
#' enviro_par <- make_enviropar(
#' replace = list(
#' T_air = set_units(c(293.15, 298.15), K)
#' )
#' )
#' tleaves(leaf_par, enviro_par, constants)
#'
#' @export
#'
tleaves <- function(leaf_par, enviro_par, constants, progress = TRUE,
quiet = FALSE, set_units = TRUE, parallel = FALSE) {
if (set_units) {
leaf_par %<>% tealeaves::leaf_par()
enviro_par %<>% tealeaves::enviro_par()
constants %<>% tealeaves::constants()
} else {
if (!quiet) warning("tleaves: units have not been checked prior to solving")
}
pars <- c(leaf_par, enviro_par)
tsky_function <- NULL
if (is.function(pars$T_sky)) {
tsky_function <- pars$T_sky
pars$T_sky <- NULL
}
par_units <- purrr::map(pars, units) %>%
magrittr::set_names(names(pars))
pars %<>%
purrr::map_if(~ inherits(.x, "units"), drop_units) %>%
make_parameter_sets()
if (!is.null(tsky_function)) {
# For this to work, tsky_function has to work with arg pars
# check with function like check_parameter_function(.f)
# also, this seems pretty inefficient, but I don't have a better solution
tidyr::crossing(
par = names(pars[[1]]),
i = 1:length(pars)
) %>%
purrr::pwalk(~ {
glue::glue("units(pars[[{i}]]${par}) <<- par_units${par}",
par = ..1, i = ..2) %>%
parse(text = .) %>%
eval()
})
pars <- pars %>%
purrr::map(~ {
.x$T_sky <- tsky_function(.x)
.x
})
par_units$T_sky <- units(pars[[1]]$T_sky)
pars <- pars %>%
purrr::map(~ {purrr::map_if(.x, ~ inherits(.x, "units"), drop_units)})
}
soln <- find_tleaves(pars, constants, progress, quiet, parallel)
pars %<>% purrr::map_dfr(purrr::flatten_dfr)
colnames(pars) %>%
glue::glue("units(pars${x}) <<- par_units${x}", x = .) %>%
parse(text = .) %>%
eval()
pars %<>% dplyr::bind_cols(soln)
pars$T_leaf %<>% set_units(K)
pars$E %<>% set_units(mol / m ^ 2 / s)
pars$R_abs %<>% set_units(W / m ^ 2)
pars$S_r %<>% set_units(W / m ^ 2)
pars$H %<>% set_units(W / m ^ 2)
pars$L %<>% set_units(W / m ^ 2)
pars
}
#' \code{tleaf}: find leaf temperatures for a single parameter set
#' @rdname tleaves
#' @export
tleaf <- function(leaf_par, enviro_par, constants, quiet = FALSE,
set_units = TRUE) {
if (set_units) {
leaf_par %<>% tealeaves::leaf_par()
enviro_par %<>% tealeaves::enviro_par()
constants %<>% tealeaves::constants()
} else {
if (!quiet) warning("tleaf: units have not been checked prior to solving")
}
if (is.function(enviro_par$T_sky)) {
enviro_par$T_sky <- enviro_par$T_sky(enviro_par)
}
# Drop units for solving
ulp <- purrr::map_if(leaf_par, ~ inherits(.x, "units"), drop_units)
uep <- purrr::map_if(enviro_par, ~ inherits(.x, "units"), drop_units)
ucs <- purrr::map_if(constants, ~ inherits(.x, "units"), drop_units)
soln <- find_tleaf(ulp, uep, ucs, quiet)
# Check results -----
if (soln$convergence == 1) {
"stats::uniroot did not converge, NA returned. Inspect parameters carefully." %>%
crayon::red() %>%
message()
}
# Energy balance components -----
components <- suppressWarnings(
if (is.na(soln$T_leaf)) {
list(R_abs = NA, S_r = NA, H = NA, L = NA, E = NA)
} else {
soln$T_leaf %>%
energy_balance(leaf_par = ulp, enviro_par = uep,
constants = ucs, quiet = TRUE, components = TRUE,
set_units = FALSE) %>%
magrittr::use_series("components")
}
)
soln %<>% dplyr::bind_cols(components)
if (is.na(soln$T_leaf)) {
soln %<>% dplyr::bind_cols(data.frame(Ar = NA, Gr = NA, Re = NA))
} else {
soln %<>% dplyr::bind_cols(Ar(soln$T_leaf, c(ulp, uep, ucs),
unitless = TRUE))
}
soln$T_leaf %<>% set_units(K)
# Adding g_bw to output
if (is.na(soln$T_leaf)) {
soln$g_bw <- NA
} else {
soln$g_bw <- .get_gbw(
T_leaf = drop_units(soln$T_leaf),
surface = "lower",
pars = c(ulp, uep, ucs),
unitless = TRUE
) +
.get_gbw(
T_leaf = drop_units(soln$T_leaf),
surface = "upper",
pars = c(ulp, uep, ucs),
unitless = TRUE
)
}
# Return -----
soln
}
#' Calculate leaf energy balance
#'
#' @inheritParams tleaves
#'
#' @param tleaf Leaf temperature in Kelvin. If input is numeric, it will be automatically converted to \code{units}.
#'
#' @param quiet Logical. Should a message appear about conversion from \code{numeric} to \code{units}? Useful for finding leaf temperature that balances heat transfer using \code{\link[stats]{uniroot}}.
#'
#' @param components Logical. Should leaf energy components be returned? Transpiration (in mol / (m^2 s)) also returned.
#'
#' @return A numeric value in W / m^2. Optionally, a named list of energy balance components in W / m^2 and transpiration in mol / (m^2 s).
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#' ep$T_sky <- ep$T_sky(ep)
#'
#' T_leaf <- set_units(298.15, K)
#'
#' energy_balance(T_leaf, lp, ep, cs, FALSE, TRUE, TRUE)
#'
#' @export
#'
energy_balance <- function(tleaf, leaf_par, enviro_par, constants,
quiet = FALSE, components = FALSE, set_units = FALSE) {
# Checks -----
stopifnot(length(set_units) == 1L & is.logical(set_units))
if (set_units) {
tleaf %<>%
set_units(K) %>%
drop_units()
leaf_par %<>% leaf_par()
enviro_par %<>% enviro_par()
constants %<>% constants()
}
stopifnot(length(quiet) == 1L & is.logical(quiet))
stopifnot(length(components) == 1L & is.logical(components))
pars <- c(leaf_par, enviro_par, constants) %>%
purrr::map_if(~ inherits(.x, "units"), drop_units)
# R_abs: total absorbed radiation (W m^-2) -----
R_abs <- .get_Rabs(pars, unitless = TRUE)
# S_r: longwave re-radiation (W m^-2) -----
S_r <- .get_Sr(tleaf, pars)
# H: sensible heat flux density (W m^-2) -----
H <- .get_H(tleaf, pars, unitless = TRUE)
# L: latent heat flux density (W m^-2) -----
L <- .get_L(tleaf, pars, unitless = TRUE)
# Return -----
if (components) {
# E: transpiration (mol / (m^2 s))
E <- .get_gtw(tleaf, pars, unitless = TRUE) *
.get_dwv(tleaf, pars, unitless = TRUE)
ret <- list(
energy_balance = R_abs - (S_r + H + L),
components = list(R_abs = R_abs, S_r = S_r, H = H, L = L, E = E)
)
ret$energy_balance %<>% set_units(W / m ^ 2)
ret$components$R_abs %<>% set_units(W / m ^ 2)
ret$components$S_r %<>% set_units(W / m ^ 2)
ret$components$H %<>% set_units(W / m ^ 2)
ret$components$L %<>% set_units(W / m ^ 2)
ret$components$E %<>% set_units(mol / m ^ 2 / s)
return(ret)
}
R_abs - (S_r + H + L)
}
#' R_abs: total absorbed radiation (W / m^2)
#'
#' @param pars Concatenated parameters (\code{leaf_par}, \code{enviro_par}, and \code{constants})
#'
#' @param unitless Logical. Should function use parameters with \code{units}? The function is faster when FALSE, but input must be in correct units or else results will be incorrect without any warning.
#'
#' @inheritParams tleaves
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' The following treatment follows Okajima et al. (2012):
#'
#' \deqn{R_\mathrm{abs} = \alpha_\mathrm{s} (1 + r) S_\mathrm{sw} + \alpha_\mathrm{l} \sigma (T_\mathrm{sky} ^ 4 + T_\mathrm{air} ^ 4)}{R_abs = \alpha_s (1 + r) S_sw + \alpha_l \sigma (T_sky ^ 4 + T_air ^ 4)}
#'
#' The incident longwave (aka thermal infrared) radiation is modeled from sky and air temperature \eqn{\sigma (T_\mathrm{sky} ^ 4 + T_\mathrm{air} ^ 4)}{\sigma (T_sky ^ 4 + T_air ^ 4)} where \eqn{T_\mathrm{sky}}{T_sky} is function of the air temperature and incoming solar shortwave radiation:
#'
#' \deqn{T_\mathrm{sky} = T_\mathrm{air} - 20 S_\mathrm{sw} / 1000}{T_sky = T_air - 20 S_sw / 1000}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\alpha_\mathrm{s}}{\alpha_s} \tab \code{abs_s} \tab absorbtivity of shortwave radiation (0.3 - 4 \eqn{\mu}m) \tab none \tab 0.80\cr
#' \eqn{\alpha_\mathrm{l}}{\alpha_l} \tab \code{abs_l} \tab absorbtivity of longwave radiation (4 - 80 \eqn{\mu}m) \tab none \tab 0.97\cr
#' \eqn{r} \tab \code{r} \tab reflectance for shortwave irradiance (albedo) \tab none \tab 0.2 \cr
#' \eqn{\sigma} \tab \code{s} \tab Stefan-Boltzmann constant \tab W / (m\eqn{^2} K\eqn{^4}) \tab 5.67e-08 \cr
#' \eqn{S_\mathrm{sw}}{S_sw} \tab \code{S_sw} \tab incident short-wave (solar) radiation flux density \tab W / m\eqn{^2} \tab 1000 \cr
#' \eqn{S_\mathrm{lw}}{S_lw} \tab \code{S_lw} \tab incident long-wave radiation flux density \tab W / m\eqn{^2} \tab calculated \cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15 \cr
#' \eqn{T_\mathrm{sky}}{T_sky} \tab \code{T_sky} \tab sky temperature \tab K \tab calculated
#' }
#'
#' @references
#'
#' Okajima Y, H Taneda, K Noguchi, I Terashima. 2012. Optimum leaf size predicted by a novel leaf energy balance model incorporating dependencies of photosynthesis on light and temperature. Ecological Research 27: 333-46.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#' ep$T_sky <- ep$T_sky(ep)
#'
#' tealeaves:::.get_Rabs(c(cs, ep, lp), FALSE)
#'
.get_Rabs <- function(pars, unitless) {
if (unitless) {
R_abs <- pars$abs_s * (1 + pars$r) * pars$S_sw + pars$abs_l * pars$s * (pars$T_sky ^ 4 + pars$T_air ^ 4)
} else {
R_abs <- pars$abs_s * (set_units(1) + pars$r) * pars$S_sw + pars$abs_l * pars$s * (pars$T_sky ^ 4 + pars$T_air ^ 4)
R_abs %<>% set_units(W / m ^ 2)
}
R_abs
}
#' S_r: longwave re-radiation (W / m^2)
#'
#' @inheritParams .get_Rabs
#' @param T_leaf Leaf temperature in Kelvin
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' \deqn{S_\mathrm{r} = 2 \sigma \alpha_\mathrm{l} T_\mathrm{air} ^ 4}{S_r = 2 \sigma \alpha_l T_air ^ 4 }
#'
#' The factor of 2 accounts for re-radiation from both leaf surfaces (Foster and Smith 1986). \cr
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\alpha_\mathrm{l}}{\alpha_l} \tab \code{abs_l} \tab absorbtivity of longwave radiation (4 - 80 \eqn{\mu}m) \tab none \tab 0.97\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{\sigma} \tab \code{s} \tab Stefan-Boltzmann constant \tab W / (m\eqn{^2} K\eqn{^4}) \tab 5.67e-08
#' }
#'
#' Note that leaf absorbtivity is the same value as leaf emissivity
#'
#' @references
#'
#' Foster JR, Smith WK. 1986. Influence of stomatal distribution on transpiration in low-wind environments. Plant, Cell & Environment 9: 751-9.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_Sr(T_leaf, c(cs, ep, lp))
#'
.get_Sr <- function(T_leaf, pars) 2 * pars$s * pars$abs_l * T_leaf ^ 4
#' H: sensible heat flux density (W / m^2)
#'
#' @inheritParams tleaves
#' @inheritParams .get_Rabs
#' @inheritParams .get_Sr
#'
#' @return Value in W / m\eqn{^2} of class \code{units}
#'
#' @details
#'
#' \deqn{H = P_\mathrm{a} c_p g_\mathrm{h} (T_\mathrm{leaf} - T_\mathrm{air})}{H = P_a c_p g_h * (T_leaf - T_air)}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{c_p} \tab \code{c_p} \tab heat capacity of air \tab J / (g K) \tab 1.01\cr
#' \eqn{g_\mathrm{h}}{g_h} \tab \code{g_h} \tab boundary layer conductance to heat \tab m / s \tab \link[=.get_gh]{calculated}\cr
#' \eqn{P_\mathrm{a}}{P_a} \tab \code{P_a} \tab density of dry air \tab g / m^3 \tab \link[=.get_Pa]{calculated}\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @seealso \code{\link{.get_gh}}, \code{\link{.get_Pa}}
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_H(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_H <- function(T_leaf, pars, unitless) {
# Density of dry air
P_a <- .get_Pa(T_leaf, pars, unitless)
# Boundary layer conductance to heat
g_h <- sum(.get_gh(T_leaf, "lower", pars, unitless),
.get_gh(T_leaf, "upper", pars, unitless))
H <- P_a * pars$c_p * g_h * (T_leaf - pars$T_air)
H
}
#' P_a: density of dry air (g / m^3)
#'
#' @inheritParams .get_H
#'
#' @return Value in g / m\eqn{^3} of class \code{units}
#'
#' @details
#'
#' \deqn{P_\mathrm{a} = P / (R_\mathrm{air} (T_\mathrm{leaf} - T_\mathrm{air}) / 2)}{P_a = P / (R_air (T_leaf - T_air) / 2)}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{P} \tab \code{P} \tab atmospheric pressure \tab kPa \tab 101.3246\cr
#' \eqn{R_\mathrm{air}}{R_air} \tab \code{R_air} \tab specific gas constant for dry air \tab J / (kg K) \tab 287.058\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_Pa(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_Pa <- function(T_leaf, pars, unitless) {
P_a <- pars$P / (pars$R_air * (pars$T_air + T_leaf) / 2)
if (unitless) {
P_a %<>% magrittr::multiply_by(1e6)
} else {
P_a %<>% set_units(g / m ^ 3)
}
P_a
}
#' g_h: boundary layer conductance to heat (m / s)
#'
#' @inheritParams .get_H
#' @param surface Leaf surface (lower or upper)
#'
#' @return Value in m/s of class \code{units}
#'
#' @details
#'
#' \deqn{g_\mathrm{h} = D_\mathrm{h} Nu / d}{g_h = D_h Nu / d}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{h}}{D_h} \tab \code{D_h} \tab diffusion coefficient for heat in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{Nu} \tab \code{Nu} \tab Nusselt number \tab none \tab \link[=.get_nu]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gh(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_gh <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
# Calculate diffusion coefficient to heat
D_h <- .get_Dx(pars$D_h0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
# Calculate Nusselt numbers
Nu <- .get_nu(T_leaf, surface, pars, unitless)
D_h * Nu / pars$leafsize
}
#' D_x: Calculate diffusion coefficient for a given temperature and pressure
#'
#' @param D_0 Diffusion coefficient at 273.15 K (0 °C) and 101.3246 kPa
#' @param Temp Temperature in Kelvin
#' @param eT Exponent for temperature dependence of diffusion
#' @param P Atmospheric pressure in kPa
#' @inheritParams .get_Rabs
#' @inheritParams tleaves
#'
#' @return Value in m\eqn{^2}/s of class \code{units}
#'
#' @details
#'
#' \deqn{D = D_\mathrm{0} (T / 273.15) ^ {eT} (101.3246 / P)}{D = D_0 [(T / 273.15) ^ eT] (101.3246 / P)}
#' \cr
#' According to Montieth & Unger (2013), eT is generally between 1.5 and 2. Their data in Appendix 3 indicate \eqn{eT = 1.75} is reasonable for environmental physics.
#'
#' @references
#'
#' Monteith JL, Unsworth MH. 2013. Principles of Environmental Physics. 4th edition. Academic Press, London.
#'
#' @examples
#'
#' tealeaves:::.get_Dx(
#' D_0 = set_units(2.12e-05, m^2/s),
#' Temp = set_units(298.15, K),
#' eT = set_units(1.75),
#' P = set_units(101.3246, kPa),
#' unitless = FALSE
#' )
#'
.get_Dx <- function(D_0, Temp, eT, P, unitless) {
# See Eq. 3.10 in Monteith & Unger ed. 4
if (unitless) {
Dx <- D_0 * (Temp / 273.15) ^ eT * (101.3246 / P)
} else {
Dx <- D_0 *
drop_units((set_units(Temp, K) / set_units(273.15, K))) ^ drop_units(eT) *
drop_units((set_units(101.3246, kPa) / set_units(P, kPa)))
}
Dx
}
#' Gr: Grashof number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' \deqn{Gr = t_\mathrm{air} G d ^ 3 |T_\mathrm{v,leaf} - T_\mathrm{v,air}| / D_\mathrm{m} ^ 2}{Gr = t_air G d ^ 3 abs(Tv_leaf - Tv_air) / D_m ^ 2}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{m}}{D_m} \tab \code{D_m} \tab diffusion coefficient of momentum in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{G} \tab \code{G} \tab gravitational acceleration \tab m / s\eqn{^2} \tab 9.8\cr
#' \eqn{t_\mathrm{air}}{t_air} \tab \code{t_air} \tab coefficient of thermal expansion of air \tab 1 / K \tab 1 / Temp \cr
#' \eqn{T_\mathrm{v,air}}{Tv_air} \tab \code{Tv_air} \tab virtual air temperature \tab K \tab \link[=.get_Tv]{calculated}\cr
#' \eqn{T_\mathrm{v,leaf}}{Tv_leaf} \tab \code{Tv_leaf} \tab virtual leaf temperature \tab K \tab \link[=.get_Tv]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gr(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_gr <- function(T_leaf, pars, unitless) {
# Calculate virtual temperature
# Assumes inside of leaf is 100% RH
Tv_leaf <- .get_Tv(T_leaf, .get_ps(T_leaf, pars$P, unitless), pars$P,
pars$epsilon, unitless)
Tv_air <- .get_Tv(pars$T_air, pars$RH * .get_ps(pars$T_air, pars$P, unitless), pars$P,
pars$epsilon, unitless)
D_m <- .get_Dx(pars$D_m0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
if (unitless) {
Gr <- (1 / pars$T_air) * pars$G * pars$leafsize ^ 3 *
abs(Tv_leaf - Tv_air) / D_m ^ 2
} else {
Gr <- (set_units(1) / pars$T_air) * pars$G * pars$leafsize ^ 3 *
abs(Tv_leaf - Tv_air) / D_m ^ 2
}
Gr
}
#' Calculate virtual temperature
#'
#' @inheritParams .get_Dx
#' @param p water vapour pressure in kPa
#' @param epsilon ratio of water to air molar masses (unitless)
#'
#' @return Value in K of class \code{units}
#'
#' @details
#'
#' \deqn{T_\mathrm{v} = T / [1 - (1 - \epsilon) (p / P)]}{T_v = T / [1 - (1 - epsilon) (p / P)]}
#'
#' Eq. 2.35 in Monteith & Unsworth (2013) \cr
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{\epsilon} \tab \code{epsilon} \tab ratio of water to air molar masses \tab unitless \tab 0.622 \cr
#' \eqn{p} \tab \code{p} \tab water vapour pressure \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{P} \tab \code{P} \tab atmospheric pressure \tab kPa \tab 101.3246
#' }
#'
#' @references
#'
#' Monteith JL, Unsworth MH. 2013. Principles of Environmental Physics. 4th edition. Academic Press, London.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#' p <- ep$RH * tealeaves:::.get_ps(T_leaf, ep$P, FALSE)
#' tealeaves:::.get_Tv(T_leaf, p, ep$P, cs$epsilon, FALSE)
#'
.get_Tv <- function(Temp, p, P, epsilon, unitless) {
if (unitless) {
Tv <- Temp / (1 - (1 - epsilon) * (p / P))
} else {
Tv <- set_units(Temp, K) /
(set_units(1) - (set_units(1) - epsilon) *
(set_units(p, kPa) / set_units(P, kPa)))
}
Tv
}
#' Saturation water vapour pressure (kPa)
#'
#' @inheritParams .get_Dx
#'
#' @return Value in kPa of class \code{units}
#'
#' @details
#'
#' Goff-Gratch equation (see http://cires1.colorado.edu/~voemel/vp.html) \cr
#' \cr
#' This equation assumes P = 1 atm = 101.3246 kPa, otherwise boiling temperature needs to change \cr
#'
#' @references \url{http://cires1.colorado.edu/~voemel/vp.html}
#'
#' @examples
#'
#' T_leaf <- set_units(298.15, K)
#' P <- set_units(101.3246, kPa)
#' tealeaves:::.get_ps(T_leaf, P, FALSE)
#'
.get_ps <- function(Temp, P, unitless) {
# Goff-Gratch equation (see http://cires1.colorado.edu/~voemel/vp.html)
# This assumes P = 1 atm = 101.3246 kPa, otherwise boiling temperature needs to change
# This returns p_s in hPa
if (unitless) {
P %<>% magrittr::multiply_by(10)
} else {
Temp %<>% set_units(K) %>% drop_units()
P %<>% set_units(hPa) %>% drop_units()
}
p_s <- 10 ^ (-7.90298 * (373.16 / Temp - 1) +
5.02808 * log10(373.16 / Temp) -
1.3816e-7 * (10 ^ (11.344 * (1 - Temp / 373.16) - 1)) +
8.1328e-3 * (10 ^ (-3.49149 * (373.16 / Temp - 1)) - 1) +
log10(P))
# Convert to kPa
if (unitless) {
p_s %<>% magrittr::multiply_by(0.1)
} else {
p_s %<>%
set_units(hPa) %>%
set_units(kPa)
}
p_s
}
#' Re: Reynolds number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' \deqn{Re = u d / D_\mathrm{m}}{Re = u d / D_m}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{m}}{D_m} \tab \code{D_m} \tab diffusion coefficient of momentum in air \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{u} \tab \code{wind} \tab windspeed \tab m / s \tab 2
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_re(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_re <- function(T_leaf, pars, unitless) {
D_m <- .get_Dx(pars$D_m0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
Re <- pars$wind * pars$leafsize / D_m
Re
}
#' Nu: Nusselt number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' The Nusselt number depends on a combination how much free or forced convection predominates. For mixed convection: \cr
#' \cr
#' \deqn{Nu = (a Re ^ b) ^ {3.5} + (c Gr ^ d) ^ {3.5}) ^ {1 / 3.5}}{Nu = (a Re ^ b) ^ 3.5 + (c Gr ^ d) ^ 3.5) ^ (1 / 3.5)}
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{a, b, c, d} \tab \code{a, b, c, d} \tab empirical coefficients \tab none \tab \link[=make_constants]{calculated}\cr
#' \eqn{Gr} \tab \code{Gr} \tab Grashof number \tab none \tab \link[=.get_gr]{calculated}\cr
#' \eqn{Re} \tab \code{Re} \tab Reynolds number \tab none \tab \link[=.get_re]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_nu(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_nu <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
Gr <- .get_gr(T_leaf, pars, unitless)
Re <- .get_re(T_leaf, pars, unitless)
# Archemides number
# Ar <- Gr / Re ^ 2
cons <- pars$nu_constant(Re, "forced", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_forced <- cons$a * Re ^ cons$b
} else {
Nu_forced <- cons$a * drop_units(Re) ^ cons$b
}
cons <- pars$nu_constant(Re, "free", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_free <- cons$a * Gr ^ cons$b
} else {
Nu_free <- cons$a * drop_units(Gr) ^ cons$b
}
Nu <- (Nu_forced ^ 3.5 + Nu_free ^ 3.5) ^ (1 / 3.5)
if (!unitless) Nu %<>% set_units()
Nu
}
#' L: Latent heat flux density (W / m^2)
#'
#' @inheritParams .get_H
#'
#' @return Value in W / m^2 of class \code{units}
#'
#' @details
#'
#' \deqn{L = h_\mathrm{vap} g_\mathrm{tw} d_\mathrm{wv}}{L = h_vap g_tw d_wv}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d_\mathrm{wv}}{d_wv} \tab \code{d_wv} \tab water vapour gradient \tab mol / m ^ 3 \tab \link[=.get_dwv]{calculated} \cr
#' \eqn{h_\mathrm{vap}}{h_vap} \tab \code{h_vap} \tab latent heat of vaporization \tab J / mol \tab \link[=.get_hvap]{calculated} \cr
#' \eqn{g_\mathrm{tw}}{g_tw} \tab \code{g_tw} \tab total conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab \link[=.get_gtw]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_L(T_leaf, c(cs, ep, lp), FALSE)
#'
.get_L <- function(T_leaf, pars, unitless) {
h_vap <- .get_hvap(T_leaf, unitless)
g_tw <- .get_gtw(T_leaf, pars, unitless)
d_wv <- .get_dwv(T_leaf, pars, unitless)
L <- h_vap * g_tw * d_wv
if (!unitless) L %<>% set_units(W / m ^ 2)
L
}
#' d_wv: water vapour gradient (mol / m ^ 3)
#'
#' @inheritParams .get_H
#'
#' @return Value in mol / m^3 of class \code{units}
#'
#' @details
#'
#' \bold{Water vapour gradient:} The water vapour pressure differential from inside to outside of the leaf is the saturation water vapor pressure inside the leaf (p_leaf) minus the water vapor pressure of the air (p_air):
#'
#' \deqn{d_\mathrm{wv} = p_\mathrm{leaf} / (R T_\mathrm{leaf}) - RH p_\mathrm{air} / (R T_\mathrm{air})}{d_wv = p_leaf / (R T_leaf) - RH p_air / (R T_air)}
#'
#' Note that water vapor pressure is converted from kPa to mol / m^3 using ideal gas law. \cr
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{p_\mathrm{air}}{p_air} \tab \code{p_air} \tab saturation water vapour pressure of air \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{p_\mathrm{leaf}}{p_leaf} \tab \code{p_leaf} \tab saturation water vapour pressure inside the leaf \tab kPa \tab \link[=.get_ps]{calculated}\cr
#' \eqn{R} \tab \code{R} \tab ideal gas constant \tab J / (mol K) \tab 8.3144598\cr
#' \eqn{\mathrm{RH}}{RH} \tab \code{RH} \tab relative humidity \tab \% \tab 0.50\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' # Water vapour gradient:
#'
#' leaf_par <- make_leafpar()
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' pars <- c(leaf_par, enviro_par, constants)
#' T_leaf <- set_units(300, K)
#' T_air <- set_units(298.15, K)
#' p_leaf <- set_units(35.31683, kPa)
#' p_air <- set_units(31.65367, kPa)
#'
#' d_wv <- p_leaf / (pars$R * T_leaf) - pars$RH * p_air / (pars$R * T_air)
#'
.get_dwv <- function(T_leaf, pars, unitless) {
# Water vapour differential converted from kPa to mol m ^ -3 using ideal gas law
d_wv <- .get_ps(T_leaf, pars$P, unitless) / (pars$R * T_leaf) -
pars$RH * .get_ps(pars$T_air, pars$P, unitless) / (pars$R * pars$T_air)
if (unitless) {
d_wv %<>% magrittr::multiply_by(1e3)
} else {
d_wv %<>% set_units(mol / m ^ 3)
}
d_wv
}
#' g_tw: total conductance to water vapour (m/s)
#'
#' @inheritParams .get_H
#'
#' @return Value in m/s of class \code{units}
#'
#' @details
#'
#' \bold{Total conductance to water vapor:} The total conductance to water vapor (\eqn{g_\mathrm{tw}}{g_tw}) is the sum of the parallel lower (abaxial) and upper (adaxial) conductances:
#'
#' \deqn{g_\mathrm{tw} = g_\mathrm{w,lower} + g_\mathrm{w,upper}}{g_tw = gw_lower + gw_upper}
#'
#' The conductance to water vapor on each surface is a function of parallel stomatal (\eqn{g_\mathrm{sw}}{g_sw}) and cuticular (\eqn{g_\mathrm{uw}}{g_uw}) conductances in series with the boundary layer conductance (\eqn{g_\mathrm{bw}}{g_bw}). The stomatal, cuticular, and boundary layer conductance on the lower surface are:
#'
#' \deqn{g_\mathrm{sw,lower} = g_\mathrm{sw} (1 - sr) R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{gsw_lower = g_sw (1 - sr) R (T_leaf + T_air) / 2}
#' \deqn{g_\mathrm{uw,lower} = g_\mathrm{uw} / 2 R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{guw_lower = g_uw / 2 R (T_leaf + T_air) / 2}
#' \cr
#' See \code{\link{.get_gbw}} for details on calculating boundary layer conductance. The equations for the upper surface are:
#'
#' \deqn{g_\mathrm{sw,upper} = g_\mathrm{sw} sr R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{gsw_upper = g_sw sr R (T_leaf + T_air) / 2}
#' \deqn{g_\mathrm{uw,upper} = g_\mathrm{uw} / 2 R (T_\mathrm{leaf} + T_\mathrm{air}) / 2}{guw_upper = g_uw / 2 R (T_leaf + T_air) / 2}
#' \cr
#' Note that the stomatal and cuticular conductances are given in units of (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) (see \code{\link{make_leafpar}}) and converted to m/s using the ideal gas law. The total leaf stomatal (\eqn{g_\mathrm{sw}}{g_sw}) and cuticular (\eqn{g_\mathrm{uw}}{g_uw}) conductances are partitioned across lower and upper surfaces. The stomatal conductance on each surface depends on stomatal ratio (sr); the cuticular conductance is assumed identical on both surfaces.
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{g_\mathrm{sw}}{g_sw} \tab \code{g_sw} \tab stomatal conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab 5\cr
#' \eqn{g_\mathrm{uw}}{g_uw} \tab \code{g_uw} \tab cuticular conductance to H2O \tab (\eqn{\mu}mol H2O) / (m\eqn{^2} s Pa) \tab 0.1\cr
#' \eqn{R} \tab \code{R} \tab ideal gas constant \tab J / (mol K) \tab 8.3144598\cr
#' \eqn{\mathrm{logit}(sr)}{logit(sr)} \tab \code{logit_sr} \tab stomatal ratio (logit transformed) \tab none \tab 0 = logit(0.5)\cr
#' \eqn{T_\mathrm{air}}{T_air} \tab \code{T_air} \tab air temperature \tab K \tab 298.15\cr
#' \eqn{T_\mathrm{leaf}}{T_leaf} \tab \code{T_leaf} \tab leaf temperature \tab K \tab input
#' }
#'
#' @examples
#'
#' # Total conductance to water vapor
#'
#' ## Hypostomatous leaf; default parameters
#' leaf_par <- make_leafpar(replace = list(logit_sr = set_units(-Inf)))
#' enviro_par <- make_enviropar()
#' constants <- make_constants()
#' pars <- c(leaf_par, enviro_par, constants)
#' T_leaf <- set_units(300, K)
#'
#' ## Fixing boundary layer conductance rather than calculating
#' gbw_lower <- set_units(0.1, m / s)
#' gbw_upper <- set_units(0.1, m / s)
#'
#' # Lower surface ----
#' ## Note that pars$logit_sr is logit-transformed! Use stats::plogis() to convert to proportion.
#' gsw_lower <- set_units(pars$g_sw * (set_units(1) - stats::plogis(pars$logit_sr)) * pars$R *
#' ((T_leaf + pars$T_air) / 2), "m / s")
#' guw_lower <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
#' gtw_lower <- 1 / (1 / (gsw_lower + guw_lower) + 1 / gbw_lower)
#'
#' # Upper surface ----
#' gsw_upper <- set_units(pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
#' ((T_leaf + pars$T_air) / 2), m / s)
#' guw_upper <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
#' gtw_upper <- 1 / (1 / (gsw_upper + guw_upper) + 1 / gbw_upper)
#'
#' ## Lower and upper surface are in parallel
#' g_tw <- gtw_lower + gtw_upper
#'
.get_gtw <- function(T_leaf, pars, unitless) {
# Lower surface ----
gbw_lower <- .get_gbw(T_leaf, "lower", pars, unitless)
# Convert stomatal and cuticular conductance from molar to 'engineering' units
# See email from Tom Buckley (July 4, 2017)
if (unitless) {
gsw_lower <- (pars$g_sw * (1 - stats::plogis(pars$logit_sr)) * pars$R *
((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
guw_lower <- (pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
} else {
gsw_lower <- set_units(pars$g_sw * (set_units(1) - stats::plogis(pars$logit_sr)) *
pars$R * ((T_leaf + pars$T_air) / 2), m / s)
guw_lower <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
}
gtw_lower <- 1 / (1 / (gsw_lower + guw_lower) + 1 / gbw_lower)
# Upper surface ----
gbw_upper <- .get_gbw(T_leaf, "upper", pars, unitless)
# Convert stomatal and cuticular conductance from molar to 'engineering' units
# See email from Tom Buckley (July 4, 2017)
if (unitless) {
gsw_upper <- (pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
guw_upper <- (pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2)) %>%
magrittr::divide_by(1e6)
} else {
gsw_upper <- set_units(pars$g_sw * stats::plogis(pars$logit_sr) * pars$R *
((T_leaf + pars$T_air) / 2), m / s)
guw_upper <- set_units(pars$g_uw * 0.5 * pars$R * ((T_leaf + pars$T_air) / 2), m / s)
}
gtw_upper <- 1 / (1 / (gsw_upper + guw_upper) + 1 / gbw_upper)
# Lower and upper surface are in parallel
g_tw <- gtw_lower + gtw_upper
g_tw
}
#' h_vap: heat of vaporization (J / mol)
#'
#' @inheritParams .get_H
#'
#' @return Value in J/mol of class \code{units}
#'
#' @details
#'
#' \bold{Heat of vaporization:} The heat of vaporization (\eqn{h_\mathrm{vap}}{h_vap}) is a function of temperature. I used data from on temperature and \eqn{h_\mathrm{vap}}{h_vap} from Nobel (2009, Appendix 1) to estimate a linear regression. See Examples.
#'
#' @examples
#'
#' # Heat of vaporization and temperature
#' ## data from Nobel (2009)
#'
#' T_K <- 273.15 + c(0, 10, 20, 25, 30, 40, 50, 60)
#' h_vap <- 1e3 * c(45.06, 44.63, 44.21, 44.00,
#' 43.78, 43.35, 42.91, 42.47) # (in J / mol)
#'
#' fit <- lm(h_vap ~ T_K)
#'
#' ## coefficients are 56847.68250 J / mol and 43.12514 J / (mol K)
#'
#' coef(fit)
#'
#' T_leaf <- 298.15
#' h_vap <- set_units(56847.68250, J / mol) -
#' set_units(43.12514, J / mol / K) * set_units(T_leaf, K)
#'
#' ## h_vap at 298.15 K is 43989.92 [J/mol]
#'
#' set_units(h_vap, J / mol)
#'
#' tealeaves:::.get_hvap(set_units(298.15, K), FALSE)
#'
#' @references Nobel PS. 2009. Physicochemical and Environmental Plant Physiology. 4th Edition. Academic Press.
#'
.get_hvap <- function(T_leaf, unitless) {
# Equation from Foster and Smith 1986 seems to be off:
# h_vap <- 4.504e4 - 41.94 * T_leaf
# Instead, using regression based on data from Nobel (2009, 4th Ed, Appendix 1)
# T_K <- 273.15 + c(0, 10, 20, 25, 30, 40, 50, 60)
# h_vap <- 1e3 * c(45.06, 44.63, 44.21, 44, 43.78, 43.35, 42.91, 42.47) # (in J / mol)
# fit <- lm(h_vap ~ T_K)
if (unitless) {
h_vap <- 56847.68250 - 43.12514 * T_leaf
} else {
h_vap <- set_units(56847.68250, J / mol) -
set_units(43.12514, J / mol / K) * set_units(T_leaf, K)
h_vap %<>% set_units(J / mol)
}
h_vap
}
#' g_bw: Boundary layer conductance to water vapour (m / s)
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return Value in m / s of class \code{units}
#'
#' @details
#'
#' \deqn{g_\mathrm{bw} = D_\mathrm{w} Sh / d}{g_bw = D_w Sh / d}
#'
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{d} \tab \code{leafsize} \tab Leaf characteristic dimension in meters \tab m \tab 0.1\cr
#' \eqn{D_\mathrm{w}}{D_w} \tab \code{D_w} \tab diffusion coefficient for water vapour \tab m\eqn{^2} / s \tab \link[=.get_Dx]{calculated}\cr
#' \eqn{Sh} \tab \code{Sh} \tab Sherwood number \tab none \tab \link[=.get_sh]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_gbw(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_gbw <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
D_w <- .get_Dx(pars$D_w0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
# Calculate Sherwood numbers
Sh <- .get_sh(T_leaf, surface, pars, unitless)
D_w * Sh / pars$leafsize
}
#' Sh: Sherwood number
#'
#' @inheritParams .get_H
#' @inheritParams .get_gh
#'
#' @return A unitless number of class \code{units}
#'
#' @details
#'
#' The Sherwood number depends on a combination how much free or forced convection predominates. For mixed convection: \cr
#' \cr
#' \deqn{Sh = (a Re ^ b) ^ {3.5} + (c Gr ^ d) ^ {3.5}) ^ {1 / 3.5}}{Sh = (a Re ^ b) ^ 3.5 + (c Gr ^ d) ^ 3.5) ^ (1 / 3.5)}
#' \cr
#' \tabular{lllll}{
#' \emph{Symbol} \tab \emph{R} \tab \emph{Description} \tab \emph{Units} \tab \emph{Default}\cr
#' \eqn{a, b, c, d} \tab \code{a, b, c, d} \tab empirical coefficients \tab none \tab \link[=make_constants]{calculated}\cr
#' \eqn{Gr} \tab \code{Gr} \tab Grashof number \tab none \tab \link[=.get_gr]{calculated}\cr
#' \eqn{Re} \tab \code{Re} \tab Reynolds number \tab none \tab \link[=.get_re]{calculated}
#' }
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' tealeaves:::.get_sh(T_leaf, "lower", c(cs, ep, lp), FALSE)
#'
.get_sh <- function(T_leaf, surface, pars, unitless) {
surface %<>% match.arg(c("lower", "upper"))
Gr <- .get_gr(T_leaf, pars, unitless)
Re <- .get_re(T_leaf, pars, unitless)
# Archemides number
# Ar <- Gr / Re ^ 2
D_h <- .get_Dx(pars$D_h0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
D_w <- .get_Dx(pars$D_w0, (pars$T_air + T_leaf) / 2, pars$eT, pars$P, unitless)
cons <- pars$nu_constant(Re, "forced", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_forced <- cons$a * Re ^ cons$b
Sh_forced <- Nu_forced * (D_h / D_w) ^ pars$sh_constant("forced")
} else {
Nu_forced <- cons$a * drop_units(Re) ^ cons$b
Sh_forced <- Nu_forced * drop_units(D_h / D_w) ^ pars$sh_constant("forced")
}
cons <- pars$nu_constant(Re, "free", pars$T_air, T_leaf, surface, unitless)
if (unitless) {
Nu_free <- cons$a * Gr ^ cons$b
Sh_free <- Nu_free * (D_h / D_w) ^ pars$sh_constant("free")
} else {
Nu_free <- cons$a * drop_units(Gr) ^ cons$b
Sh_free <- Nu_free * drop_units(D_h / D_w) ^ pars$sh_constant("free")
}
Sh <- (Sh_forced ^ 3.5 + Sh_free ^ 3.5) ^ (1 / 3.5)
if (!unitless) Sh %<>% set_units()
Sh
}
make_parameter_sets <- function(pars) {
pars %<>%
names() %>%
glue::glue("{x} = pars${x}", x = .) %>%
stringr::str_c(collapse = ", ") %>%
glue::glue("tidyr::crossing({x})", x = .) %>%
parse(text = .) %>%
eval() %>%
purrr::transpose()
pars
}
find_tleaves <- function(par_sets, constants, progress, quiet, parallel) {
# This function is not intended to be called directly by users.
# It assumes all parameters have correct units
if (!quiet) {
glue::glue("\nSolving for T_leaf from {n} parameter set{s}...",
n = length(par_sets),
s = dplyr::if_else(length(par_sets) > 1, "s", "")) %>%
crayon::green() %>%
message(appendLF = FALSE)
}
if (parallel) future::plan("multisession")
if (progress & !parallel) pb <- dplyr::progress_estimated(length(par_sets))
soln <- suppressWarnings(
par_sets %>%
furrr::future_map_dfr(~{
ret <- tleaf(.x, .x, constants, quiet = TRUE, set_units = FALSE)
if (progress & !parallel) pb$tick()$print()
ret
}, .progress = progress)
)
soln
}
find_tleaf <- function(leaf_par, enviro_par, constants, quiet) {
# This function is not intended to be called directly by users.
# It assumes all parameters have correct units
# Balance energy fluxes -----
if (!quiet) {
"\nSolving for T_leaf ..." %>%
crayon::green() %>%
message(appendLF = FALSE)
}
.f <- function(T_leaf, ...) {
eb <- energy_balance(T_leaf, ...)
eb
}
fit <- safely_uniroot(
f = .f,
leaf_par = leaf_par, enviro_par = enviro_par, constants = constants,
quiet = TRUE, set_units = FALSE,
lower = enviro_par$T_air - 30,
upper = enviro_par$T_air + 30
)
if (is.null(fit$result)) {
fit <- list(root = NA, f.root = NA, convergence = 1)
} else {
fit <- fit$result
}
soln <- data.frame(T_leaf = fit$root, value = fit$f.root,
convergence = dplyr::if_else(is.null(fit$convergence), 0, 1))
if (!quiet) {
" done" %>%
crayon::green() %>%
message()
}
soln
}
#' Evaporation (mol / (m^2 s))
#'
#' @inheritParams .get_Rabs
#' @inheritParams .get_Sr
#' @inheritParams tleaves
#'
#' @return
#' \code{unitless = TRUE}: A value in units of mol / (m ^ 2 / s) number of class \code{numeric}
#' \code{unitless = FALSE}: A value in units of mol / (m ^ 2 / s) of class \code{units}
#'
#' @details
#' The leaf evaporation rate is the product of the total conductance to water vapour (m / s) and the water vapour gradient (mol / m^3):
#'
#' \deqn{E = g_\mathrm{tw} D_\mathrm{wv}}{E = g_tw D_wv}
#'
#' If \code{unitless = TRUE}, \code{T_leaf} is assumed in degrees K without checking.
#'
#' @examples
#'
#' library(tealeaves)
#'
#' cs <- make_constants()
#' ep <- make_enviropar()
#' lp <- make_leafpar()
#'
#' T_leaf <- set_units(298.15, K)
#'
#' E(T_leaf, c(cs, ep, lp), FALSE)
#'
#' @export
#'
E <- function(T_leaf, pars, unitless) {
if (unitless & is(T_leaf, "units")) T_leaf %<>% drop_units()
if (!unitless) T_leaf %<>% set_units(K)
E <- .get_gtw(T_leaf, pars, unitless) * .get_dwv(T_leaf, pars, unitless)
if (!unitless) E %<>% set_units(mol / m ^ 2 / s)
E
}
Try the tealeaves package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.