Computers & Geosciences 26 (2000) 319±327

RETMCL: Incorporating maximum-likelihood estimation
principles in the RETC soil hydraulic parameter estimationp
code
K.J. Hollenbeck a, b, J. SÏimunek c,*, M.Th van Genuchten caDepartment of Hydrodynamics and Water Resources, Technical University of Denmark, Lyngby, Denmarkb
Danish Hydraulic Institute, Agern AlleÂ 5, DK-2970, Hùrsholm, Denmarkc
US Salinity Laboratory, USDA-ARS, Riverside, CA, USA
Received 3 May 1999; received in revised form 13 August 1999; accepted 13 August 1999

Abstract
RETC is a public domain computer code for estimating parameters of the water retention curve and hydraulicconductivity functions of unsaturated soils. RETC was developed at the US Salinity Laboratory and is now usedworld-wide with thousands of copies distributed. Evaluation of the ®nal estimation results in the code has been
improved to yield a new version, RETCML, based on the maximum-likelihood theory for the special case ofweighted least-squares estimators. This paper ®rst explains the theory of maximum-likelihood and introduces theprinciples of model adequacy and parameter uncertainty on a formal basis. Next, this paper presents a user guidefor the code. RETCML is also free and has been programmed to be almost fully compatible with the original
RETC input ®les, thus making it easy to re-analyze previous data. The output of RETCML includes a thoroughevaluation of estimation results. # 2000 Elsevier Science Ltd. All rights reserved.
Keywords: Maximum-likelihood; Estimation; Soil hydraulic properties

1. Introduction
The RETC (RETention Curve) code is a widely used
computer program developed at the US Salinity Laboratory for estimating parameters of the retention
curve and hydraulic conductivity functions of unsaturated soils (van Genuchten et al., 1991). While theCodeavailableat
http://www.iamg.org/CGEditor/
index.htm and http://www.ussl.ars.usda.gov/models/models.html.
* Corresponding author. Tel.: +1-909-369-4865.
E-mail addresses: kjh@dhi.dk (K.J. Hollenbeck), jsimunek@ussr.ars.usda.gov (J. SÏimunek).
pretention curve (often also called the soil water characteristic curve) characterizes the energy status of the soil
water, the unsaturated hydraulic conductivity function
describes the ability of the porous medium to conduct
water. The RETC program uses the parametric modelsof Brooks and Corey (1964) and van Genuchten(1980a) to represent the soil water retention curve, andthe theoretical pore-size distribution models of Mualem (1976) and Burdine (1953) to either predict theunsaturated hydraulic conductivity function fromobserved soil water retention data or to use those datain the ®tting procedure also. The van Genuchten retention function (van Genuchten, 1980a) has been very
popular in the ®elds of soil physics and vadose zone

0098-3004/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved.
PII: S 0 0 9 8 - 3 0 0 4 ( 9 9 ) 0 0 1 3 8 - 7

320

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

hydrology. The suggested S-shape function and its
combination with the pore distribution model of Mualem (1976) is used not only for the description of thesoil water characteristic curve (retention curve) but
also for functions describing oil retention and transport properties.
The retention and hydraulic conductivity functions
constitute properties needed for simulating ¯ow andsolute transport in the unsaturated zone. Computermodels using analytical descriptions of soil hydraulicproperties are now routinely used in both research andmanagement to predict the movement of water andchemicals in the vadose zone between the soil surface
and the groundwater table. Interest in the vadose zone
has dramatically increased in recent years because of
growing evidence that the quality of the subsurface environment is being adversely aected by industrial, municipal, and agricultural activities (van Genuchten etal., 1991).
The public domain RETC code, and its earlier version SOHYP (SOil HYdraulic Properties) (van Genuchten, 1978), have been freely distributed by several
agencies, including the US Salinity Laboratory (USSL)of the US Department of Agriculture, Agriculture
Research Service (USDA, ARS), the US Environmental Protection Agency (EPA), the International
Groundwater Modeling Center (IGWMC), and others.
Hundreds of copies of RETC have been distributed by
USSL and IGWMC, thousands by EPA, while
uncounted numbers of copies have been downloaded
from internet sites of these and other agencies. Becauseof the wide distribution of RETC (and its descendant/
successor programs), the implementation of the estimation problem in RETC, including the statisticalanalysis of estimated parameters, was carried over tonumerous other programs, including CFITM (vanGenuchten, 1980b), CFITIM (van Genuchten, 1981),
ONESTEP (Kool et al., 1985), SFIT (Kool et al.,1987), MULTISTEP (van Dam et al., 1994), 3DADE(Leij and Bradford, 1994), CXTFIT (Toride et al.,1995), HYDRUS-1D (SÏimunek et al., 1998) andHYDRUS-2D (SÏimunek et al., 1996), and probably
many others unknown to the authors. While these programs all correctly implement the estimation problem,the statistical analysis concerning the ®nal estimatedparameters is limited because of its imperfect implementation.
We here present an improved version of RETC,
called RETCML (RETention Curve Ð MaximumLikelihood), for which the uncertainty analysis isbased on maximum-likelihood theory for the special
case of weighted least-squares estimators. In
RETCML, two theorems for evaluating the estimation
results are exploited. The ®rst theorem allows an evaluation of whether or not a selected parametric model is``adequate''. Adequacy here signi®es that the model

