Soil & Tillage Research 47 (1998) 27±36

Parameter estimation of unsaturated soil hydraulicproperties from transient ¯ow processes1
JirÏÂõ SÏimuÊneka,*, Martinus Th. van Genuchtena, Molly M. Gribbb, Jan W. HopmanscbaU.S. Salinity Laboratory, USDA, ARS, 450 W. Big Springs Road, Riverside, CA 92507, USA
Department of Civil and Environmental Engineering, University of South Carolina, Columbia, SC 29208, USAc
Hydrology Program, Department LAWR, University of California, Davis, CA 95616, USA

Abstract
Three ®eld methods recently proposed for estimating the soil hydraulic properties by numerical inversion of the Richards'
equation are discussed. The ®rst method involves the use of tension disc permeameter data, while the second method uses datacollected with modi®ed cone penetrometer. The third method involves the use of a multiple-step ®eld extraction device.
Experimental data for each of the above three methods were analyzed by using the HYDRUS-2D code coupled with theLevenberg±Marquardt parameter estimation algorithm. Advantages and disadvantages of the three methods are discussed.
# 1998 Elsevier Science B.V. All rights reserved.
Keywords: Parameter estimation; Soil hydraulic properties; Tension disc permeameter; Cone penetrometer

1. Introduction
Numerical models are increasingly used to study
variably-saturated ¯ow and chemical transport processes between the soil surface and the groundwatertable. Effective use of these models depends to a large
extent on the accuracy of available soil water retention
and unsaturated conductivity data. Accurate measurement of the soil hydraulic properties is laborious
because of the nonlinear nature of these properties,
and the extreme spatial heterogeneity of the subsurface environment. While many laboratory and ®eld
methods exist to determine the soil water retention and*Corresponding author. Tel.: (909) 369 4865; fax: (909) 342
9964; e-mail: JSIMUNEK@USSC.ARS.USPA.GOV
1
Submitted to Soil Technology on September 11, 1997.
0167-1987/98/$19.00 # 1998 Elsevier Science B.V. All rights reserved.
PII S0167-1987(98)00069-5

unsaturated hydraulic conductivity curves (Klute andDirksen, 1986; Green et al., 1986), most methods
remain relatively expensive and too cumbersome forapplications to larger areas of land. Hence, relatively
cost-effective methods are needed for rapid determination of the hydraulic properties in the ®eld.
In this paper, we discuss three ®eld methods
recently proposed for estimating the soil hydraulicproperties by numerical inversion of the Richards'
equation. One method currently being developed
involves the use of tension disc permeameter data(SÏimuÊnek and van Genuchten, 1996, 1997), while asecond method uses data collected with a modi®ed
cone penetrometer (Gribb, 1996; Gribb et al., 1998).
The third method involves the use of a ®eld multiple
extraction device (Inoue et al., 1998). Experimental
data for each of the above three methods were ana-

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±36

28lyzed by coupling the HYDRUS-2D code (SÏimuÊnek
et al., 1996) with the Levenberg±Marquardt(Marquardt, 1963) parameter estimation algorithm.

m(1ÿ1/n)[-] are empirical parameters. The hydraulic characteristics de®ned by Eqs. (2) and (3) contain 5
unknown parameters: r , s, , n, and Ks.

2. Theory

2.2. Inverse solution

Inverse methods are typically based upon the minimization of a suitable objective function which
expresses the discrepancy between the observed
values and the predicted system response. Soil hydraulic properties are assumed to be described by ananalytical model with unknown parameter values,
e.g. by Eqs. (2) and (3) given below. The system
response is represented by a numerical solution ofthe ¯ow equation augmented with the parameterized
hydraulic functions and suitable initial and boundary
conditions. Initial estimates of the optimized systemparameters are then iteratively improved during theminimization process until a desired degree of precision is obtained.

The objective function  to be minimized during
the parameter estimation process may be de®ned as

