The soilhypfit Package

Description

This is a summary of the features and functionality of soilhypfit, a package in R for parameteric modelling of soil water retention and hydraulic conductivity data.

Details

Estimation approach

The main function, fit_wrc_hcc, estimates parameters of models for soil water retention and hydraulic conductivity by maximum likelihood (ml, default), maximum posterior density (mpd) estimation (Stewart et al., 1992) or nonlinear weighted least squares (wls) from data on volumetric soil water content, \boldsymbol{\theta}^\mathrm{T} =(\theta_1, \theta_2, ..., \theta_{n_\theta}), and/or soil hydraulic conductivity, \boldsymbol{K}^\mathrm{T} =(K_1, K_2, ..., K_{n_K}), both measured at given capillary pressure head, \boldsymbol{h}^\mathrm{T} =(h_1, h_2, ...).

For mpd and ml estimation, the models for the measurements are

\theta_i = \theta(h_i;\boldsymbol{\mu}, \boldsymbol{\nu}) + \varepsilon_{\theta,i}, \ \ i=1, 2, ..., n_\theta,

\log(K_j) = \log(K(h_j;\boldsymbol{\mu}, \boldsymbol{\nu})) + \varepsilon_{K,j}, \ \ j=1, 2, ..., n_K,

where \theta(h_i; \boldsymbol{\mu}, \boldsymbol{\nu}) and K(h_j; \boldsymbol{\mu}, \boldsymbol{\nu}) denote modelled water content and conductivity, \boldsymbol{\mu}^\mathrm{} and \boldsymbol{\nu}^\mathrm{} are the conditionally linear and nonlinear model parameters (see below and Bates and Watts, 1988, sec. 3.3.5), and \varepsilon_{\theta,i} and \varepsilon_{K,j} are independent, normally distributed errors with zero means and variances \sigma^2_\theta and \sigma^2_K, respectively.

Let

Q_{\theta}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{h}) =\left( \boldsymbol{\theta} - \boldsymbol{\theta}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} ) \right)^\mathrm{T} \boldsymbol{W}_\theta \left( \boldsymbol{\theta} - \boldsymbol{\theta}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} ) \right),

and

Q_{K}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{K}, \boldsymbol{h}) = \left( \log(\boldsymbol{K}) - \log(\boldsymbol{K}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} )) \right)^\mathrm{T} \boldsymbol{W}_K \left( \log(\boldsymbol{K}) - \log(\boldsymbol{K}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} )) \right),

denote the (possibly weighted) residual sums of squares between measurements and modelled values. \boldsymbol{\theta}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} ) ) and \boldsymbol{K}( \boldsymbol{h}; \boldsymbol{\mu}, \boldsymbol{\nu} ) are vectors with modelled values of water content and hydraulic conductivity, and \boldsymbol{W}_\theta and \boldsymbol{W}_K are optional diagonal matrices with weights w_{\theta,i} and w_{K,j}, respectively. The weights are products of case weights w^\prime_{\theta,i} and w^\prime_{K,j} and variable weights w_\theta, w_K, hence w_{\theta,i} = w_\theta\, w^\prime_{\theta,i} and w_{K,j} = w_K\, w^\prime_{K,j} .

The objective function for ml and mpd estimation is equal to (Stewart et al., 1992, eqs 14 and 15, respectively)

Q(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) = \frac{\kappa_\theta}{2} \log( Q_{\theta}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{h})) + \frac{\kappa_K}{2} \log( Q_{K}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{K}, \boldsymbol{h})),

where \kappa_v = n_v + 2 for mpd and \kappa_v = n_v for ml estimation, v \in (\theta, K), and weights w_{\theta,i} = w_{K,j} = 1. Note that Q(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h}) is equal to the negative logarithm of the concentrated loglikelihood or the concentrated posterior density, obtained by replacing \sigma^2_\theta and \sigma^2_K by their conditional maximum likelihood or maximum density estimates \widehat{\sigma}^2_\theta(\boldsymbol{\mu}), \boldsymbol{\nu}) and \widehat{\sigma}^2_K(\boldsymbol{\mu}) respectively (Stewart et al., 1992, p. 642).

For wls the objective function is equal to

Q(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) = Q_{\theta}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{h}) + Q_{K}(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{K}, \boldsymbol{h}).

If either water content or conductivity data are not available, then the respective terms are omitted from Q(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) .

The function fit_wrc_hcc does not attempt to estimate the parameters by minimising
Q(\boldsymbol{\mu}, \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h}) directly with respect to \boldsymbol{\mu}^\mathrm{} and \boldsymbol{\nu}^\mathrm{} . Rather, it exploits the fact that for given nonlinear parameters \boldsymbol{\nu}^\mathrm{} , the conditionally linear parameters \boldsymbol{\mu}^T = (\theta_r, \theta_s, \log(K_0)) can be estimated straightforwardly by minimising the conditional residual sums of squares

Q_\theta^*(\theta_r, \theta_s; \boldsymbol{\theta}, \boldsymbol{h}, \boldsymbol{\nu}) =\left( \boldsymbol{\theta} - [ \boldsymbol{1}, \boldsymbol{S}( \boldsymbol{h}; \boldsymbol{\nu} )] \, \left[ \begin{array}{c} \theta_r \\ \theta_s - \theta_r \end{array} \right] \right)^\mathrm{T} \boldsymbol{W}_\theta \left( \boldsymbol{\theta} - [ \boldsymbol{1}, \boldsymbol{S}( \boldsymbol{h}; \boldsymbol{\nu} )] \, \left[ \begin{array}{c} \theta_r \\ \theta_s - \theta_r \end{array} \right] \right)

with respect to \theta_r and \theta_s-\theta_r and/or

Q_K^*(K_0; \boldsymbol{K}, \boldsymbol{h}, \boldsymbol{\nu})= \left( \log(\boldsymbol{K}) - \log(K_0 \, \boldsymbol{k}( \boldsymbol{h}; \boldsymbol{\nu} )) \right)^\mathrm{T} \boldsymbol{W}_K \left( \log(\boldsymbol{K}) - \log(K_0 \, \boldsymbol{k}( \boldsymbol{h}; \boldsymbol{\nu} )) \right),

with respect to \log(K_0), where \boldsymbol{1}^\mathrm{} is a vector of ones, \boldsymbol{S}( \boldsymbol{h}; \boldsymbol{\nu} )^\mathrm{T} = ( S(h_1; \boldsymbol{\nu}), ..., S(h_{n_\theta}; \boldsymbol{\nu}) ) and \boldsymbol{k}( \boldsymbol{h}; \boldsymbol{\nu} )^\mathrm{T} = ( k(h_1; \boldsymbol{\nu}), ..., k(h_{n_K}; \boldsymbol{\nu}) ) are vectors of modelled water saturation and modelled relative conductivity values, \theta_r and \theta_s are the residual and saturated water content, and K_0 is the saturated hydraulic conductivity.

Unconstrained conditional estimates, say \widehat{\theta}_r(\boldsymbol{\nu}) , \widehat{\theta}_s(\boldsymbol{\nu}) - \widehat{\theta}_r(\boldsymbol{\nu}) and \widehat{\log(K_0)}(\boldsymbol{\nu}) can be easily obtained from the normal equations of the respective (weighted) least squares problems, and quadratic programming yields conditional (weighted) least squares estimates that honour the inequality constraints 0 \le \theta_r \le \theta_s \le 1.

Let \widehat{\boldsymbol{\mu}}(\boldsymbol{\nu})^\mathrm{T} = ( \widehat{\theta}_r(\boldsymbol{\nu}), \widehat{\theta}_s(\boldsymbol{\nu}), \widehat{\log(K_0)}(\boldsymbol{\nu}) ) be the conditional estimates of the linear parameters obtained by minimising Q_\theta^*(\theta_r, \theta_s; \boldsymbol{\theta}, \boldsymbol{h}, \boldsymbol{\nu}) , and Q_K^*(K_0; \boldsymbol{K}, \boldsymbol{h}, \boldsymbol{\nu}) , respectively. fit_wrc_hcc then estimates the nonlinear parameters by minimising Q( \widehat{\boldsymbol{\mu}}(\boldsymbol{\nu}), \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h}) with respect to \boldsymbol{\nu}^\mathrm{} by a nonlinear optimisation algorithm.

For mpd and ml estimation the variances of the model errors are estimated by (Stewart et al., 1992, eq. 16)

\widehat{\sigma}_\theta^2 = \frac{Q_{\theta}( \widehat{\boldsymbol{\mu}}(\widehat{\boldsymbol{\nu}}), \widehat{\boldsymbol{\nu}}; \boldsymbol{\theta}, \boldsymbol{h})}{\kappa_\theta},

and

\widehat{\sigma}_K^2 = \frac{Q_{K}( \widehat{\boldsymbol{\mu}}(\widehat{\boldsymbol{\nu}}), \widehat{\boldsymbol{\nu}}; \boldsymbol{K}, \boldsymbol{h})}{\kappa_K}.