predictions come so close to the data that any remaining discrepancy is no larger than measurement error in
the underlying experiment. Given adequacy, the second
theorem provides a basis for evaluating uncertainty inthe estimated parameters as a function of measurement
uncertainty. To bene®t from RETCML's features, theuser thus has to specify the measurement uncertainty
in the input ®le. Users of all other above-mentioned
programs may also want to implement the suggestedchanges in their codes in order to be able to morethoroughly evaluate the model adequacy and con®dence intervals of optimized parameters in their particular applications.

2. About estimation
Parameter estimation refers to the problem where avariable y is observed as a function of a control variable x, and one has a corresponding model with parameters grouped in a vector p. The optimal parameters
are then generally (also in RETC and RETCML)
found by minimizing a sum-of-squares (SSQ) objectivefunction:
SSQ p 

N 
X
1
i1

syi

yi, meas x i  ÿ yi, model x i ;p21

where the subscripts ``meas'' and ``model'' refer tomeasured and model predicted values of y, and N is
the number of measurements. A ``model'' here means aparticular parametric expression for the data, e.g., thevan Genuchten retention expression"
c y;p y ÿ yrys ÿ yr#1=nn= 1ÿn
ÿ1
a

2

Fig. 1. Some example retention data and corresponding best®t van Genuchten model.

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

with parameters p=(yr, ys, n, a ) or any subset thereof
if some of these are known by other means. An
example of retention data and the corresponding best®t van Genuchten model is shown in Fig. 1. As thispaper focuses on the estimation problem, the reader isreferred to the literature for additional examples ofhydraulic functions and properties (van Genuchten,
1980a).
A key feature of the above SSQ estimator are the i= 1, . . ., N weights sÿ1
yi : As the notation implies, theweights must be set to re¯ect measurement error information about the data:
1
1

:
syi
measurement error standard deviation i3

It is intuitive that information about the quality of thedata should determine the quality of the parameter
estimates, and that this information hence should becontained in the SSQ estimator. This intuitive principle
has a formal mathematical representation. Withweights set as given by Eq. (3), the weighted leastsquares (WLS) estimator Eq. (1) is in fact a maximumlikelihood (ML) estimator. The equivalence between
the WLS and the more general ML estimator strictly
holds only when all normalized residuals ri  yi ,meas x i  ÿ yi , model x i ; p=syi are (i) independent and(ii) normally distributed. Condition (ii) is not reallyrestrictive because the estimator is robust as long as
the error distributions are symmetric and free of a bias
and extreme outliers (Press et al., 1992, Chap. 15.6).
While condition (i) is more dicult to meet in practice,
the case of residuals with unknown correlation is in
general dicult to treat, even by more advanced theories. Since this paper is not intended to discuss all
aspects of ML estimation, the interested reader isreferred to the literature (e.g., Bard, 1974; Beck and
Arnold, 1977). Most importantly for practical applications, the error introduced by ignoring some residual
correlation is much less severe than the error introduced by not setting the weights properly.
Once the weights are set as in Eq. (3), two theorems
for evaluating estimation results can be exploited. The®rst theorem allows an evaluation of whether or not aselected parametric model is adequate. For adequate
models, the second theorem provides a basis for evaluating uncertainty in the estimated parameters. The twotheorems and their implementation in RETCML will
be presented in the following. In addition, the dierences between RETC and RETCML will be pointed
out.
2.1. Model adequacy
An adequate model is one that explains the data tosuch a degree that it is plausible that remaining discre-

321

pancies are due only to measurement error noise. In
other words for the best ®t, the absolute values of thedeviations di  yi , meas x i  ÿ yi , model x i ; p shouldbe about as large as the measurement error standard
deviations syi , yielding normalized residuals each withj ri j 11, thus r2i 11: Furthermore, the SSQ estimatoris a sum of N squared residuals r2i , so if these are all
independent, it can be expected that at best, SSQ 1 N.
In addition, it is known that a sum of squared randomvariables that each follow a normal distribution has achi-square distribution.
The preceding intuitive argument may elucidate theprinciple behind the formal statistical test for model
adequacy. This test is based on the minimum SSQ
value, SSQmin, attained for the optimal parameters
popt, which is distributed as (Press et al., 1992, Chap.
15.1):
SSQmin 0w2NÿP :4