2.1. Governing flow equation
The governing ¯ow equation for radially symmetricisothermal Darcian ¯ow in a variably-saturated isotropic rigid porous medium is given by the followingmodi®ed form of the Richards' equation:@ 1 @
@h
@
@h
@K

rK

K

(1)
@t r @r@r@z@z@z
where  is the volumetric water content [L3Lÿ3], h is
the pressure head [L], K is the hydraulic conductivity[LTÿ1], r is a radial coordinate [L], z is a verticalcoordinate [L] positive upward, and t is time [T].
The unsaturated soil hydraulic properties in thispaper are assumed to be described by the following
expressions (van Genuchten, 1980):
 h ÿ r1s ÿ r
1  j hjn m(2)

1=m m 2
K   Ks S0:5
e 1 ÿ 1 ÿ Se  (3)Se h where Se is the effective water content [-], Ks is thesaturated hydraulic conductivity [LTÿ1], r and sdenote residual and saturated water contents[Lÿ1], n[-], and[L3Lÿ3], respectively, and b; q; p 

mq
X

vj

nqj
X

wi;j qj x; ti  ÿ qj x; ti ; b2

j1
i1
mp
npj
X
X

vj



j1



nb
X
j1

i1

 i;j pj i  ÿ pj i ; b2w

^vj bj ÿ bj 2(4)

where the ®rst term on the right-hand side represents
the deviations between the measured and calculated
space-time variables (e.g. observed pressure heads orwater contents at different locations and/or time, or theactual or cumulative in®ltration or extraction rateversus time). In this term, mq is the number of different
sets of measurements, and nqj is the number of measurements in a particular measurement set. qj x; ti 
represents speci®c measurements at time ti for the jthmeasurement set at location x(r,z), qj(x,ti,b) are thecorresponding model predictions for the vector of
optimized parameters b (e.g. r, s, , n, and Ks),
and vj and wi,j are weights associated with a particularmeasurement set or point, respectively. The second
term of Eq. (4) represents differences between independently measured and predicted soil hydraulic properties (e.g. retention, (h), hydraulic conductivity, K()or K(h), and/or diffusivity, D() or D(h), data), while i;j havethe terms mp, npj, pj (i), pj(i,b), vj and w
similar meanings as for the ®rst term but now for
the soil hydraulic properties. The last term of Eq. (4)represents a penalty function for deviations between
prior knowledge of the soil hydraulic parameters, bj ,and their ®nal estimates, bj, with nb being the numberof parameters with prior knowledge and ^vj representing pre-assigned weights. Estimates which make useof prior information (such as those used in the third
term of Eq. (4)) are known as Bayesian estimates. We
note that the covariance (weighting) matrices which
provide information about the measurement accuracy,

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±36

as well as any possible correlation between measurement errors and/or parameters, are assumed to bediagonal in this study.
Minimization of the objective function  is accomplished by using the Levenberg±Marquardt nonlinear
minimization method (Marquardt, 1963), which combines the Newton and steepest descend methods, andgenerates con®dence intervals for the optimized parameters.
3. Tension disc infiltrometer
Tension in®ltration data have, thus, far been usedprimarily for evaluating saturated and unsaturated
hydraulic conductivities, and for quantifying theeffects of macropores and preferential ¯ow paths onin®ltration (Perroux and White, 1988; Ankeny et al.,
1991; among many others). A relatively standard way
for estimating unsaturated hydraulic conductivities
from tension in®ltrometer data has been to invoke
Wooding's, 1968 analytical solution. This approachrequires steady-state in®ltration rates for two different
supply pressure heads, and assumes applicability ofGardner's, 1958 exponential function for K(h). SÏimuÊnek and van Genuchten, 1996 suggested the combined
use of transient in®ltration data obtained during asingle tension in®ltration experiment, and tensiometeror TDR data measured in the soil below the disc toestimate the unknown soil hydraulic parameters in
Eqs. (2) and (3) via parameter estimation. SÏimuÊnek
and van Genuchten, 1997 later revised this method byusing multiple tension in®ltration experiments in combination with knowledge of the initial and ®nal watercontents. This modi®cation avoided the cumbersome
use of tensiometers and TDRs. Evaluation of thenumerical stability and parameter uniqueness using
numerically generated data with superimposed stochastic and deterministic errors showed, that a combination of the multiple cumulative tension in®ltration
data, a measured ®nal water content, and an initialcondition expressed in terms of the water content,
provided the most promising parameter estimation
approach for practical applications (SÏimuÊnek andvan Genuchten, 1997).
The numerical inversion method was applied totension disc in®ltrometer data obtained in two ®eld
studies. In the ®rst study (SÏimuÊnek et al., 1998a), two29

