soilhypfit_package: The soilhypfit Package

In soilhypfit: Modelling of Soil Water Retention and Hydraulic Conductivity Data

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, \mbvec\thetaT=(\theta_1, \theta_2, ..., \theta_n_\theta), and/or soil hydraulic conductivity, \mbvecKT=(K_1, K_2, ..., K_n_K), both measured at given capillary pressure head, \mbvechT=(h_1, h_2, ...).

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

θ_i = θ(h_i; μ, ν) + ε_{θ,i}, i = 1, 2, ..., n_θ,

log(K_j) = log(K(h_j; μ, ν)) + ε_{K,j}, j = 1, 2, ..., n_K,

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

Let

Q_θ(μ, ν; θ, h) = (θ - θ(h; μ, ν))^T W_θ (θ - θ(h; μ, ν)),

and

Q_K(μ, ν; K, h) = (log(K) - log(K(h; μ, ν)))^T W_K (log(K) - log(K(h; μ, ν))),

denote the (possibly weighted) residual sums of squares between measurements and modelled values. θ(h; μ, ν)) and K(h; μ, ν) are vectors with modelled values of water content and hydraulic conductivity, and W_θ and W_K are optional diagonal matrices with weights w_{θ,i} and w_{K,j}, respectively. The weights are products of case weights w^\prime_{θ,i} and w^\prime_{K,j} and variable weights w_θ, w_K, hence w_{θ,i} = w_θ * w^\prime_{θ,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(μ, ν; θ, K, h) = κ_θ / 2 Q_θ(μ, ν; θ, h) + κ_K / 2 Q_K(μ, ν; K, h),

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

For wls the objective function is equal to

Q(μ, ν; θ, K, h) = Q_θ(μ, ν; θ, h) + Q_K(μ, ν; K, h).

If either water content or conductivity data are not available, then the respective terms are omitted from Q(μ, ν; θ, K, h).

The function fit_wrc_hcc does not attempt to estimate the parameters by minimising
Q(μ, ν; θ, K, h) directly with respect to \mbvec\mu and \mbvec\nu. Rather, it exploits the fact that for given nonlinear parameters \mbvec\nu, the conditionally linear parameters μ^T = (θ_r, θ_s, log(K_0)) can be estimated straightforwardly by minimising the conditional residual sums of squares

Q*_θ(θ_r, θ_s; θ, h, ν) = (θ - [1, S(h;ν)] [θ_r, θ_s - θ_r]^T])^T W_θ (θ - [1, S(h;ν)] [θ_r, θ_s - θ_r]^T])

with respect to θ_r and θ_s-θ_r and/or

Q*_K(K_0; K, h, ν) = (log(K) - log(K_0 k(h; ν)))^T W_K (log(K) - log(K_0 k(h; ν))),

with respect to log(K_0), where \mbvec1 is a vector of ones, S(h;ν)^T = (S(h_1;ν), S(h_2;ν), ..., S(h_{n_θ};ν)) and k(h;ν)^T = (k(h_1;ν), k(h_2;ν), ..., k(h_{n_K};ν)) are vectors of modelled water saturation and modelled relative conductivity values, θ_r and θ_s are the residual and saturated water content, and K_0 is the saturated hydraulic conductivity.

Unconstrained conditional estimates, say hatθ_r(ν), hatθ_s(ν) - hatθ_r(ν) and hatlog(K_0)(ν) 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 ≤ θ_r ≤ θ_s ≤ 1.

Let hatμ(ν)^T = (hatθ_r(ν), hatθ_s(ν), hatlog(K_0)(ν)) be the conditional estimates of the linear parameters obtained by minimising Q*_θ(θ_r, θ_s; θ, h, ν), and Q*_K(K_0; K, h, ν), respectively. fit_wrc_hcc then estimates the nonlinear parameters by minimising Q(hatμ(ν), ν; θ, K, h) with respect to \mbvec\nu 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)

hatσ^2_θ = Q_θ(hatμ(hatν), hatν; θ, h) / κ_θ,

and

hatσ^2_K= Q_K(hatμ(hatν), hatν; K, h) / κ_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(hatμ(ν), ν; θ, K, h) at the solution hatν (Stewart and Sørensen, 1981), i.e.

Cov[hatν, hatν^T] ≈ A^{-1},

where

[A]_{k,l} = \partial^2 / (\partialν_k, \partialν_l) Q(hatμ(ν), ν; θ, K, h)|_(ν = hatν).

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 - 1/n, see wc_model and hc_model. This model has the following parameters:

• \mbvec\mu

  T= (\theta_r, \theta_s, K_0) (see above) and

• \mbvec\nu

  T= (\alpha, n, \tau) where α is the inverse air entry pressure, n the shape and τ 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 \mbvec\mu, nonlinear optimisation algorithms with box constraints for estimating \mbvec\nu) 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_c of stage-I evaporation from a soil (Lehmann et al., 2008). L_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 α 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_t of L_c for given n and tortuosity parameter τ, see evaporative-length. fit_wrc_hcc allows to constrain estimates of the nonlinear parameters \mbvec\nu by defining permissible lower and upper bounds for the ratio L_c/L_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 \mbvec\nu, fit_wrc_hcc minimises Q(hatμ(ν), ν; θ, K, 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_c/L_t
   (settings = "uglobal" or settings = "sce"),

2. global optimisation with inequality constraints for the ratio L_c/L_t
   (settings = "cglobal"),

3. local optimisation without constraints for the ratio L_c/L_t
   (settings = "ulocal"),

4. local optimisation with inequality constraints for the ratio L_c/L_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 \mbvec\nu are indispensable for successful estimation. fit_wrc_hcc allows to compute initial values of α and n for the Van Genuchten model from water retention data by the following procedure:

1. Smooth the water retention data, (θ_i, y_i = log(h_i)), i=1,2, ... n_θ, 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 hatm, of S* = (1 + 1/m)^{-m}. One obtains the right-hand side of this equation by solving d^2/dy^2[S_VG(exp(y); ν)] = 0 for y and plugging the result into the expression for S_VG(exp(y); ν), see wc_model.

4. Compute hatn = 1 /(1-hatm) and hatα = 1 / exp(y*) * (1/hatm)^{1-hatm}. The second expression is again a result of solving d^2/dy^2[S_VG(exp(y); ν)] = 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_c/L_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 [e-mail].

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.

soilhypfit documentation built on Aug. 31, 2022, 5:09 p.m.