The number of degrees-of-freedom is NÿP and not N,
because there is a loss of one degree-of-freedom foreach of the P parameters.
The expected value (the mean) of a w 2 function is its
number of degrees-of-freedom (Abramowitz and Stegun, 1965 , Eq. 26.4.33):
Ew2NÿP   N ÿ P:

5Hence, from Eq. (4), it can be expected that
SSQmin=NÿP. As any estimation problem is bestde®ned when the number of data signi®cantly exceedsthat of the number of parameters, a good rule ofthumb is that for an adequate model, SSQmin shouldbe about the same as the number of data Ð just as
was expected from the introductory intuitive argument.
If the model is inadequate, however (e.g., the situ-

Fig. 2. ``Fitting'' S-shaped van Genuchten retention curve tosynthetic data representing straight line. R 2 value is high evenfor this inadequate model, whereas padeq correctly reveals mis®t.

322

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

ation shown in Fig. 2), SSQmin will be signi®cantly larger than NÿP. This is so because the model-data deviations di will (by absolute value) on average be larger
than the measurement errors syi ; hence, the squared residuals will on average be larger than 1.
For any observed value of SSQmin, the probabilityof model adequacy is:padeq  1 ÿ Q SSQmin , N ÿ P 

6

where Q() is the w 2 cumulative density function.
The probability of model adequacy, padeq, in practical applications is relatively sensitive to outliers in thedata. Even a single outlier, giving one large squared residual r2outlier , can result in a high SSQmin. Hence, a rare
relatively low padeq is not necessarily reason for concern.
More often, the reason for a low padeq is the assignment of too high weights sÿ1
yi , as the SSQ grows withthe square of the weights. Too high weights sÿ1
yi resultif the user's speci®cation for the measurement errors
syi is too optimistic (too small). Conversely, if theseerrors are speci®ed too pessimistically (too large),
SSQmin will be very small, resulting in an often unrealistic probability of model adequacy essentially equal
to 1.
Also the original RETC allows for speci®cation of
the weights for the data. On p. 35 in the manual (vanGenuchten et al., 1991), the recommendation of setting
the weights as inverses of observation errors is given,
but its importance is not particularly emphasized.
Accordingly, the probability of model adequacy is notcomputed in RETC, but only the correlation coecientR 2 between model predictions and data. However, theR 2 value does not follow any distribution, so one cannot evaluate quantitatively how good a particular value
is. In general, the R 2 values found with RETC are allvery close to 1, even for the clearly inadequate model.
The advantage of examining model adequacy rather
than a correlation coecient is illustrated for a ``nonsense'' example (Fig. 2). The synthetic data used for
estimating the van Genuchten retention parameters
were set to describe a straight line, which is not aretention curve one can observe in nature. Nevertheless, the R 2 value output by RETC is 0.98, close to theperfect value of 1. In contrast, the probability ofmodel adequacy, padeq, as computed in RETCML,
unequivocally reveals the poor ®t. In general, ``good''®ts are correctly indicated by both R 2 and padeq (forthe ®t in Fig. 1, R 2=0.99710 and padeq > 0.99999).
Apparently, the statistic R 2 thus has little power fordiscriminating between good and poor ®ts. (Note thatmeasurement error standard deviations were taken to
be as large as the crosshairs on the data in Fig. 2 whencalculating padeq.)

2.2. Parameter con®dence regions
A con®dence region for parameter estimates isde®ned via a maximum allowable change in the SSQ
from its minimum value. The relevant theorem has
been implemented in RETCML, but only incorrectly
so in the original RETC.
From theory (Press et al., 1992, Chap. 15.6), it is
known that the increment DSSQ(p)=SSQ(p)ÿSSQmin
follows a chi-square distribution with P degrees-offreedom:
DSSQ p0w2P :7

Thus, the appropriate value for DSSQconf for a desired
level of con®dence pconf is:
DSSQconf  Qÿ1 pconf , P 

8