tension disc experiments were carried out on a ®nesandy loam in Riverside, California. The tension dischad a diameter of 20 cm, and the experiment was
executed with three consecutive tensions of 20, 10,
and 3 cm. Each tension was maintained for about 1 h.
Pressure transducers were used to record transientin®ltration rates in a setup similar to that describedby Ankeny et al., 1988. The soil hydraulic properties
were also independently measured with standard
laboratory methods on ®ve soil samples. Fig. 1 showsthe experimental and ®tted cumulative in®ltration
curves versus time for the two runs, including deviations between the best-®t and measured values. The
soil hydraulic parameters determined by parameter
estimation and by Wooding's analysis are given in
Table 1.
Excellent agreement between the measured and®tted cumulative ®eld in®ltration curves was obtained
when the ®ve major soil hydraulic parameters in van
Genuchten's model were optimized (Fig. 1). Hydraulic conductivities (both saturated and unsaturated) as
determined by the numerical inversion compared wellwith those obtained using Wooding's analytical solution (Table 1). Agreement between retention curves
obtained by numerical inversion of the ®eld in®ltration
experiment and by direct laboratory measurement,
however, was relatively poor. The entire laboratory
retention curve shifted by about 0.10 units towards
higher water contents, while also the parameter n had
different values ( -values were much closer for bothanalyses). The in®ltration experiment could not be
closely reproduced when the mean laboratory soil
hydraulic parameters were used and validity of theclosed-form van Genuchten±Mualem model was
assumed (SÏimuÊnek et al., 1998a).
The second experiment was used to estimate the soil
hydraulic characteristics of a two-layered soil system
involving a crusted soil in a Sahel region (SÏimuÊnek etal., 1998b). Here, we will report only results for the
sandy subsoil obtained with a tension disc diameter of
25 cm and with supply tensions of 11.5, 9, 6, 3, 1, and
0.1 cm. Fig. 2 shows measured and optimized cumulative in®ltration curves and their differences. The
small breaks in the cumulative in®ltration curve(Fig. 2(a)) were caused by brief removal of the in®ltrometer from the soil surface to resupply it with waterand to adjust the tension for a new time interval. Veryclose agreement between the measured and optimized

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±3630Fig. 1. Measured and fitted cumulative infiltration curves (a) and their differences (b) for two tension disc infiltrometer experiments carriedout in Riverside, California.

Table 1Comparison of results obtained with Wooding's analysis and using parameter estimation for two multiple tension infiltration experiments
Parameter

Run 1

Run 2

Wooding's analysis(cmÿ1)nrs
Ks (cm dÿ1)
K(ÿ6.5 cm) (cm dÿ1)
K(-15 cm) (cm dÿ1)
a

32.0
6.81

14.0
9.65

Parameter
estimation
0.0412
3.61
0.0
0.274
9.24 (ÿ34%) a
8.80 (ÿ8.8%)
4.95 (ÿ27%)

Wooding's analysis

20.4
3.74

12.3
7.26

Parameter
estimation
0.0406
2.76
0.0
0.238
6.57 (ÿ47%)
5.34 (ÿ26%)
2.49 (ÿ33%)

Mean lab.
curve

0.0401
1.55
0.107
0.347
272

