evaporative-length {soilhypfit}

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)

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.

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.

Value

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

Author(s)

Andreas Papritz [e-mail].

References

Lehmann, P., Assouline, S., Or, D. (2008) Characteristic lengths affecting evaporative drying of porous media. Physical Review E, 77, 056309, 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, 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, doi:10.1016/j.jhydrol.2017.01.004.

See Also

soilhypfitIntro for a description of the models and a brief summary of the parameter estimation approach;

fit_wrc_hcc for (constrained) estimation of parameters of models for soil water retention and hydraulic conductivity data;

control_fit_wrc_hcc for options to control fit_wrc_hcc;

soilhypfitmethods for common S3 methods for class fit_wrc_hcc;

vcov for computing (co-)variances of the estimated nonlinear parameters;

prfloglik_sample for profile loglikelihood computations;

wc_model and hc_model for currently implemented models for soil water retention curves and hydraulic conductivity functions;

Examples

# 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))