where Q ÿ1( pconf, P ) is the inverse w 2 cumulative density function.
In practice it would be dicult to compute the SSQ
for ``all'' parameter vectors p and determine for each
one whether DSSQ(p) is either larger or smaller than
the allowed DSSQconf value. (Although contouring the
SSQ is often a method for analyzing how well-de®ned
estimation problems are (e.g., Toorman et al., 1992).)
Fortunately, an approximate method exists based on
the Cramer±Rao theorem (Schweppe, 1973, Chap. 12).
This theorem gives a lower bound for the parameter covariance matrix S p as:
S p rHÿ19

where H is the Hessian matrix. H can be computedfrom the sensitivity of the model to its parameters, i.e.,from the partial derivatives of the objective function(Press et al., 1992, Chap. 15.5), as
Hkl 



N
X
1 @ yi p @ yi pi1s2y, i@ pk@ plk, l  1, . . . , P:

10

With the parameter covariance matrix known, a geometric description of the con®dence region is given
through:p ÿ popt  0 S ÿ1p p ÿ popt   DSSQconf :11

Eq. (11) describes an ellipsoid in the parameter space,and all parameter vectors p within the ellipsoid arepart of the con®dence region.
Eqs. (10) and (9) suggest that smaller weights imply
a larger (less desirable) con®dence region. This makes
sense because smaller weights express larger measurement errors. Now recall that smaller weights also
imply a higher (more desirable) probability of model
adequacy. Hence, the magnitudes of the weights

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

impact the two measures of the quality of the estimation result in opposite directions. Inadequate
models can have attractively small parameter con®dence regions. The user should not be tempted to set
the weights too large to achieve a small con®dence
region, but at the expense of an unrealistically small
probability of model adequacy.
A multivariate con®dence region in P dimensions isdicult to visualize, but pseudo-univariate ``parameterstandard errors'' can be taken from the diagonal of theparameter covariance matrix as:p
sp, i  S p, ii i  1, . . . , P:
12
A conservative way to ®nd the actual con®dence interval for a parameter is to ®nd the projections of thecon®dence ellipsoid on the parameter axes. The halflength of the symmetric con®dence interval for the i-th
parameter is (Press et al., 1992, Eq. 15.6.4):p
cp, i  Qÿ1 pconf , P sp, i :13An example of a con®dence ellipse and the projectedcon®dence intervals computed using Eq. (13) is shownin Fig. 3. This schematic example is for two genericparameters here called ``p1'' and ``p2''.
It is important to be aware of the approximate
nature of the Cramer±Rao theorem when choosing acon®dence level for the parameters. Only for linear
models is the con®dence region an ellipsoid, with its
size proportional to the norm of the covariance matrix.
The uncertainty region for nonlinear models is often
odd-shaped, and larger than the Cramer±Rao ellipsoid.
The discrepancy grows the higher the con®dence level,i.e., the larger the region. For example, while none ofthe van Genuchten parameters can be negative, for a

Fig. 3. Schematic representation of two-dimensional con®dence ellipse with projections delineating con®dence intervals
for parameters ``p1'' and ``p2'' around optimal (estimated) parameters.

323

suciently large desiredcon®dencep
 level (a suciently
large multiplierQÿ1 pconf , P ), the Cramer±Rao
ellipsoid will grow into the negative parts of parameter
space. Hence, it is not advisable to ``stretch'' the theorem by computing con®dence limits at too high a level.
With this in mind, RETCML is coded to output 90%
rather than 95% (RETC) con®dence limits. Exact andCramer±Rao con®dence regions for the soil hydraulicdata are shown in Hollenbeck and Jensen (1998, Fig.
4) . The approximation can be judged admissible.
In RETC, univariate parameter con®dence halflengths are apparently computed incorrectly. The formula used in the code is:
cp, i  tÿ1 pconf , N ÿ P r
SSQminsp, i
NÿPi  1, . . . ,14