The number in parentheses shows by how many percent the hydraulic conductivity obtained by parameter estimation differs fromcorresponding values obtained using Wooding's analysis.

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±3631Fig. 3. Unsaturated hydraulic conductivities calculated using
Wooding's analytical solution for particular pressure heads, andthe complete function obtained with numerical inversion.

Fig. 2. Measured and optimized cumulative infiltration curves (a)and their differences (b) for a tension disc infiltrometer experimentcarried out in the Sahel region.

cumulative in®ltration curves was obtained; the largest deviations were generally less than 60 ml, which
was only about 0.5% of the total in®ltration volume.
Fig. 3 shows a comparison of parameter estimation
results against results obtained with Wooding's analysis. Both methods give almost identical unsaturated
hydraulic conductivities for pressure heads in theinterval between ÿ2 and ÿ10.25 cm. However, thehydraulic conductivity at the highest pressure headinterval was overestimated by a factor of two using
Wooding's analysis. In Fig. 4, we further compare thenumerical inversion results with hydraulic properties
estimated from available soil textural information
using a neural-network-based pedotransfer function
approach (Schaap et al., 1998). Relatively good agree-

Fig. 4. Soil hydraulic functions obtained with numerical inversion
and using neural network predictions.

32

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±36

ment between the inverse and neural network predictions was obtained.
4. Cone penetrometer
While tension in®ltrometer experiments provide
relatively quick estimates of the hydraulic properties,
they can be used only at the soil surface. By comparison, a new cone penetrometer method currently underdevelopment (Gribb, 1996; Gribb et al., 1998; KodesÏovaÂ et al., 1998) can be used at depth. Cone penetrometers were originally used to obtain soil strengthcharacteristics by measuring the tip resistance andsleeve friction during penetration at a constant rate.
To obtain the hydraulic properties, a modi®ed cone
penetrometer, instrumented with a porous ®lter close
to the penetrometer tip and two tensiometer rings 5
and 9 cm above the ®lter, is used (Fig. 5). The device
is pushed into a soil to the desired depth, and aconstant head is applied to the 5 cm ®lter. The volume
of water imbibed into the soil is monitored, as are
tensiometer ring readings registering the advancementof the wetting front for a short period of time (300±500 s).
Gribb, 1996 gave a detailed numerical analysis of
this experiment, including a study of the identi®abilityof the soil hydraulic parameters. She showed that theinverse solution was least sensitive to n and s, andFig. 5. Schematic of the modified cone penetrometer (Gribb,
1996).

most sensitive to Ks and . The method was recently
used to estimate the hydraulic parameters of a sandy
soil in a laboratory aquifer system measuring
553 m (Gribb et al., 1998; KodesÏovaÂ et al., 1998).
Fig. 6 shows observed ¯ow data, as well as resultsof the numerical inversion. Excellent agreement
between measured and optimized values was obtained
for the inverse solution with four optimized parameters. Fig. 7 presents the retention curves obtainedwith selected laboratory methods and the parameter
estimation technique. Optimized curve is close to thewetting curves determined in the laboratory. The
overprediction of s for the inverse solution could
be caused by anisotropy of the soil. Gribb et al.,1998 showed that the estimated saturated hydraulicconductivities were similar to those obtained with

Fig. 6. Comparison of observed and optimized cumulative infiltration curves and tensiometer readings for the modified cone penetrometer test
( , n, s, and Ks are optimized, r is fixed).

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±3633Fig. 7. Comparison of soil-water characteristic curves obtained bydirect methods and that obtained by inversion of cone penetrometerdata (CR ± capillary rise test, CA ± computer automated wetting
curve test, PP ± pressure plate test, wc ± wetting curve, dc ± drying
curve).