Furthermore, for mpd and ml estimation, the covariance matrix of the estimated nonlinear parameters may be approximated by the inverse Hessian matrix of Q(\widehat{\boldsymbol{\mu}}(\boldsymbol{\nu}), \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h}) at the solution \widehat{\boldsymbol{\nu}} (Stewart and Sørensen, 1981), i.e.

\mathrm{Cov}[ \widehat{\boldsymbol{\nu}}, \widehat{\boldsymbol{\nu}}^T] \approx \boldsymbol{A}^{-1},

where

[\boldsymbol{A}]_{kl} = \frac{\partial^2}{\partial\nu_k\, \partial\nu_l} \left. Q(\widehat{\boldsymbol{\mu}}(\boldsymbol{\nu}), \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h})\right|_{\nu=\hat{\nu}}.

Details on parameter estimation

Models for water retention curves and hydraulic conductivity functions

Currently, fit_wrc_hcc allows to estimate the parameters of the simplified form of the Van Genuchten-Mualem (VGM) model (Van Genuchten, 1980) with the restriction m = 1 - \frac{1}{n}, see wc_model and hc_model. This model has the following parameters:

• \boldsymbol{\mu}^\mathrm{T} = (\theta_r, \theta_s, K_0) (see above) and

• \boldsymbol{\nu}^\mathrm{T} = (\alpha, n, \tau) where \alpha is the inverse air entry pressure, n the shape and \tau the tortuosity parameter.

Any of these parameters can be either estimated from data or kept fixed at the specified initial values (see arguments param and fit_param of fit_wrc_hcc).

Imposing physical constraints on the estimated parameters

Parameters of models for the water retention curve and the hydraulic conductivity function may vary only within certain bounds (see wc_model, hc_model and param_boundf for allowed ranges). fit_wrc_hcc either estimates transformed parameters that vary over the whole real line and can therefore be estimated without constraints (see param_transf), or it uses algorithms (quadratic programming for estimating \boldsymbol{\mu}^\mathrm{} , nonlinear optimisation algorithms with box constraints for estimating \boldsymbol{\nu}^\mathrm{} ) that restrict estimates to permissible ranges, see Details section of control_fit_wrc_hcc.

In addition, for natural soils, the parameters of the VGM model cannot vary independently from each other over the allowed ranges. Sets of fitted parameters should always result in soil hydraulic quantities that are physically meaningful. One of these quantities is the characteristic length L_\mathrm{c} of stage-I evaporation from a soil (Lehmann et al., 2008). L_\mathrm{c} can be related to the parameters of the VGM model, see Lehmann et al. (2008, 2020) and evaporative-length.

Using several soil hydrological databases, Lehmann et al. (2020) analysed the mutual dependence of VGM parameters and proposed regression equations to relate the inverse air entry pressure \alpha and the saturated hydraulic K_0 to the shape parameter n, which characterises the width of the pore size distribution of a soil. Using these relations, Lehmann et al. (2020) then computed the expected value (“target”) L_\mathrm{t} of L_\mathrm{c} for given n and tortuosity parameter \tau, see evaporative-length. fit_wrc_hcc allows to constrain estimates of the nonlinear parameters \boldsymbol{\nu}^\mathrm{} by defining permissible lower and upper bounds for the ratio L_\mathrm{c}/L_\mathrm{t}, see arguments ratio_lc_lt_bound of fit_wrc_hcc and settings of control_fit_wrc_hcc.

Choice of optimisation algorithm for estimating the nonlinear parameters

To estimate \boldsymbol{\nu}^\mathrm{} , fit_wrc_hcc minimises Q( \widehat{\boldsymbol{\mu}}(\boldsymbol{\nu}), \boldsymbol{\nu}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h}) either by a nonlinear optimisation algorithm available in the library NLopt (Johnson, see nloptr) or by the Shuffled Complex Evolution (SCE) optimisation algorithm (Duan et al., 1994, see SCEoptim). The choice of the algorithm is controlled by the argument settings of the function control_fit_wrc_hcc:

1. global optimisation without constraints for the ratio L_\mathrm{c}/L_\mathrm{t}
   (settings = "uglobal" or settings = "sce"),

2. global optimisation with inequality constraints for the ratio L_\mathrm{c}/L_\mathrm{t}
   (settings = "cglobal"),

3. local optimisation without constraints for the ratio L_\mathrm{c}/L_\mathrm{t}
   (settings = "ulocal"),

4. local optimisation with inequality constraints for the ratio L_\mathrm{c}/L_\mathrm{t}
   (settings = "clocal").

The settings argument also sets reasonable default values for the termination (= convergence) criteria for the various algorithms, see NLopt documentation, section Termination conditions. The NLopt documentation contains a very useful discussion of (constrained) optimisation problems in general, global vs. local optimisation and gradient-based vs. derivative-free algorithms.

Note that besides the settings argument of control_fit_wrc_hcc, the arguments nloptr and sce along with the functions control_nloptr and control_sce allow to fully control the nonlinear optimisation algorithms, see control_fit_wrc_hcc for details.

Computing initial values of parameters

For local optimisation algorithms “good” initial values of \boldsymbol{\nu}^\mathrm{} are indispensable for successful estimation. fit_wrc_hcc allows to compute initial values of \alpha and n for the Van Genuchten model from water retention data by the following procedure:

1. Smooth the water retention data, (\theta_i, y_i= \log(h_i)), i=1,2, ... n_\theta, by an additive model.

2. Determine the saturation, S^*, and the logarithm of capillary pressure head, y^* = \log(h^*), at the inflection point of the additive model fit.

3. Find the root, say \widehat{m}, of S^* = (1 + 1/m)^{-m}. One obtains the right-hand side of this equation by solving \frac{\partial^2}{\partial y^2} \left[ S_\mathrm{VG}(\exp(y); \boldsymbol{\nu}) \right] = 0 for y and plugging the result into the expression for S_\mathrm{VG}(\exp(y); \boldsymbol{\nu}), see wc_model.

4. Compute \widehat{n} = 1 / (1 - \widehat{m}) and \widehat{\alpha} = 1 / \exp(y^*) \,(1/\widehat{m})^{1-\widehat{m}}. The second expression is again a result of solving \frac{\partial^2}{\partial y^2} \left[ S_\mathrm{VG}(\exp(y); \boldsymbol{\nu}) \right] = 0.

Initial values for local optimisation algorithms can of course also be obtained by first estimating the parameters by a global algorithm. These estimates can be “refined” in a second step by a local unconstrained algorithm, possibly followed by a third call of fit_wrc_hcc to constrain the estimated parameters by the ratio L_\mathrm{c}/L_\mathrm{t}. The method coef.fit_wrc_hcc can be used to extract the estimated parameters from an object of class fit_wrc_hcc and to pass them as initial values to fit_wrc_hcc, see fit_wrc_hcc for examples.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Bates, D. M., and Watts, D. G. (1988) Nonlinear Regression Analysis and its Applications. John Wiley & Sons, New York doi:10.1002/9780470316757.

Duan, Q., Sorooshian, S., and Gupta, V. K. (1994) Optimal use of the SCE-UA global optimisation method for calibrating watershed models, Journal of Hydrology 158, 265–284, doi:10.1016/0022-1694(94)90057-4.

Johnson, S.G. The NLopt nonlinear-optimisation package. https://github.com/stevengj/nlopt.

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.

Stewart, W.E., Caracotsios, M. Sørensen, J.P. 1992. Parameter estimation from multiresponse data. AIChE Journal, 38, 641–650, doi:10.1002/aic.690380502.

Stewart, W.E. and Sørensen, J.P. (1981) Bayesian estimation of common parameters from multiresponse data with missing observations. Technometrics, 23, 131–141,
doi:10.1080/00401706.1981.10486255.