P
where t ÿ1( pconf, NÿP ) is the inverse cumulative density of the t distribution. With this formulation, theweight factors almost cancel out, because as seen from
Eqs. (9), (10), and (12), the squared weights are contained in the parameter standard error sp,i, whilefromp
Eq. (1), their reciprocals are contained in SSQmin :The net result is that dierent multiples of weights
result in the same computed standard errors. Considering the de®nitions more closely, the weights enter as
multipliers for the deviations, which as a product aresquared in SSQmin, while in the Hessian, the weights
themselves are squared. Hence, the two sums de®ning
SSQmin and Hÿ1
kl are not exactly reciprocals of eachother when not all weights are scaled by the same factor. The use of dierent weights within a data set (i.e.,
not a dierent multiple for all) will impact the estimates to a limited extent. For example, in RETC one
can assign larger weights for ``wet'' retention observations and smaller ones for some of the ``dry'' data,thus placing more emphasis on the ``wet'' data as compared to the ``dry'' data. The overall magnitude of theestimation error, however, will be determined only byNÿP. Recalling the discussion of model adequacy, theestimation errors as computed from RETC are always
about equal to that of an adequate model, irrespectiveof measurement uncertainty [cf Eq. (5)].
RETC makes use of t statistics in two instances.
Comparing RETC's Eq. (14) with RETCML's Eq.
(13), it can be seen that the scale factor for the projections of the parameter con®dence interval diers.
RETC's tÿ1
NÿP statistic is only correct for a univariate
estimation problem, but not for the generally encountered multivariate case. Furthermore, RETC's con®dence analysis output contains the t value for testing
the hypothesis H: pi=0. It is not obvious why one
should test whether a particular parameter is equal tozero (some, such as van Genuchten's n, are restricted

324

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

to be greater than zero or even one, depending upon
the selected soil-hydraulic parametric model). Hence,
the t values have not been carried over from RETC's
to RETCML's output.
As a closing remark for this Section, it should be
mentioned that maximum-likelihood theorems for evaluating estimation results are strictly valid only if thenumber of data clearly exceeds the number of parameters to be estimated. This restriction does not constrain RETCML relative to RETC, because as seen in
Eq. (14), also RETC makes reference to a maximumlikelihood theorem (the Cramer±Rao theorem). A``reasonably'' large number of data is required for solving any estimation problem in a meaningful way.
2.3. Assigning proper weights
Although Eq. (3) is a formally correct and complete
guideline on how to set weight factors in RETCML,
the user might encounter situations where moredetailed directions can be helpful. These are given inthis Section.
2.3.1. Errors in both x and y or in x only
According to Eq. (3), the SSQ estimator Eq. (1),allows for speci®cations of errors in y only, when y isa function of x. In practice, however, data may sometimes be available in the form x( y ), with errors known
for x. For example, in RETC and RETCML, theretention curve is ®tted as y(c ), not c(y ). A related
situation often encountered in practice is that instrument errors aect not only the measured variable y,but also the independent variable x.
To still use RETC and RETCML in the case oferrors in both y and x, or only in x, we suggest the following approach. One can compute a lumped total
error in y as:s
2
@ysyi, total  s2yi 15
sx i :@x i
The ®rst term, syi , is the error in y itself, and the second represents the change/error in y caused by a random change/error in x. For the unusual case of errorfree y data, the second term under the square root
guarantees a non-zero lumped error (otherwise theweight would be ``1/0=1''). The total error information syi, total then replaces syi in the objective function
Eq. (1).
Note that the scheme expressed in Eq. (15) is notpart of RETCML. The user has to compute the proper
weights separately and then make them part of the
input ®le for RETCML. For application of Eq. (15), it
is useful to assume a locally linearized functional
dependence of y on x, and to take ®nite-dierences

between consecutive data as approximations of the derivatives.
2.3.2. Log-transformed data
RETC's option to log-transform conductivity anddiusivity data has been retained in RETCML. The internal transformation of the input weights performedin RETC has, however, not been maintained. Rather,
the weights entering the SSQ are the input weights,
because the general philosophy with RETCML is thatthe user should assess his/her data quality and havefull control over the estimation process. The only
reason for using the log-transform feature should be tobetter ful®ll the normality assumption for the measurement errors (i.e., one should only log-transform if the
log-transformed error is closer to being normally distributed than the raw error). A wide spread in themagnitudes of data in itself is no reason to log-transform. There is also nothing wrong with specifying
weights of widely dierent magnitudes.
2.3.3. Error magnitude unknown
Often in practical applications, the measurement
error magnitude is unknown, and not even an educated guess can be made. RETCML can be used toperform ordinary-least-squares (OLS) estimation forthis case, essentially emulating RETC. Parameter
uncertainty is computed under the assumption of having an adequate model [by setting SSQmin=NÿP internally, cf Eq. (5)]. The output value of SSQmin is,
however, that for unit weights. This value can be useful for comparing dierent parametric models for thesame data in a relative sense (e.g., the lower SSQmin,the more appropriate the model). The probability ofmodel adequacy loses its meaning once SSQmin has
been scaled according to the quality of the ®t.
2.3.4. Mixed estimator
RETC allows for estimation from retention andunsaturated conductivity data simultaneously. In
essence, this is performed with a mixed sum-of-squares
estimator which uses a relative weight W1.
SSQ p  W1 SSQk p  SSQy p

16