other test methods, such as Guelph permeameter, slug
tests, and laboratory constant head tests. Issues related
to soil disturbance due to cone penetration are underinvestigation (Gribb et al., 1998).
5. Multiple step extraction
While the tension in®ltrometer and cone penetrometer methods provide information about wetting
branches of the soil hydraulic properties, a multiple
step extraction device can be used to obtain draining
branches. The device consists of a ceramic soil solution sampler (Fig. 8) which is inserted into an initially
wet soil pro®le and subjected to a series of vacuum
extraction pressures. The cumulative amount of soil
solution extracted during an experiment, as well aspressure heads at various locations near the extraction
device, are monitored during the experiment andsubsequently used in an objective function for the
nonlinear minimization problem. Inoue et al., 1998
®rst evaluated the feasibility of the vacuum extraction
technique using numerically generated data, and concluded that the method is well suited for loamytextured soils, but not necessarily for sandy soils.
They tested the method in the laboratory and in the®eld. Here we brie¯y discuss one ®eld application.
The experiment (Inoue et al., 1998) was carried out
on a Yolo silt loam with a clay content of about 22%.

Fig. 8. Schematic of the in situ multistep extraction experiment(Inoue et al., 1998).

The center of the ceramic ring (with a radius of 3 cm,
and a length of 3 cm) was located 10 cm below the soil
surface. The tensiometers were installed at the following positions: T1(r,z)(4, ÿ10 cm), T2(6, ÿ15),
T3(20, ÿ40), and T4(20, ÿ55). Readings of theforth tensiometer (T4) were used as the bottom boundary condition. The initially saturated soil was allowed
to drain for a period of 46.5 h, after which threevacuum extraction steps were applied: hexÿ195 cm hfor 46.5<t<71.2 h, hexÿ415 cm for 71.2<t<93.0 h,
and hexÿ685 cm for 93.0<t<120 h. Experimental
data, as well as the ®nal results of two numericalinversions are presented in Fig. 9. The saturated
hydraulic conductivity of the ceramic ring was optimized simultaneously with van Genuchten's hydraulicparameters. Only the cumulative extraction volume
versus time and transient tensiometer measurements
were included in the objective function for the ®rst
optimization. Three additionally measured (h) data
points were used for the second optimization. Agreement between measured and calculated values inFig. 9 was relatively good. Fig. 10 compares the
optimized soil hydraulic functions with independently
measured data obtained with the instantaneous pro®lemethod. While correspondence between the measured

34

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±36

Fig. 9. Comparison of measured and optimized cumulative extraction and pressure head values. The objective function included cumulative
extraction volumes, tensiometer measurements (a and b), and additionally 3 measured (h) data points (c and d).

and estimated hydraulic conductivities seems good forboth optimizations, only the second optimization gave
acceptable results for the retention curve (Inoue et al.,1998).
6. Conclusions
We brie¯y reviewed three ®eld methods which may
be used for estimating the soil hydraulic properties byparameter estimation. A common feature of all threemethods is the coupling of a numerical solution of the
Richards equation with the Levenberg±Marquardt
parameter estimation algorithm. Of the three methods,
the tension disc in®ltrometer and cone penetrometer
experiments provide information of the wetting

branches of the soil hydraulic properties, while themultiple extraction device leads to draining branches.
Use of the tension disc in®ltrometer and multiple
extraction device is limited to the soil surface or thenear surface. Although tested up to now only in thelaboratory, cone permeameter can potentially be used
to depths of 30 m or more in accordance with established cone penetration technology. All three methods,
except perhaps the tension in®ltrometer approach, arestill being perfected and require further research to test
their applicability in the ®eld.
Parameter estimation approaches as discussed here
provide unique tools for evaluating the design andperformance of these and related experimental
approaches, and for optimally analyzing data collectedwith the different ®eld instruments.

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±3635Fig. 10. Estimated soil water retention (a) and hydraulic conductivity (b) functions for the two optimizations of multi-step extraction datacompared with independently determined retention and unsaturated hydraulic conductivity data.

