evaporative-length: Evaporative Characteristic Length

Description

The functions lc and lt compute the characteristic length $L_\mathrm{c}$ of stage-I evaporation from a soil and its ``target'' (expected) value $L_\mathrm{t}$, used to constrain estimates of nonlinear parameters of the Van Genuchten-Mualem (VGM) model for water retention curves and hydraulic conductivity functions.

Usage

lc(alpha, n, tau, k0, e0, c0 = NULL, c3 = NULL, c4 = NULL)
lt(n, tau, e0, c0, c1, c2, c3, c4)

Value

A numeric scalar with the characteristic evaporative length (lc) or its expected value (lt).

Arguments

alpha

parameter $\alpha$ (inverse air entry pressure) of the VGM model, see wc_model and hc_model. For consistency with other quantities, the unit of $\alpha$ should be 1/meter [$\mathrm{m}^{-1}$].

n

parameter $n$ (shape parameter) of the VGM model, see wc_model and hc_model.

tau

parameter $\tau$ (tortuosity parameter) of the VGM model, see hc_model.

k0

saturated hydraulic conductivity $K_0$, see hc_model. If k0 is missing or equal to NA in calls of lc then k0 is approximated by the same relation as used for lt, see Details. For consistency with other quantities, the unit of $K_0$ should be meter/day [$\mathrm{m}\,\mathrm{d}^{-1}$].

e0

a numeric scalar with the stage-I rate of evaporation $E_0$ from a soil, see Details and soilhypfitIntro. For consistency with other quantities, the unit of $E_0$ should be meter/day [$\mathrm{m}\,\mathrm{d}^{-1}$].

c0, c1, c2, c3, c4

numeric constants to approximate the parameter $\alpha$ and the saturated hydraulic conductivity $K_0$ when computing $L_\mathrm{t}$, see Details and control_fit_wrc_hcc. For consistency with other quantities, the following units should be used for the constants:

c1: $\mathrm{m}^{-1}$,
c3: $\mathrm{m}\,\mathrm{d}^{-1}$.

The remaining constants are dimensionless.

Author

Andreas Papritz papritz@retired.ethz.ch.

Details

The characteristic length of stage-I evaporation $L_\mathrm{c}$ (Lehmann et al., 2008, 2020) is defined by $$ L_\mathrm{c}(\boldsymbol{\nu}, K_0, E_0) = \frac{ \frac{1}{\alpha \,§ n} \left(\frac{2n-1}{n-1}\right)^\frac{2n-1}{n} }{1 + \frac{E_0}{K_\mathrm{eff}}} $$ where $\boldsymbol{\nu}^\mathrm{T} = (\alpha, n, \tau)$ are the nonlinear parameters of the VGM model, $K_0$ is the saturated hydraulic conductivity, $E_0$ the stage-I evaporation rate and $K_\mathrm{eff} = 4\, K_\mathrm{VGM}(h_\mathrm{crit}; K_0, \boldsymbol{\nu})$ is the effective hydraulic conductivity at the critical pressure $$ h_\mathrm{crit} = \frac{1}{\alpha} \, \left(\frac{n-1}{n}\right)^\frac{1-2n}{n}, $$ see hc_model for the definition of $K_\mathrm{VGM}(h; K_0, \boldsymbol{\nu})$.

The quantity $L_\mathrm{t}$ is the expected value (``target'') of $L_\mathrm{c}$ for given shape ($n$) and tortuosity ($\tau$) parameters. To evaluate $L_\mathrm{t}$, the parameters $\alpha$ and $K_0$ are approximated by the following relations $$ \widehat{\alpha} = g_\alpha(n; c_0, c_1, c_2) = c_1 \, \frac{n - c_0}{1 + c_2 \, (n - c_0)}, $$ $$ \widehat{K}_0 = g_{K_0}(n; c_0, c_3, c_4) = c_3 \, (n - c_0)^{c_4}. $$ The default values for $c_0$ to $c_4$ (see argument approximation_alpha_k0 of control_fit_wrc_hcc) were estimated with data on African desert regions of the database ROSETTA3 (Zhang and Schaap, 2017), see Lehmann et al. (2020) for details.