where SSQk (p) is a conductivity (or diusivity) objective function and SSQy (p) analogously a retention
objective function, both de®ned as in Eq. (1).
In the preceding sections it was emphasized thateach datum should be given its own weight, re¯ecting
its reliability (the inverse of measurement error).
Hence, a factor W1 is formally not only super¯uous,
but can distort estimation results. Imagine the very
realistic case where each sub-SSQ in the estimator Eq.
(16) would lead to dierent optimal parameter vectors
p. By manipulating W1, one can shift the ®nal estimate

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

from one that is optimal for the conductivity data (W14 1) to one that is optimal for the retention data
(W1=0). The pragmatic suggestion in RETC's manual(p. 36) to set W1 within limits of 0.1 and 1 makes the
problem of arbitrary results less extreme, but it still
exists. No theory is available for setting W1, and as
shown here, the factor is not even needed.
Notwithstanding the above discussion, RETCML's
input also includes the factor W1. The main reason forretaining W1 was to achieve compatibility with the old
input ®le format. W1 is generally ignored in the computation, with a single exception. When maximum-likelihood estimation is impossible anyway, i.e., when OLS
estimation has been invoked, W1 can be used according to Eq. (16). In this fashion, the conductivity/diusivity data as a group can be weighed relative toretention data as a group. Still, the user should beaware of possibly distorted results, depending on howeach group determines the overall estimate. Therefore,
when using mixed estimators, it is recommended tocheck the output ®le for how each group contributed
to the optimal SSQ (see ``weighted sum of squares ofobserved versus ®tted values''). A W1 $1 in WLS analyses or when estimating from retention data only
results in a warning note in the output. For a more indepth discussion of the mixed-estimator problem, seeHollenbeck and Jensen (1998).

3. Using RETCML
3.1. Input
RETCML is compatible with RETC in terms of the
input ®le format (.in) and the control ®le format andname (still retc.ctl). Hence the user can easily reanalyze existing data. Likewise, the RETC manual(van Genuchten et al., 1991) is still generally valid for
RETCML. However, two cases exist where the input isinterpreted dierently.
One dierence concerns the way in which the OLS
estimation is carried out (see Section 2.3.3 for conditions where OLS is appropriate and how to interpretits results). In RETC, an OLS analysis could be
invoked by setting all weights to 1. This feature cannot
be carried over to RETCML because unit weights
could also, just by chance, re¯ect the true measurement
error. Hence, if an OLS analysis is desired with
RETCML, at least one weight has to be set to 0. A
zero weight is interpreted as missing, and because it
does not make sense to conduct a WLS estimation ifeven only a single weight is missing, one or more zeroweight(s) will invoke an OLS analysis. In this context,
the user ought to be careful to not accidentally introduce zero weights that can result from application of325

Eq. (15), i.e., when the gradient between a pair of data
points y(x ) is in®nite and sy=0. In this situation, theuser should remove one point of the pair (one line)from the input ®le. A data point with zero weightshould not aect the estimation results anyway,
because the apparent in®nite measurement error at thisdata point implies that it lacks any reliability.
The second dierence in input interpretation
between RETC and RETCML concerns the weight W1(position 16 in retc.ctl, weight of k/d datawith respect to water content data). As discussed in Section 2.3.4, this weight has no meaning
within the context of maximum-likelihood estimation
where each data point already has an appropriate
weight. Hence, unless OLS for two types of data isconducted, RETCML overrides any input, sets W1=1,
the neutral value, and prints a warning message. Notethat W1 in RETCML has not exactly the same eect
as in RETC, where in addition it is manipulated internally relative to data magnitude. This latter internalmanipulation is not transparent to the user and has
therefore been abandoned.
3.2. Output
RETCML preserves RETC's output of data andmodel predictions. The output for the estimation
results is, however, dierent. The individual items are
given below in the sequence they appear. Items in parentheses are not always part of the output, but dependon the particular estimation result.
1. Banner: Estimation results are announced with a
new banner printout. Unlike for RETC, the printout comes before the correlation matrix, which in
RETCML is interpreted as an estimation result,too.
2. (Warning) if at least one weight is missing (i.e.,
speci®ed as zero) and an OLS rather than a WLS
estimation has been performed (under the assumption of an adequate model; scaling the residuals as
discussed in Section 2.3.3).
3. (Warning) if there are few data (cut-o set to NÿP< 20), in which case the theorems used for interpreting estimation uncertainty do not strictly hold(cf Section 2.2).
4. (Note) if a relative weight for retention versus conductivity/diusivity SSQ's W1 $1 was used as input,
unless OLS is conducted for two data types. In allother cases, this factor is no longer supported as discussed in Sections 2.3.4 and 3.1.
5. Correlation matrix: Same as in RETC. This matrix
can help the user to identify ill-de®ned estimation
problems (for which the o-diagonal elements areclose to 21).
6. Probability of model adequacy: The value padeq as326

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