Van Genuchten, M. Th. (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44, 892–898, doi:10.2136/sssaj1980.03615995004400050002x.

See Also

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Controlling fit_wrc_hcc

Description

This page documents options to control fit_wrc_hcc. It describes the arguments of the functions control_fit_wrc_hcc, param_boundf, param_transf, fwd_transf, dfwd_transf, bwd_transf, control_nloptr, control_sce and control_pcmp, which all serve to steer fit_wrc_hcc.

Usage

control_fit_wrc_hcc(
    settings = c("uglobal", "ulocal", "clocal", "cglobal", "sce"),
    method = c("ml", "mpd", "wls"), hessian,
    nloptr = control_nloptr(), sce = control_sce(),
    wrc_model = "vg", hcc_model = "vgm",
    initial_param = c(alpha = 2., n = 1.5, tau = 0.5),
    approximation_alpha_k0 =
        c(c0 = 1, c1 = 5.55, c2 = 1.204, c3 = 2.11, c4 = 1.71),
    variable_weight = c(wrc = 1, hcc = 1),
    gam_k = 6, gam_n_newdata = 101, precBits = 256,
    min_nobs_wc = 5, min_nobs_hc = 5,
    keep_empty_fits = FALSE,
    param_bound = param_boundf(), param_tf = param_transf(),
    fwd_tf = fwd_transf(), deriv_fwd_tfd = dfwd_transf(), bwd_tf = bwd_transf(),
    pcmp = control_pcmp())

param_boundf(alpha = c(0 + 10 * sqrt(.Machine$double.eps), 500),
    n = c(1 + 10 * sqrt(.Machine$double.eps), 20), tau = c(-2, 20),
    thetar = c(0, 1), thetas = c(0, 1), k0 = c(0, Inf))

param_transf(alpha = c("log", "identity"), n = c("log1", "identity"),
    tau = c("logitlu", "identity"), k0 = c("log", "identity"))

fwd_transf(...)

dfwd_transf(...)

bwd_transf(...)

control_nloptr(...)

control_sce(reltol = 1.e-8, maxtime = 20, ...)

control_pcmp(ncores = detectCores() - 1L,
    fork = !identical(.Platform[["OS.type"]], "windows"))

Arguments

Details

Enforcing bounds on the estimated parameters

Parameters of models for the water retention curve and the hydraulic conductivity function may vary only within certain bounds (see param_boundf for allowed ranges). fit_wrc_hcc uses two mechanisms to constrain parameter estimates to permissible ranges:

1. Parameter transformations

   If a local algorithm is used for nonlinear optimisation (settings = "ulocal" or settings = "clocal") and a transformation not equal to "identity" is specified in param_tf for any of the nonlinear parameters \boldsymbol{\nu}^\mathrm{} , then the elements of \boldsymbol{\nu}^\mathrm{} are transformed by the functions given in param_tf. The values of the transformed parameters vary then over the whole real line, and an unconstrained algorithm can be used for nonlinear optimisation.

   Note that the linear parameters \theta_r (residual) and \theta_s (saturated water content) are never transformed and for the saturated hydraulic conductivity, K_0, only "log" (default) or "identity" can be chosen. Quadratic programming (see solve.QP) is employed to enforce the box constraints specified in the argument param_bound for \theta_r and \theta_s. Quadratic programming is also used to enforce the positivity constraint on K_0 if K_0 is not log-transformed ("identity"). Otherwise, the logarithm of K_0 is estimated unconstrainedly, see soilhypfitIntro for further details.

2. Box constraints

   If a global algorithm is used for the optimisation (settings equal to "uglobal" "cglobal" or "sce") or if "identity" transformations are specified for all elements of \boldsymbol{\nu}^\mathrm{} , then an optimisation algorithm is deployed that respects the box constraints given in param_bound. If parameters are estimated for several soil samples in a single call of fit_wrc_hcc and if sample-specific box constraints should be used then the lower and upper bounds of the box-constraints must be specified by the arguments lower_param and upper_param of the function fit_wrc_hcc, see explanations there.

   Further note that the transformations specified by param_tf for the nonlinear parameters \boldsymbol{\nu}^\mathrm{} are ignored when a global optimisation algorithm is used.

Parameter transformations

The arguments param_tf, fwd_tf, deriv_fwd_tfd, bwd_tf define how the model parameters are transformed for estimation by local optimisation algortihms (see above and soilhypfitIntro). The following transformations are currently available:

"log":

\log(x),

"log1":

\log(x-1),

"logitlu":

\log((x - l) / (u - x)) with l and u the allowed lower and upper bounds for a parameter, see param_boundf,

"identity":

no transformation.

These are the possible values that the various arguments of the function param_transf accept (as quoted character strings), and these are the names of the list components returned by fwd_transf, dfwd_transf and bwd_transf.

Additional transformations can be implemented by:

1. Extending the function definitions by arguments like

   fwd_tf = fwd_transf(my_fun = function(x) your transformation),
   deriv_fwd_tfd = dfwd_transf(my_fun = function(x) your derivative),
   bwd_tf = bwd_transf(my_fun = function(x) your back-transformation),

2. Assigning to a given argument of param_transf the name of the new function, e.g.
   alpha = "my_fun".

Note that the values given to the arguments of param_transf must match the names of the functions returned by fwd_transf, dfwd_transf and bwd_transf.

High-precision numerical computations

Estimation of \log(K_0) is somewhat delicate for large values of the shape parameter n and/or small values of \alpha. The water saturation and the relative conductivity are then close to zero for capillary pressure head exceeding 1/\alpha. To avoid numerical problems caused by limited accuracy of double precision numbers, fit_wrc_hcc uses the function mpfr of the package Rmpfr for high-accuracy computations. The argument precBits of control_fit_wrc_hcc controls the accuracy. Increase its value if computation fails with a respective diagnostic message.

Options to choose the approach for nonlinear optimisation

The argument settings defines sets of default options to control the optimisers. The following settings are currently defined:

"uglobal":

unconstrained optimisation by any of the global algorithms (named "NLOPT_G...") of the NLopt library.

"cglobal":

constrained optimisation by the global algorithm "NLOPT_GN_ISRES" of NLopt that allows for inequality constraints.

"ulocal":

unconstrained optimisation by any of the local algorithms
(named "NLOPT_L...") of NLopt.

"clocal":

constrained optimisation by any of the local algorithms
("NLOPT_LN_COBYLA", "NLOPT_LN_AUGLAG", "NLOPT_LD_AUGLAG", "NLOPT_LD_SLSQP",
"NLOPT_LD_MMA"), "NLOPT_LD_CCSAQ") of NLopt that allow for inequality constraints.

"sce":

unconstrained optimisation by the global algorithm implemented in SCEoptim.

The functions control_nloptr and control_sce allow finer control of the optimisers. control_nloptr and control_sce take any argument available to pass controlling options to the optimisers nloptr (by its argument opts) and SCEoptim (by its argument control), respectively.

Controlling nloptr

The function nloptr.print.options prints all options available to control nloptr by its argument opts. Detailed information on the options can be found in the NLopt documentation.

The function control_fit_wrc_hcc sets meaningful defaults for opts in dependence of the chosen optimisation approach as specified by the argument settings, and it checks the consistency of the arguments of control_nloptr if they are explicitly passed to fit_wrc_hcc.