References
Ankeny, M.D., Kaspar, T.C., Horton, R., 1988. Design for an
automated tension infiltrometer. Soil Sci. Sco. Am. J. 52, 893±
896.
Ankeny, M.D., Ahmed, M., Kaspar, T.C., Horton, R., 1991. Simple
field method for determining unsaturated hydraulic conductivity. Soil Sci. Soc. Am. J. 55, 467±470.
Gardner, W.R., 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporationfrom a water table. Soil Sci. 85, 228±232.

Green, R.E., Ahuja, L.R., Chong. S.K., 1986. Hydraulic conductivity, diffusivity, and sorptivity of unsaturated soils: Field
methods. In: Klute, A. (Ed.), Methods of Soil Analysis, Part 1 ±Physical and Mineralogical Methods, 2nd ed. Agronomy
Monograph 9, pp. 771±798.
Gribb, M.M., 1996. Parameter estimation for determining hydraulicproperties of a fine sand from transient flow measurements.
Water Resour. Res. 32(7), 1965±1974.
Gribb, M.M., SÏimuÊnek, J., Leonard, M.F., 1998. Use of cone
penetrometer method to determine soil hydraulic properties, J.
Geotech. and Geoenviron. Eng., in press.

36

J. SÏimuÊnek et al. / Soil & Tillage Research 47 (1998) 27±36

Inoue, H., SÏimuÊnek, J., Hopmans, J.W., and Clausnitzer, V., 1998.
In-situ estimation of soil hydraulic functions using a multi-stepwater extraction technique, Water Resour. Res., submitted.
Klute, A., Dirksen, C., 1986. Hydraulic conductivity anddiffusivity: Laboratory methods. In: Klute, A. (Ed.), Methodsof Soil Analysis, Part 1 ± Physical and Mineralogical Methods,
2nd ed. Agronomy Monograph 9, pp. 687±734.
KodesÏovaÂ, R., Gribb, M.M., SÏimuÊnek, J., 1998. A new CPT method
for estimating soil hydraulic properties, Proc. of ISC 98,
Atlanta, GA, submitted.
Marquardt, D.W., 1963. An algorithm for least-squares estimationof nonlinear parameters. SIAM J. Appl. Math. 11, 431±441.
Perroux, K.M., White, I., 1988. Design for disc permeameters. Soil
Sci. Soc. Am. J. 52, 1205±1215.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 1998. Neuralnetwork analysis for hierarchical prediction of soil water
retention and saturated hydraulic conductivity. Soil Sci. Soc.
Am. J., in press.
SÏimuÊnek, J., SÏejna, M., van Genuchten, M.Th., 1996. The
HYDRUS-2D software package for simulating water flow
and solute transport in two-dimensional variably saturated
media. Version 1.0. IGWMC ± TPS ± 53, International Ground

Water Modeling Center, Colorado School of Mines, Golden,
Colorado.
SÏimuÊnek, J., van Genuchten, M.Th., 1996. Parameter estimation ofsoil hydraulic properties from the tension disc infiltrometer
experiment by numerical inversion. Water Resour. Res. 32(9),2683±2696.
SÏimuÊnek, J., van Genuchten, M.Th., 1997. Parameter estimation ofsoil hydraulic properties from multiple tension disc infiltrometer data. Soil Sci. 162(6), 383±398.
SÏimuÊnek, J., Wang, D., Shouse, P.J., van Genuchten, M.Th., 1998a.
Analysis of a field tension disc infiltrometer experiment byparameter estimation. Inter. Agrophysics, submitted.
SÏimuÊnek, J., Angulo-Jaramillo, R., Schaap, M., Vandervaere, J.-P.,
van Genuchten, M.Th., 1998b. Using an inverse method toestimate the hydraulic properties of crusted soils from tension
disc infiltrometer data. Geoderma, in press.
van Genuchten, M.Th., 1980. A closed-form equation forpredicting the hydraulic conductivity of unsaturated soils. Soil
Sci. Soc. Am. J. 44, 892±898.
Wooding, R.A., 1968. Steady infiltration from large shallow
circular pond. Water Resour. Res. 4, 1259±1273.