de®ned in Eq. (6). This quantity takes the place of
RETC's R 2 value for reasons explained in Section
2.1. Note that if OLS is conducted, the residuals are
scaled to give an adequate model [cf Eq. (5)]; hence,
the now meaningless padeq is not output.
7. (Warning) if padeq is very small (cut-o set to 10ÿ4).
Con®dence limits on parameters make little sense
for a very likely inadequate model, even though forthe sake of completeness, they still are presented.
The possible error leading to a very small padeq, i.e.,
too large (optimistically speci®ed) weights, ispointed out.
8. (Warning) if padeq is unrealistically large [cut-o set
to SSQmin/(NÿP ) < 0.01; cf Eq. (6) and Section2.1]. Computation of con®dence limits in this casebecomes numerically unstable [the Hessian has avery high condition number, so its inverse is illde®ned, cf Eq. (9)]. Con®dence limits likely become
incorrectly small. The possible error leading to avery large padeq, i.e., too small (pessimistically speci®ed) weights, is pointed out.
9. Parameter estimates and con®dence: Similar to
RETC's output, also including the pseudo-univariate parameter standard errors ``std.err.'' as de®ned
in Eq. (12). Con®dence limits are, however,obtained as the optimal parameter 2 the revised(ellipsoid projection) half-lengths as de®ned in Eq.
(13). The chosen con®dence limit is 90% as opposed
to RETC's 95% for reasons discussed in Section
2.2. The t-values are no longer given as output forreasons also discussed in that Section.
10. Observed and ®tted data: This part of the output isnearly identical to that of RETC, with the exception of the normalized residuals ri (labeledhvariable-nameiRES) replacing the deviations
di (labeled hvariablenameiDEV) when WLS
estimation has been performed. The normalized residuals ri should be close to 21 on average (cf Section 2.1). Larger (in absolute terms) values of theresiduals may indicate outliers in the data. For
OLS, normalized residuals are unde®ned for lackof a measurement error weight; hence, the deviations di are given as output, just as in RETC. A
note to that extent is printed out.

4. Practical information
The source code of RETCML is available under
http://www.iamg.org/CGEditor/index.htm. The installation ®le of the Windows version is available under
http://www.ussl.ars.usda.gov/models/models.html.
RETCML has been veri®ed by comparing estimation results with those obtained with MATLAB

code written for the same purpose, but using
MATLAB's own optimization routine constr. Only
small numerical errors due to dierent machine precisions were observed; hence, RETCML should bemathematically correct. All features implemented in
RETCML (the output listed in Section 3.2) have also
been tested. Note that the original RETC can crashwith ¯oating point errors when the estimation problem
is very ill-de®ned (e.g., when one tries to estimate theresidual saturation but only supplies the near-saturated
part of the retention curve). This problem still exists in
RETCML, because the optimization algorithm has notbeen changed. To report any other errors possibly stillin the code, please send an email to Jirka SÏimunek,
simunek@ussl.ars.usda.gov.

4.1. Programming details
RETCML was developed based on the DOS version of RETC. It has been compiled with Microsoft
Developer Studio, Fortran Power Station 4.0, on a
Windows NT 4.0 station. In the source code, themodi®cations made are in small letters, whereas
RETC is written in ALL CAPS. All modi®cations
are also preceded by comment lines to make thecode easier to understand. Some changes in theexisting RETC code were necessary for compilation.
These are: the standard input unit has been set to
5 instead of RETC's 0, an argument to dmax1 has
been made a double rather than a single precision
number (line 132 now), and the terminal (screen)
input unit for a character is  instead of 5 (line 60
now).
The original optimization algorithm was not changed. It was considered to introduce analytical objectivefunction gradient evaluation, but comparison with parallel MATLAB code revealed only very small improvements. The required programming eort to introduce
analytical gradients for all hydraulic parameter
schemes in RETC seemed unwarranted for the rather
small gain in accuracy.
Several new (with respect to RETC) functions/subroutines are included at the end of the RETCML
source code. These implement the incomplete gammafunction needed for computing the w 2 statistic, which
in turn is required for computing the probability ofmodel adequacy [see Eq. (6)]. The functions are called
gammq and gammln, the subroutines gser and gcf.
All these were taken from Press et al. (1992, Chap.
6.2). The function Q ÿ1() [see Eq. (8)] is not coded,
but rather contained in the form of a lookup table for
1±6 degrees-of-freedom at 90% con®dence (lines 35±37of the code).

K.J. Hollenbeck et al. / Computers & Geosciences 26 (2000) 319±327