References

Lehmann, P., Assouline, S., Or, D. (2008) Characteristic lengths affecting evaporative drying of porous media. Physical Review E, 77, 056309, tools:::Rd_expr_doi("10.1103/PhysRevE.77.056309").

Lehmann, P., Bickel, S., Wei, Z., Or, D. (2020) Physical Constraints for Improved Soil Hydraulic Parameter Estimation by Pedotransfer Functions. Water Resources Research 56, e2019WR025963, tools:::Rd_expr_doi("10.1029/2019WR025963").

Zhang, Y. , Schaap, M. G. 2017. Weighted recalibration of the Rosetta pedotransfer model with improved estimates of hydraulic parameter distributions and summary statistics (Rosetta3). Journal of Hydrology, 547, 39-53, tools:::Rd_expr_doi("10.1016/j.jhydrol.2017.01.004").

Examples

Run this code

# \donttest{
# use of \donttest{} because execution time exceeds 5 seconds

# estimate parameters of 4 samples of the Swiss forest soil dataset
# that have water retention (theta, all samples), saturated hydraulic conductivity
# (ksat) and optionally unsaturated hydraulic conductivity data
# (ku, samples "CH2_4" and "CH3_1")

data(swissforestsoils)

# select subset of data
sfs_subset <- droplevels(
  subset(
    swissforestsoils,
    layer_id %in% c("CH2_3", "CH2_4", "CH2_6", "CH3_1")
  ))

# extract ksat measurements
ksat <- sfs_subset[!duplicated(sfs_subset$layer_id), "ksat", drop = FALSE]
rownames(ksat) <- levels(sfs_subset$layer_id)
colnames(ksat) <- "k0"

# define number of cores for parallel computations
if(interactive()) ncpu <- parallel::detectCores() - 1L else ncpu <- 1L


# unconstrained estimation (global optimisation algorithm NLOPT_GN_MLSL)
# k0 fixed at measured ksat values
rsfs_uglob <- fit_wrc_hcc(
  wrc_formula = theta ~ head | layer_id,
  hcc_formula = ku ~ head | layer_id,
  data = sfs_subset,
  param = ksat,
  fit_param = default_fit_param(k0 = FALSE),
  control = control_fit_wrc_hcc(
    settings = "uglobal", pcmp = control_pcmp(ncores = ncpu)))
summary(rsfs_uglob)
coef(rsfs_uglob, lc = TRUE, gof = TRUE)


# constrained estimation by restricting ratio Lc/Lt to [0.5, 2]
# (global constrained optimisation algorithm NLOPT_GN_MLSL)
# k0 fixed at measured ksat values
rsfs_cglob <- update(
  rsfs_uglob,
  control = control_fit_wrc_hcc(
    settings = "cglobal", nloptr = control_nloptr(ranseed = 1),
    pcmp = control_pcmp(ncores = ncpu)))
summary(rsfs_cglob)
coef(rsfs_cglob, lc = TRUE, gof = TRUE)

# get initial parameter values from rsfs_cglob
ini_param <- cbind(
  coef(rsfs_cglob)[, c("alpha", "n")],
  ksat
)
# constrained estimation by restricting ratio Lc/Lt to [0.5, 2]
# (local constrained optimisation algorithm NLOPT_LD_CCSAQ)
# k0 fixed at measured ksat values
rsfs_cloc <- update(
  rsfs_uglob,
  param = ini_param,
  control = control_fit_wrc_hcc(
    settings = "clocal", nloptr = control_nloptr(ranseed = 1),
    pcmp = control_pcmp(ncores = ncpu)))
summary(rsfs_cloc)
coef(rsfs_cloc, lc = TRUE, gof = TRUE)

op <- par(mfrow = c(4, 2))
plot(rsfs_uglob, y = rsfs_cglob)
on.exit(par(op))

op <- par(mfrow = c(4, 2))
plot(rsfs_uglob, y = rsfs_cloc)
on.exit(par(op))

# }

Run the code above in your browser using DataLab