The following defaults are set by control_fit_wrc_hcc for the argument opts of nloptr (:

1. Unconstrained, global optimisation (settings = "uglobal"):

   nloptr = control_nloptr(
     algorithm = "NLOPT_GN_MLSL_LDS",
     local_opts = list(
       algorithm = "NLOPT_LN_BOBYQA",
       xtol_rel = -1.,
       ftol_rel = 1.e-6
     ),
     xtol_rel = -1,
     ftol_rel = -1,
     maxeval = 125,
     maxtime = -1)
           

   In addition, any parameter transformations specified by param_tf are overridden and the untransformed parameters ("identity") are estimated when settings = "uglobal" is chosen.

2. Constrained, global optimisation (settings = "cglobal"):

   nloptr = control_nloptr(
     algorithm = "NLOPT_GN_ISRES",
     xtol_rel = -1,
     ftol_rel = -1,
     maxeval = 1000,
     maxtime = -1)
           

   In addition, any parameter transformations specified by param_tf are overridden and the untransformed parameters ("identity") are estimated when settings = "cglobal" is chosen.

3. Unconstrained, local optimisation (settings = "ulocal"):

   nloptr = control_nloptr(
     algorithm = "NLOPT_LN_BOBYQA",
     xtol_rel = -1,
     ftol_rel = 1.e-8,
     maxeval = 250,
     maxtime = -1)
           

4. Constrained, local optimisation (settings = "clocal"):

   nloptr = control_nloptr(
     algorithm = "NLOPT_LD_CCSAQ",
     xtol_rel = -1,
     ftol_rel = 1.e-8,
     maxeval = 1000,
     maxtime = -1)
           

   If the algorithm "NLOPT_LD_AUGLAG" is used for constrained, local optimisation then

   nloptr = control_nloptr(
     algorithm = "NLOPT_LD_AUGLAG",
     local_opts = list(
       algorithm = "NLOPT_LD_LBFGS",
       xtol_rel = -1.,
       ftol_rel = 1.e-6
     ),
     xtol_rel = -1,
     ftol_rel = 1.e-8,
     maxeval = 1000,
     maxtime = -1)
           

For other, unspecified elements of opts default values as listed by nloptr.print.options are used.

Controlling SCEoptim

The function control_sce sets meaningful defaults for the argument control of SCEoptim. Currently, the following defaults are defined:

sce = control_sce(
  reltol = 1e-08,
  maxtime = 20)
      

In addition, any parameter transformations specified by param_tf are overridden and the untransformed parameters ("identity") are estimated when settings = "sce" is chosen.

Value

control_fit_wrc_hcc creates a list with components settings, hessian, method, nloptr, sce, wrc_model, hcc_model, initial_param, approximation_alpha_k0, variable_weight, gam_k, gam_n_newdata, precBits, min_nobs_wc, min_nobs_hc, keep_empty_fits, param_bound, param_tf, fwd_tf, deriv_fwd_tfd, bwd_tf, pcmp corresponding to its arguments and some further components (delta_sat_0, grad_eps, use_derivative) that cannot be changed by the user.

control_nloptr and control_sce create lists with control parameters passed to nloptr and SCEoptim, respectively, see Details.

param_boundf generates a list with allowed lower and upper bounds of the model parameters.

param_transf generates a list with keywords that define what transformations are used for estimating the model parameters, and fwd_transf, bwd_transf and dfwd_transf return lists of functions with forward and backward transformations and the first derivatives of the forward transformations, see Details.

control_pcmp generates a list with control parameters for parallel computations.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Johnson, S.G. The NLopt nonlinear-optimisation package. https://github.com/stevengj/nlopt.

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.

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;

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

# use of \donttest{} because execution time exceeds 5 seconds
data(sim_wrc_hcc)

# estimate parameters for a single soil sample by maximizing loglikelihood ...

# ... with unconstrained, global optimisation algorithm NLOPT_GN_MLSL
coef(
  fit1 <- fit_wrc_hcc(
    wrc_formula = wc ~ head, hcc_formula = hc ~ head,
    data = sim_wrc_hcc, subset = id == 2
  ), gof = TRUE)

# ... as fit1 but fitting parameter tau as well
coef(
  fit2 <- update(fit1,
    fit_param = default_fit_param(tau = TRUE)
  ), gof = TRUE)

plot(fit1, y = fit2)

# ... with unconstrained, local optimisation algorithm NLOPT_LN_BOBYQA,
#     initial values for alpha and n are computed from data and
#     transformed nonlinear parameters are estimated without box-constraints
coef(
  fit3 <- update(
    fit2,
    control = control_fit_wrc_hcc(settings = "ulocal"),
    verbose = 2), gof = TRUE)

# estimate parameters by unconstrained, weighted least squares minimisation with
#     algorithm NLOPT_LD_LBFGS, giving larger weight to conductivity data,
#     using specified initial values for alpha and n and
#     fitting untransformed nonlinear parameters with default box constraints
#     defined by param_boundf()
#     diagnostic output directly from nloptr
coef(
  fit4 <- update(
    fit2,
    param = c(alpha = 1.7, n = 2),
    control = control_fit_wrc_hcc(
      settings = "ulocal", method = "wls",
      variable_weight = c(wrc = 1, hcc = 2),
      nloptr = control_nloptr(algorithm = "NLOPT_LD_LBFGS", print_level = 3),
      param_tf = param_transf(alpha = "identity", n = "identity", tau = "identity")
    ), verbose = 0), gof = TRUE)

# ... as fit4 but giving even larger weight to conductivity data
coef(
  fit5 <- update(
    fit4,
    control = control_fit_wrc_hcc(
      settings = "ulocal", method = "wls",
      variable_weight = c(wrc = 1, hcc = 5),
      nloptr = control_nloptr(algorithm = "NLOPT_LD_LBFGS", print_level = 3),
      param_tf = param_transf(alpha = "identity", n = "identity", tau = "identity")
    ), verbose = 0), gof = TRUE)

plot(fit4, y = fit5)

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

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 papritz@retired.ethz.ch.

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

Parametric Modelling of Soil Hydraulic Properties

Description

The function fit_wrc_hcc estimates parameters of models for the soil water retention curve and/or soil hydraulic conductivity function from respective measurements by nonlinear regression methods, optionally subject to physical constraints on the estimated parameters. fit_wrc_hcc uses optimisation algorithms of the NLopt library (Johnson, see nloptr-package) and the Shuffled Complex Evolution (SCE) algorithm (Duan et al., 1994) implemented in the function SCEoptim.

Usage

fit_wrc_hcc(
    wrc_formula, hcc_formula, data,
    subset = NULL, wrc_subset = subset, hcc_subset = subset,
    weights = NULL, wrc_weights = weights, hcc_weights = weights,
    na.action, param = NULL, lower_param = NULL, upper_param = NULL,
    fit_param = default_fit_param(),
    e0 = 2.5e-3, ratio_lc_lt_bound = c(lower = 0.5, upper = 2),
    control = control_fit_wrc_hcc(), verbose = 0)

default_fit_param(
    alpha = TRUE, n = TRUE, tau = FALSE,
    thetar = TRUE, thetas = TRUE, k0 = TRUE)

Arguments

Details

Estimating model parameters of water retention curves and/or hydraulic conductivity functions

The function fit_wrc_hcc estimates model parameters of water retention curves and/or hydraulic conductivity functions by maximum likelihood (ml, default) maximum posterior density (mpd) or nonlinear least squares (wls), see argument method of control_fit_wrc_hcc. Measurements must be available for at least one of the two material functions. If one type of data is missing then the respective formula, subset and weights arguments should be omitted.

If both types of data are available then measurements are weighed when computing the residual sum of squares for wls, but unit weights are used by default for mpd and ml estimation, see soilhypfitIntro. The wls weights are the product of the case weights w^\prime_{\theta,i} and w^\prime_{K,i} given by weights, wrc_weights and hcc_weights, respectively, and the variable weights as specified by the argument variable_weight of control_fit_wrc_hcc. Note that the variable_weights are not used “as is”, but they are multiplied by the inverse variances of the water content or the (log-transformed) conductivity data per sample to obtain the variable weights w_\theta, w_K mentioned in soilhypfitIntro.

Estimating model parameters for a single or multiple soil samples

When parameters are estimated for a single soil sample only, then model formulae are of the form wc ~ head and/or hc ~head, and there is no need to specify a sample id. In this case param, lower_param, upper_param, fit_param and ratio_lc_lt_bound are vectors and e0 is a scalar.

fit_wrc_hcc allows to fit models to datasets containing data for multiple soil samples. The model formula must then be of the form wc ~ head|id and/or hc ~ head|id where the factor id identifies the soil samples. If param, lower_param, upper_param, fit_param and ratio_lc_lt_bound are vectors and e0 is a scalar then this information is used for processing all the soil samples. If specific information per sample should be used then param, lower_param, upper_param, fit_param and ratio_lc_lt_bound must be matrices (or data frames) with as many rows as there are soil samples. The matrices (or data.frames) must have rownames matching the levels of the factor id defining the soil samples. e0 must be a named vector with as many elements as there are soil samples and the names must again match the levels of id.

By default, fit_wrc_hcc processes data of multiple soil samples in parallel, see control_pcmp for options to control parallel computing. Note that diagnostic output (except warning messages) is suppressed in parallel computations. If computations fail for a particular soil sample, the id of the sample with the failed fit can be extracted by the utility function select_failed_fits and the respective error messages by extract_error_messages.

Controlling fit_wrc_hcc

The argument control is used to pass a list of options to fit_wrc_hcc that steer the function, see soilhypfit-package and control_fit_wrc_hcc for full details.

Value

fit_wrc_hcc returns an object of class fit_wrc_hcc, which is a list with the following components:

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Duan, Q., Sorooshian, S., and Gupta, V. K. (1994) Optimal use of the SCE-UA global optimisation method for calibrating watershed models, Journal of Hydrology 158, 265–284, doi:10.1016/0022-1694(94)90057-4.

Johnson, S.G. The NLopt nonlinear-optimisation package. https://github.com/stevengj/nlopt.

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.

See Also

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

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

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

data(sim_wrc_hcc)

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


# estimate parameters for single sample ...

#  ... from wrc and hcc data
plot(rfit_wrc_hcc <- fit_wrc_hcc(
  wrc_formula = wc ~ head, hcc_formula = hc ~ head,
  data = sim_wrc_hcc, subset = id == 1))
coef(rfit_wrc_hcc, gof = TRUE)

#  ... from wrc data
plot(rfit_wrc <- fit_wrc_hcc(
  wrc_formula = wc ~ head,
  data = sim_wrc_hcc, subset = id == 1))
coef(rfit_wrc, gof = TRUE)

#  ... from  hcc data
plot(rfit_hcc <- fit_wrc_hcc(
  hcc_formula = hc ~ head,
  data = sim_wrc_hcc, subset = id == 1))
coef(rfit_hcc, gof = TRUE)

# ... from wrc and hcc data
#     keeping some parameters fixed
plot(rfit_wrc_hcc_fixed <- fit_wrc_hcc(
  wrc_formula = wc ~ head, hcc_formula = hc ~ head,
  data = sim_wrc_hcc, subset = id == 1,
  param = c(alpha = 2.1, thetar = 0.1),
  fit_param = default_fit_param(alpha = FALSE, thetar = FALSE)),
  y = rfit_wrc_hcc)
coef(rfit_wrc_hcc, gof = TRUE)
coef(rfit_wrc_hcc_fixed, gof = TRUE)


# estimate parameters for 3 samples simultaneously by ...

# ... unconstrained, global optimisation algorithm NLOPT_GN_MLSL (default)
rfit_uglob <- fit_wrc_hcc(
  wrc_formula = wc ~ head | id, hcc_formula = hc ~ head | id,
  data = sim_wrc_hcc,
  control = control_fit_wrc_hcc(pcmp = control_pcmp(ncores = ncpu)))
summary(rfit_uglob)
op <- par(mfrow = c(3, 2))
plot(rfit_uglob)
on.exit(par(op))

# ... unconstrained, global optimisation algorithm SCEoptim
rfit_sce <- update(
  rfit_uglob,
  control = control_fit_wrc_hcc(
    settings = "sce", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_sce, gof = TRUE, lc = TRUE)
convergence_message(2, sce = TRUE)
op <- par(mfrow = c(3, 2))
plot(rfit_sce, y = rfit_uglob)
on.exit(par(op))

# ... unconstrained, local optimisation algorithm NLOPT_LN_BOBYQA,
#     logging iteration results to console
rfit_uloc <- update(
  rfit_uglob,
  param = as.matrix(coef(rfit_uglob, what = "nonlinear")),
  control = control_fit_wrc_hcc(
    settings = "ulocal", pcmp = control_pcmp(ncores = 1L)),
  verbose = 2)
coef(rfit_uloc, gof = TRUE, lc = TRUE)

# ... constrained, global optimisation algorithm NLOPT_GN_ISRES
rfit_cglob <- update(
  rfit_uglob,
  ratio_lc_lt_bound = c(lower = 0.8, upper = 1.2),
  control = control_fit_wrc_hcc(
    settings = "cglobal", nloptr = control_nloptr(ranseed = 1),
    pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cglob, gof = TRUE, lc = TRUE)

# ... constrained, local optimisation algorithm NLOPT_LD_CCSAQ
#     starting from unconstrained, locally fitted initial values
rfit_cloc_1 <- update(
  rfit_uglob,
  param = coef(rfit_uloc, what = "nonlinear"),
  ratio_lc_lt_bound = c(lower = 0.8, upper = 1.2),
  control = control_fit_wrc_hcc(
    settings = "clocal", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cloc_1, gof = TRUE, lc = TRUE)
op <- par(mfrow = c(3, 2))
plot(x = rfit_uloc, y = rfit_cloc_1)
on.exit(par(op))

# ... constrained, local optimisation algorithm NLOPT_LD_CCSAQ
#     starting from constrained, globally fitted initial values,
#     using sample-specific evaporation rates and limits for ratio lc/lt
rfit_cloc_2 <- update(
  rfit_uglob,
  param = as.matrix(coef(rfit_cglob, what = "nonlinear")),
  e0 = c("1" = 0.002, "2" = 0.0025, "3" = 0.003),
  ratio_lc_lt_bound = rbind(
    "1" = c(lower = 0.7, upper = 2),
    "2" = c(lower = 0.8, upper = 1.4),
    "3" = c(lower = 0.8, upper = 2)
  ),
  control = control_fit_wrc_hcc(
    settings = "clocal", pcmp = control_pcmp(ncores = ncpu)))
coef(rfit_cloc_2, gof = TRUE, lc = TRUE, e0 = TRUE)

# ... global optimisation algorithm NLOPT_GN_MLSL
#     with sample-specific box-constraints

rfit_uglob_bc <- update(
  rfit_uglob,
  lower_param = rbind(
    "1" = c(alpha = 2.4, n = 1.11, thetar = 0.2, thetas = 0.6),
    "2" = c(alpha = 1.2, n = 1.12, thetar = 0.2, thetas = 0.6),
    "3" = c(alpha = 1.2, n = 1.13, thetar = 0.2, thetas = 0.6)
  ),
  upper_param = rbind(
    "1" = c(alpha = 20.1, n = 2.51, thetar = 0.6, thetas = 0.6),
    "2" = c(alpha = 20.2, n = 1.5,  thetar = 0.6, thetas = 0.6),
    "3" = c(alpha = 1.3, n = 2.53,  thetar = 0.6, thetas = 0.6)
  )
)
coef(rfit_uglob, gof = TRUE)
coef(rfit_uglob_bc, gof = TRUE)

Models for Soil Hydraulic Conductivity Functions

Description

The functions hc_model and hcrel_model compute, for given capillary pressure head h, the hydraulic conductivity K(h) and the relative hydraulic conductivity k(h) respectively, of a soil by parametrical models.

Usage

hcrel_model(h, nlp, precBits = NULL, hcc_model = "vgm")

hc_model(h, nlp, lp, precBits = NULL, hcc_model = "vgm")

Arguments

Details

The functions hcrel_model and hc_model currently model soil hydraulic conductivity functions only by the simplified form of the Van Genuchten-Mualem model (Van Genuchten, 1980) with the restriction m = 1 - \frac{1}{n}, i.e. by

k_\mathrm{VGM}(h; \boldsymbol{\nu}) = S_\mathrm{VG}(h; \boldsymbol{\nu})^\tau \, \left[ 1 - \left( 1 - S_\mathrm{VG}(h; \boldsymbol{\nu})^\frac{n}{n-1} \right)^\frac{n-1}{n} \right]^2,

K_\mathrm{VGM}(h; \mbox{$\mu$}, \boldsymbol{\nu}) = K_0 \, k_\mathrm{VGM}(h; \boldsymbol{\nu}),

where \mu = K_0 is the saturated hydraulic conductivity (K_0 > 0), \boldsymbol{\nu}^\mathrm{T} = (\alpha, n, \tau) are the inverse air entry pressure (\alpha > 0), the shape (n > 1) and tortuosity parameter (\tau > -2), and S_\mathrm{VG}(h; \boldsymbol{\nu}) is the the volumetric water saturation, see sat_model for an expression.

Note that \mu and \boldsymbol{\nu}^\mathrm{} are passed to the functions by the arguments lp and nlp, respectively.

Value

A numeric vector with values of (relative) hydraulic conductivity.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Van Genuchten, M. Th. (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44, 892–898, doi:10.2136/sssaj1980.03615995004400050002x.

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

## define capillary pressure head (unit meters)
h <- c(0.01, 0.1, 0.2, 0.3, 0.5, 1., 2., 5.,10.)

## compute (relative) hydraulic conductivity
hcrel <- hcrel_model(h, nlp = c(alpha = 1.5, n = 2, tau = 0.5))
hc <- hc_model(h, nlp = c(alpha = 1.5, n = 2, tau = 0.5), lp = c(k0 = 5))

## display hydraulic conductivity function
op <- par(mfrow = c(1, 2))
plot(hcrel ~ h, log = "xy", type = "l")
plot(hc ~ h, log = "xy", type = "l")
on.exit(par(op))

Internal Functions of Package soilhypfit

Description

The unexported internal functions

• d1ineq_constraint_nlp

• d1l0

• d1lc

• estimate_lp

• estimate_nlp

• fit_wrc_hcc_fit

• gradient_nlp

• ineq_constraint_nlp

• model_names_nlp_lp

• model_fit_param_consistent

• model_param_tf_nlp_identity

• objective_nlp

• stop_cluster

are not intended for direct use. However, as any unexported function, they can be accessed by typing soilhypfit:::function-name.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Profile Loglikelihoods and Confidence Intervals for Parametric Modelling of Soil Hydraulic Properties

Description

The function prfloglik_sample computes for a single soil sample a loglikelihood profile as a function of the specified values for subsets of model parameters of soil water retention and/or soil hydraulic conductivity functions. The function confint_prfloglik_sample computes a confidence interval of one model parameter based on the likelihood ratio test for a single soil sample, and the S3 method confint computes confidence intervals of nonlinear model parameters for multiple soil samples.

Usage

prfloglik_sample(object, values, soil_sample,
    ncores = min(detectCores() - 1L, NROW(values)), verbose = 0)

confint_prfloglik_sample(object, parm = names(default_fit_param()),
    soil_sample, level = 0.95, test = c("F", "Chisq"),
    denominator_df = c("nonlinear", "all"), param_bound = NULL,
    root_tol = .Machine$double.eps^0.25, froot_tol = sqrt(root_tol),
    verbose = 0)

## S3 method for class 'fit_wrc_hcc'
confint(object,
    parm = names(object[["control"]][["initial_param"]]), level = 0.95,
    subset = NULL, type = c("loglik", "normal"), test = c("F", "Chisq"),
    denominator_df = c("nonlinear", "all"),
    root_tol = .Machine$double.eps^0.25, froot_tol = sqrt(root_tol),
    ncores = detectCores() - 1L, verbose = 0, ...)

Arguments

Details

Computing likelihood profiles

The function prfloglik_sample computes the loglikelihood profile for a subset of the nonlinear (\boldsymbol{\nu}^\mathrm{} ) (and linear \boldsymbol{\mu}^\mathrm{} ) model parameters. We denote the profiling parameters of interest by \boldsymbol{\phi} and the nuisance parameters that are estimated when computing the profile by \boldsymbol{\psi}, so that ( \boldsymbol{\phi}^\mathrm{T}, \boldsymbol{\psi}^\mathrm{T} )^\mathrm{T} is a rearranged version of the parameter vector ( \boldsymbol{\mu}^\mathrm{T}, \boldsymbol{\nu}^\mathrm{T} )^\mathrm{T} , see soilhypfitIntro. prfloglik_sample computes the estimates \widehat{\boldsymbol{\psi}}(\boldsymbol{\phi}) of \boldsymbol{\psi} and the profile loglikelihood Q(\boldsymbol{\phi}, \widehat{\boldsymbol{\psi}}(\boldsymbol{\phi}); \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) as a function of \boldsymbol{\phi}.

To compute the estimates, \boldsymbol{\psi} is partitioned into the nonlinear and conditionally linear parameters, see soilhypfitIntro for details. prfloglik_sample uses the packages parallel and snowfall for parallelized computation of loglikelihood profiles.

Computing likelihood based confidence intervals

The function confint_prfloglik_sample computes the confidence interval of a single model parameter \phi \in ( \boldsymbol{\mu}^\mathrm{T}, \boldsymbol{\nu}^\mathrm{T} )^\mathrm{T} based on the likelihood ratio test. The interval is computed by assuming either

• that the test statistic T(\phi) = \Delta Q(\phi) = 2( Q(\widehat{\phi}, \widehat{\boldsymbol{\psi}}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) - Q(\phi, \widehat{\boldsymbol{\psi}}(\phi); \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} )) follows a \chi^2_1-distribution with 1 degrees of freedom (possible choice when both water retention and/or hydraulic conductivity data were used to estimate the parameters), or

• that the transformed test statistic T(\phi) = \left(\exp(\Delta Q(\phi) / n) - 1 \right) (n - p) follows an F(1,n-p)-distribution where n is the number of water content or hydraulic conductivity measurements, and p is the number of estimated parameters (see Uusipaikka, 2008, p. 115, for details). denominator_df allows to control how p is chosen. Note that this test distribution can only be chosen (and is then the default) when only water retention or only hydraulic conductivty data were used to estimate the parameters.

confint_prfloglik_sample computes profile loglikelihoods Q(\phi, \widehat{\boldsymbol{\psi}}(\phi); \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) by
prfloglik_sample and then uses uniroot to search the roots of the equation

f(\phi) = q_T(\gamma)- T(\phi)

in the interval defined by param_bound. q_T(\gamma) is the \gamma-quantile of the chosen test distribution for T.

Computing confidence intervals for several soil samples

The confint method computes 2 types of confidence intervals (in dependence of type) of only the nonlinear parameters \boldsymbol{\nu}^\mathrm{} , possibly for multiple soil samples at the same time:

• intervals based on the asymptotic normal distribution of maximum likelihood estimates with standard errors computed by the vcov method for class fit_wrc_hcc (see vcov.fit_wrc_hcc,

• intervals based on the likelihood ratio test by confint_prfloglik_sample. The intervals for several soil samples are computed in parallel.

Requirements for computing loglikelihood profiles

The parameters contained in object must be estimated by maximum likelihood (method = "ml", see soilhypfitIntro and control_fit_wrc_hcc. Use of other estimation methods results an error.

Preferably an unconstrained local algorithm (settings = "ulocal", see soilhypfitIntro and control_fit_wrc_hcc)) is used for minimizing the negative loglikelihood when generating object. Use of other algorithms generates a warning message.

Value

The function prfloglik_sample returns a data.frame with the columns of values (parameters of interest \boldsymbol{\phi}), a column loglik with the maximized profile loglikelihood Q(\phi, \widehat{\boldsymbol{\psi}}(\phi); \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) , columns with the estimated parameters \widehat{\boldsymbol{\psi}}(\boldsymbol{\phi}) , columns with the gradient of the loglikelihood with respect to the estimated nonlinear parameters (missing if all nonlinear parameters were fixed) and a column converged, indicating whether convergence has occurred (see convergence_message) when estimating the nonlinear parameters (equal to NA when all nonlinear parameters are fixed).

The function confint_prfloglik_sample returns a numeric vector of length 2 with the lower and upper limits of the confidence interval. If no roots were found then NA is returned. The returned result further contains as attribute list prfloglik the parameter estimate \widehat{\phi} (param_estimate), the maximized loglikelihood Q(\widehat{\phi}, \widehat{\boldsymbol{\psi}}; \boldsymbol{\theta}, \boldsymbol{K}, \boldsymbol{h} ) (loglik), the quantile of the test distribution q_{\mathrm{test}}(\gamma) (qtest), the type of test distribution used (test), the significance level, the number of water content (nobs_wrc) and conductivity measurements (nobs_hcc) and the function values of f(\phi) evaluated at the roots returned by uniroot.

The method confint returns a dataframe with the lower and upper limits of the confidence intervals for the estimated nonlinear parameters.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Uusipaikka, E. (2008) Confidence Intervals in Generalized Regression Models. Chapman & Hall/CRC Press, Boca Raton doi:10.1201/9781420060386.

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;

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

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

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

if(!requireNamespace("lattice", quietly = TRUE)){
  install.packages("lattice", repos = "https://cloud.r-project.org/")
  requireNamespace("lattice", quietly = TRUE)
}

data(sim_wrc_hcc)

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

# estimate parameters for 3 samples simultaneously by ...

# ... unconstrained, global optimisation algorithm NLOPT_GN_MLSL (default)
rfit_uglob <- fit_wrc_hcc(
  wrc_formula = wc ~ head | id, hcc_formula = hc ~ head | id,
  data = sim_wrc_hcc,
  control = control_fit_wrc_hcc(param_bound = param_boundf(
      alpha = c(0.00001, 50), n = c(1.0001, 7), tau = c(-1, 5)
    ), pcmp = control_pcmp(ncores = ncpu)))

# ... unconstrained, local optimisation algorithm NLOPT_LN_BOBYQA,
rfit_uloc <- update(
  rfit_uglob,
  param = as.matrix(coef(rfit_uglob, what = "nonlinear")),
  control = control_fit_wrc_hcc(
    settings = "ulocal", param_tf = param_transf(
      alpha = "identity", n = "identity", tau = "identity"
    ), param_bound = param_boundf(
      alpha = c(0.00001, 50), n = c(1.0001, 7), tau = c(-1, 5)
    ), pcmp = control_pcmp(ncores = ncpu)))

# extract estimated parameters for sample id == "1"
coef_id_1 <- unlist(coef(rfit_uloc, gof = TRUE, se = TRUE, subset = "1"))

# compute loglikelihood profile of parameter alpha for sample id == "1"
rprfloglik_alpha <- prfloglik_sample(
  rfit_uloc, values = data.frame(alpha = seq(1.5, 3.0, length = 40L)),
  soil_sample = "1", ncores = ncpu)
# plot loglikelihood profile along with 95% confidence intervals
plot(loglik ~ alpha, rprfloglik_alpha, type = "l")
abline(v = coef_id_1["alpha"])
# 95% confidence intervall based on likelihood ratio test
abline(h = -coef_id_1["obj"] - 0.5 * qchisq(0.95, 1), lty = "dashed")
# 95% confidence intervall based on asymptotic normal distribution
segments(
  x0 = coef_id_1["alpha"] + qnorm(0.025) * coef_id_1["se.alpha"],
  x1 = coef_id_1["alpha"] + qnorm(0.975) * coef_id_1["se.alpha"],
  y0 = min(rprfloglik_alpha$loglik)
)

# compute loglikelihood profile of parameter n for sample id == "1"
rprfloglik_n <- prfloglik_sample(
  rfit_uloc, values = data.frame(n = seq(1.7, 2.25, length = 40L)),
  soil_sample = "1", ncores = ncpu
)
# plot loglikelihood profile along with 95% confidence intervals
plot(loglik ~ n, rprfloglik_n, type = "l")
abline(v = coef_id_1["n"])
# 95% confidence intervall based on likelihood ratio test
abline(h = -coef_id_1["obj"] - 0.5 * qchisq(0.95, 1), lty = "dashed")
# 95% confidence intervall based on asymptotic normal distribution
segments(
  x0 = coef_id_1["n"] + qnorm(0.025) * coef_id_1["se.n"],
  x1 = coef_id_1["n"] + qnorm(0.975) * coef_id_1["se.n"],
  y0 = min(rprfloglik_n$loglik)
)

# compute loglikelihood profile of parameters alpha and n for sample id == "1"
rprfloglik_alpha_n <- prfloglik_sample(
  rfit_uloc, values = expand.grid(
    alpha = seq(1.5, 3.0, length = 40L), n = seq(1.7, 2.25, length = 40L)),
  soil_sample = "1", ncores = ncpu
)

# joint confidence region of alpha and n based on likelihood ratio test
lattice::levelplot(loglik ~ alpha + n, rprfloglik_alpha_n,
  panel = function(x, y, z, subscripts, ...){
    lattice::panel.levelplot(x, y, z, subscripts, ...)
    lattice::panel.levelplot(x, y, z, subscripts, region = FALSE, contour = TRUE,
      at = -coef_id_1["obj"] - 0.5 * qchisq(0.95, 2),
      lty = "solid"
    )
    lattice::panel.levelplot(x, y, z, subscripts, region = FALSE, contour = TRUE,
      at = -coef_id_1["obj"] - 0.5 * qchisq(0.9, 2),
      lty = "dashed"
    )
    lattice::panel.levelplot(x, y, z, subscripts, region = FALSE, contour = TRUE,
      at = -coef_id_1["obj"] - 0.5 * qchisq(0.8, 2),
      lty = "dotted"
    )
    lattice::panel.points(coef_id_1["alpha"], coef_id_1["n"], pch = 16)
    lattice::panel.lines(
      x = rprfloglik_alpha[, c("alpha", "n")], col = "blue"
    )
    lattice::panel.lines(
      x = rprfloglik_n[, c("alpha", "n")], col = "magenta"
    )
  }, key = list(
    corner = c(1, 0.97),
    lines = list(
      lty = c(rep("solid", 3), "dashed", "dotted"),
      col = c("blue", "magenta", rep("black", 3))
    ),
    text = list(c(
        "estimate of n as a function of fixed alpha",
        "estimate of alpha as a function of fixed n",
        "95% joint confidence region",
        "90% joint confidence region",
        "80% joint confidence region"
      ))
  ))

# compute 95%-confidence interval
(ci.alpha <- confint_prfloglik_sample(
  rfit_uloc, parm = "alpha", soil_sample = "1"
))
# use limits to draw loglikelihood profile for alpha
rprfloglik_alpha <- prfloglik_sample(
  rfit_uloc, values = data.frame(
    alpha = seq(0.9 * ci.alpha[1], 1.1 * ci.alpha[2], length = 40L)),
  soil_sample = "1", ncores = ncpu)
plot(loglik ~ alpha, rprfloglik_alpha, type = "l")
lines(
  ci.alpha,
  rep(diff(unlist(attr(ci.alpha, "prfloglik")[c("qtest", "loglik")])) , 2)
)
abline(v = attr(ci.alpha, "prfloglik")["estimate"], lty = "dashed")

# 95%-confidence intervals of all nonlinear parameters based for all
# soil samples asymptotic normal distribution of maximum likelihood estimates
confint(rfit_uloc, type = "normal")

# 95%-confidence intervals of all nonlinear parameters based for all
# soil samples based on likelihood ratio test
confint(rfit_uloc, ncores = ncpu)

Simulated Soil Water Retention and Hydraulic Conductivity Data

Description

The data give simulated values of water content and hydraulic conductivity at given capillary pressure head for 3 soil samples.

Usage

data(sim_wrc_hcc)

Format

A data frame with 28 observations on the following 4 variables.

id

a factor with levels 1, 2, 3 coding the soil samples.

head

a numeric vector with values of capillary pressure head (unit \mathrm{m}).

wc

a numeric vector with values of simulated volumetric water content (unit -)

hc

a numeric vector with values of simulated hydraulic conductivity (unit \mathrm{m} \, \mathrm{d}^{-1}).

Details

The values of wc and hc were simulated by the model of Van Genuchten Mualem (Van Genuchten, 1980, see wc_model and hc_model) with the following parameters:

Normally distributed errors were added to the model values (wc: sd: 0.05; log(hc): sd 0.5).

References

Van Genuchten, M. Th. (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44, 892–898, doi:10.2136/sssaj1980.03615995004400050002x.

Examples

data(sim_wrc_hcc)

if(!requireNamespace("lattice", quietly = TRUE)){
  install.packages("lattice", repos = "https://cloud.r-project.org/")
  requireNamespace("lattice", quietly = TRUE)
}

lattice::xyplot(wc ~ head|id, type = "l", sim_wrc_hcc, 
  scales = list( x = list(log=TRUE)))
lattice::xyplot(hc ~ head|id, type = "l", sim_wrc_hcc, 
  scales = list( log = TRUE))

Common S3 Methods for Class fit_wrc_hcc

Description

This page documents the methods coef, summary, print, plot and lines for the class fit_wrc_hcc.

Usage

## S3 method for class 'fit_wrc_hcc'
coef(object, what = c("all", "nonlinear", "linear"),
    subset = NULL, residual_se = FALSE, se = FALSE, gof = FALSE, lc = FALSE,
    e0 = FALSE, bound = lc, ...)

## S3 method for class 'fit_wrc_hcc'
summary(object, what = c("all", "nonlinear", "linear"),
    subset = NULL, gof = TRUE, lc = TRUE, ...)

## S3 method for class 'fit_wrc_hcc'
print(x, ...)

## S3 method for class 'fit_wrc_hcc'
plot(x, what = c("wrc", "hcc"), y = NULL,
    subset = NULL, ylim_wc = NULL, ylim_hc = NULL,
    head_saturation = 0.01,
    beside = identical(sum(par("mfrow")), 2L), pch = 1, col_points = "black",
    col_line_x = "blue", lty_x = "solid",
    col_line_y = "orange", lty_y = "dashed",
    xlab_wc = "head [m]", ylab_wc = "water content [-]",
    xlab_hc = "head [m]", ylab_hc = "hyd. conductivity [m/d]",
    draw_legend = TRUE, draw_parameter = FALSE, cex_legend = 0.7, ...)

## S3 method for class 'fit_wrc_hcc'
lines(x, what = c("wrc", "hcc"), id = 1,
    head_saturation = 0.01, ...)

Arguments

Details

Residual standard errors, standard errors of the nonlinear parameters and confidence intervals based on the asymptotic normal distribution are computed only for mpd and ml estimates, see
soilhypfitIntro, control_fit_wrc_hcc and vcov.

The plot method for class fit_wrc_hcc displays for each sample the measurements of the water retention curve and/or the hydraulic conductivity function, along with the fitted model curve(s). Optionally, the curves of a second model fit (specified by the argument y) can be plotted for each sample.

The lines method adds the curve of a fitted model to an existing plot.

Value

The method coef returns a dataframe with the estimated parameters (and optionally standard errors), optionally the value of the objective function along with convergence code and/or information on the characteristic evaporative length.

The method summary generates a list (of class summary.fit_wrc_hcc) with the following components:

data

a named integer vector with the total number of samples (nsamp) and the number of samples with only water retention (nwrc), only hydraulic conductivity (nhcc) and both type of measurements (nwrchcc).

control

a list with a subset (settings, nloptr, sce, approximation_alpha_k0, param_bound, param_tf) of the components of object[["control"]], see fit_wrc_hcc and control_fit_wrc_hcc.

result

a dataframe with the estimated parameters and optionally the residual sum of squares along with convergence code and/or information on the characteristic evaporative length.

call

the call component of object.

Note that only a print method is available for class summary.fit_wrc_hcc.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

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;

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

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

data(sim_wrc_hcc)

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

# estimate parameters for 3 samples by unconstrained, global optimisation
# algorithm NLOPT_GN_MLSL
# sample 1: use only conductivity data
# sample 2: use only water content data
# sample 3: use both types of data
rfit_uglob <- fit_wrc_hcc(
  wrc_formula = wc ~ head | id, hcc_formula = hc ~ head | id,
  wrc_subset = id != 1, hcc_subset = id != 2,
  data = sim_wrc_hcc, fit_param = default_fit_param(tau = TRUE),
  control = control_fit_wrc_hcc(param_bound = param_boundf(
      alpha = c(0.00001, 50), n = c(1.0001, 7), tau = c(-1, 5)
    ), pcmp = control_pcmp(ncores = ncpu)))
print(rfit_uglob)
summary(rfit_uglob)
coef(rfit_uglob, what = "nonlinear")
coef(rfit_uglob, what = "linear", gof = TRUE)
coef(vcov(rfit_uglob), status = TRUE, se = FALSE)
op <- par(mfrow = c(3, 2))
plot(rfit_uglob)
on.exit(par(op))

Physical properties of Swiss forest soils

Description

The data give basic physical properties, water content and hydraulic conductivity measurements (at given capillary pressure head) for 128 soil layers (horizons) measured at 23 forest sites in Switzerland.

Usage

data(swissforestsoils)

Format

A data frame with 1373 observations on the following 21 variables.

profile_id

a factor with short labels for the 23 sites.

profile

a factor with with long labels for the 23 sites.

longitude

a numeric vector with the latitude of the site in degree.

latitude

a numeric vector with the latitude of the site in degree.

layer_id

a factor with labels for the 128 soil layer.

layer_ub, layer_lb

numeric vectors with the upper and lower depth (unit cm) of the soil layer for the measurements of the basic physical properties (particle_density, ..., ksat).

particle_density

a numeric vector with the density of the solid soil material (unit \mathrm{g} \,\mathrm{cm^{-3}}).

bulk_density

a numeric vector with soil (bulk) density (unit \mathrm{g} \,\mathrm{cm^{-3}}).

porosity

a numeric vector with the soil porosity (unit volume percentage).

clay

a numeric vector with the clay content (unit mass percentage).

silt

a numeric vector with the silt content (unit mass percentage).

sand

a numeric vector with the sand content (unit mass percentage).

ksat

a numeric vector with the saturated hydraulic conductivity (unit \mathrm{m} \,\mathrm{d^{-1}}).

head

a numeric vector with capillary pressure head at which theta (water retention curve) and/or ku (hydraulic conductivity function) were measured (unit m).

layer_ub_theta, layer_lb_theta

numeric vectors with the upper and lower depth (unit cm) of the soil layer for which the water retention curve was measured.

theta

a numeric vector with volumetric water content measurements (dimensionless) of the water retention curve.

layer_ub_ku, layer_lb_ku

a numeric vector with the upper and lower depth (unit cm) of the soil layer for which the water retention curve was measured.

ku

a numeric vector with hydraulic conductivity measurements (unit \mathrm{m} \,\mathrm{d^{-1}}) of the hydraulic conductivity function.

Details

clay, silt and sand refer to soil particles with diameter less than 2, between 2 and 50 and larger than 50 \mum.

Source

Richard, F. & Lüscher, P. 1978 – 1987. Physikalische Eigenschaften von Böden der Schweiz. Lokalformen Bände 1 – 4. Eidgenössische Anstalt für das forstliche Versuchswesen, Birmensdorf.

Examples

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

# estimate parameters using all samples (samples with water retention,
# hydraulic conductivity, or with both type of measurements)

data(swissforestsoils)

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

# unconstrained estimation (global optimisation algorithm NLOPT_GN_MLSL)
r_uglob <- fit_wrc_hcc(
  wrc_formula = theta ~ head | layer_id,
  hcc_formula = ku ~ head | layer_id,
  data = swissforestsoils,
  control = control_fit_wrc_hcc(
    settings = "uglobal", pcmp = control_pcmp(ncores = ncpu)))
summary(r_uglob)
coef(r_uglob)

Utility functions for package soilhypfit

Description

This page documents the functions convergence_message, extract_error_messages, select_failed_fits and check_param_boundf.

Usage

convergence_message(x, sce = FALSE)

extract_error_messages(object, start = 1, stop = 80, prt = TRUE)

select_failed_fits(object)

check_param_boundf(y, compare_thetar_thetas = TRUE)

Arguments

Details

The function convergence_message prints a message that explains the convergence codes, for nloptr, see NLopt return values. The function extract_error_messages extract the error messages of estimations that failed and optionally prints sub-strings of them.
The function select_failed_fits returns the ids of the soil samples for which parameter estimation failed. The function check_param_boundf checks the validity and consistecy of bounds of box constraints of model parameters.

Value

The function convergence_message and extract_error_messages return invisibly the convergence code or the error messages.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Duan, Q., Sorooshian, S., and Gupta, V. K. (1994) Optimal use of the SCE-UA global optimisation method for calibrating watershed models, Journal of Hydrology 158, 265–284, doi:10.1016/0022-1694(94)90057-4.

Johnson, S.G. The NLopt nonlinear-optimisation package. https://github.com/stevengj/nlopt.

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

convergence_message(3)
convergence_message(2, sce = TRUE)

vcov Method for Class fit_wrc_hcc

Description

This page documents the method vcov for the class fit_wrc_hcc and its coef method. vcov extracts the covariance matrices of the nonlinear parameters \widehat{\boldsymbol{\nu}} estimated by maximum likelihood or maximum posterior density.

Usage

## S3 method for class 'fit_wrc_hcc'
vcov(object, subset = NULL, grad_eps,
    bound_eps = sqrt(.Machine$double.eps), ...)

## S3 method for class 'vcov_fit_wrc_hcc'
coef(object, se = TRUE, correlation = se,
    status = FALSE, ...)

Arguments

Details

The function vcov extracts (co-)variances of the nonlinear parameters from the inverse Hessian matrix of the objective function at the solution \widehat{\boldsymbol{\nu}} for mpd and ml estimates, see soilhypfitIntro and Stewart and Sørensen (1981).

vcov checks whether the gradient at the solution is approximately equal to zero and issues a warning if this is not the case. This is controlled by the argument grad_eps which is the tolerable largest modulus of the scaled gradient (= gradient divided by the absolute value of objective function) at the solution. The function control_fit_wrc_hcc selects a default value for grad_eps in the dependence of the chosen optimisation approach (argument settings of control_fit_wrc_hcc).

vcov sets covariances equal to NA if the parameter estimates differ less than bound_eps from the boundaries of the parameter space as defined by param_boundf.

Value

The method vcov returns an object of of class vcov_fit_wrc_hcc, which is a list of covariance matrices of the estimated nonlinear parameters for the soil samples. The attribute status of the matrices qualifies the covariances.

The coef method for class vcov_fit_wrc_hcc extracts the entries of the covariances matrices, optionally computes standard errors and correlation coefficients and returns the results in a dataframe.

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Stewart, W.E. and Sørensen, J.P. (1981) Bayesian estimation of common parameters from multiresponse data with missing observations. Technometrics, 23, 131–141,
doi:10.1080/00401706.1981.10486255.

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;

prfloglik_sample for profile loglikelihood computations;

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

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

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

data(sim_wrc_hcc)

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

# estimate parameters for 3 samples by unconstrained, global optimisation
# algorithm NLOPT_GN_MLSL
# sample 1: use only conductivity data
# sample 2: use only water content data
# sample 3: use both types of data
rfit_uglob <- fit_wrc_hcc(
  wrc_formula = wc ~ head | id,
  hcc_formula = hc ~ head | id,
  wrc_subset = id != 1,
  hcc_subset = id != 2,
  data = sim_wrc_hcc,
  control = control_fit_wrc_hcc(pcmp = control_pcmp(ncores = ncpu)))
print(rfit_uglob)
summary(rfit_uglob)
coef(rfit_uglob, what = "nonlinear")
coef(rfit_uglob, what = "linear", gof = TRUE)
coef(vcov(rfit_uglob), status = TRUE, se = FALSE)
op <- par(mfrow = c(3, 2))
plot(rfit_uglob)
on.exit(par(op))

Models for Soil Water Retention Curves

Description

The functions sat_model and wc_model compute, for given capillary pressure head h, the volumetric water saturation S(h) and the volumetric water content \theta(h), respectively, of a soil by parametrical models.

Usage

sat_model(h, nlp, precBits = NULL, wrc_model = "vg")

wc_model(h, nlp, lp, precBits = NULL, wrc_model = "vg")

Arguments

Details

The functions sat_model and wc_model currently model soil water retention curves only by the simplified form of the model by Van Genuchten (1980) with the restriction m = 1 - \frac{1}{n}, i.e.

S_\mathrm{VG}(h; \boldsymbol{\nu}) = ( 1 + (\alpha \ h)^n)^\frac{1-n}{n}

\theta_\mathrm{VG}(h; \boldsymbol{\mu}, \boldsymbol{\nu}) = \theta_r + (\theta_s - \theta_r) \, S_\mathrm{VG}(h; \boldsymbol{\nu})

where \boldsymbol{\mu}^\mathrm{T} = (\theta_r, \theta_s) are the residual and saturated water content (0 \le \theta_r \le \theta_s \le 1), respectively, and \boldsymbol{\nu}^\mathrm{T} = (\alpha, n) are the inverse air entry pressure (\alpha > 0) and the shape parameter (n > 1).

Note that \boldsymbol{\mu}^\mathrm{} and \boldsymbol{\nu}^\mathrm{} are passed to the functions by the arguments lp and nlp, respectively.

Value

A numeric vector with values of volumetric water saturation (sat_model) or water content (wc_model).

Author(s)

Andreas Papritz papritz@retired.ethz.ch.

References

Van Genuchten, M. Th. (1980) A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal, 44, 892–898, doi:10.2136/sssaj1980.03615995004400050002x.

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;

evaporative-length for physically constraining parameter estimates of soil hydraulic material functions.

Examples

## define capillary pressure head (unit meters)
h <- c(0.01, 0.1, 0.2, 0.3, 0.5, 1., 2., 5.,10.)

## compute water saturation and water content
sat <- sat_model(h, nlp = c(alpha = 1.5, n = 2))
theta <- wc_model(
  h ,nlp = c(alpha = 1.5, n = 2), lp = c(thetar = 0.1, thetas = 0.5))

## display water retention curve
op <- par(mfrow = c(1, 2))
plot(sat ~ h, log = "x", type = "l")
plot(theta ~ h, log = "x", type = "l")
on.exit(par(op))