References
Abramowitz, M., Stegun, I.A., 1965. Handbook ofMathematical Functions. Dover, New York (1046 pp).
Bard, Y., 1974. Nonlinear Parameter Estimation. Academic
Press (341 pp).
Beck, J.V., Arnold, K.J., 1977. Parameter Estimation in
Engineering and Science. Wiley, New York (501 pp).
Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology Paper 3. Colorado State University,
Fort Collins, Colorado (27 pp).
Burdine, N.T., 1953. Relative permeability calculations frompore-size distribution data, Petroleum Transactions.
American Institute of Mineral Engineering 198, 71±77.
Hollenbeck, K.J., Jensen, K.H., 1998. Maximum-likelihood
estimation of unsaturated hydraulic parameters. Journal of
Hydrology 210, 192±205.
Kool, J.B., Parker, J.C., van Genuchten, M.Th, 1985.
ONESTEP: a nonlinear parameter estimation program forevaluating soil hydraulic properties from one-step out¯owexperiments.
Bulletin
85-3,
Virginia
Agricultural
Experimental Station, Blacksburg, VA (43 pp).
Kool, J.B., Parker, J.C., van Genuchten, M.Th, 1987.
Parameter estimation for unsaturated ¯ow and transportmodels Ð a review. Journal of Hydrology 91, 255±293.
Leij, F.J., Bradford, S.A., 1994. 3DADE: a computer program for evaluating three-dimensional equilibrium solutetransport in porous media. Research Report 134, US
Salinity Laboratory, US Department of Agriculture, ARS,
Riverside, CA (81 pp).
Mualem, Y., 1976. A new model for predicting the hydraulicconductivity of unsaturated porous media. Water
Resources Research 12 (3), 513±522.
Press, W.S., Teukolsky, S.A., Vetterling, W.T., Flannery,
B.P., 1992. Numerical Recipes in C, 2nd ed. Cambridge
University Press, Cambridge (994 pp).
Schweppe, F.C., 1973. Uncertain Dynamic Systems. PrenticeHall, Englewood Clis, NJ (563 pp).
SÏimunek, J., SÏejna, M., van Genuchten, M.Th, 1996. The
HYDRUS-2D software package for simulating water ¯ow
and solute transport in two-dimensional variably saturated
media, version 1.0, IGWMC-TPS-56. International

327

Ground Water Modeling Center, Colorado School ofMines, Golden, CO (167 pp).
SÏimunek, J., SÏejna, M., van Genuchten, M.Th, 1998. The
HYDRUS-1D software package for simulating water ¯ow
and solute transport in one-dimensional variably saturated
media, version 2.0, IGWMC-TPS-70. International
Ground Water Modeling Center, Colorado School ofMines, Golden, CO (178 pp).
Toorman, A.F., Wierenga, P.J., Hills, R.G., 1992. Parameter
estimation of hydraulic properties from one-step out¯owdata. Water Resources Research 28 (11), 3021±3028.
Toride, N., Leij, F.J., van Genuchten, M.Th, 1995. TheCXTFIT code for estimating transport parameters fromlaboratory or ®eld tracer experiments, version 2.0.
Research Report 137, US Salinity Laboratory, US
Department of Agriculture, ARS, Riverside, CA (121 pp).
van Dam, J.C., Stricker, J.N.M., Droogers, P., 1994. Inverse
method to determine soil hydraulic functions from multistep out¯ow experiment. Soil Science Society of America
Journal 58, 647±652.
van Genuchten, M.Th, 1978. Calculating the unsaturated
hydraulic conductivity with a new closed-form analyticalmodel. Technical Report 78-WR-08. Water Resources
Program, Department of Civil Engineering, Princeton
University, Princeton, NJ (63 pp).
van Genuchten, M.Th, 1980a. A closed form equation forpredicting the hydraulic conductivity of unsaturated soils.
Soil Science Society of America Journal 44, 892±898.
van Genuchten, M.Th, 1980b. Determining transport parameters from solute displacement experiments. Research
Report 118, US Salinity Laboratory, US Department ofAgriculture, ARS, Riverside, CA (37 pp).
van Genuchten, M.Th, 1981. Non-equilibrium transport parameters from miscible displacement experiments. Research
Report 119, US Salinity Laboratory, US Department ofAgriculture, ARS, Riverside, CA (88 pp).
van Genuchten, M.Th, Leij, F.J., Yates, S.R., 1991. The
RETC code for quantifying hydraulic functions of unsaturated soils. Technical Report IAG-DW 12933934, US
Salinity Laboratory, US Department of Agriculture,
Agricultural Research Service, Riverside, CA (83 pp).