Published online July 18, 2005

Inverse Analysis of Upward Water Flow in a Groundwater Table Lysimeter
T. J. Kelleners,* R. W. O. Soppe, J. E. Ayars, J. Šimůnek, and T. H. Skaggs

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT

scale often leads to disappointing results. Soil disturbance during sampling and spatial variability in the soil
hydraulic properties play major roles in the failure toobtain suitable field-scale parameter values for the hydraulic functions (Kool et al., 1987).
Alternatively, soil hydraulic properties may be determined in situ. This avoids soil disturbance and provides
an estimate of the hydraulic properties that integrates
soil heterogeneity. Also, if the soil hydraulic functions
are ultimately to be used to analyze field-scale processes,
in situ estimation intuitively seems more appropriate
than laboratory analysis.
Parameter optimization has emerged as an important
technique for estimating soil hydraulic parameters. Bothlaboratory and field experiments can now be run andanalyzed under a large variety of flow conditions. In
the past, these experiments were restricted to certain
well-defined conditions that allow the calculation of thehydraulic parameters by solving analytical solutions of
the relevant flow equations. Today, the hydraulic parameters can be obtained from flow experiments bycombining a numerical flow model with a parameter
estimation code. Recent advances in the computational
power of computers have made this method even more
attractive.
Numerous studies have been published on the combined use of flow experiments and parameter optimization to estimate the soil hydraulic parameters. Reviewsof these studies can be found in Kool et al. (1987),
Hopmans and Šimůnek (1997), and Hopmans et al.
(2002). A recurring issue is the problem of parameter
uniqueness. If the objective function to be minimized(usually the sum of squared differences between measured and calculated flow variables) does not have aclear global minimum, the inversion will be nonunique.
The occurrence of local minima in the objective function
may also cause nonuniqueness. The occurrence of uniqueness problems depends on the soil type, the boundary
conditions, the type of data used in the objective function, and the parameter estimation algorithm.
Several investigators have tried to infer soil hydraulicproperties by applying the inverse modeling technique
to upward infiltration experiments. The experiments
provide information on the wetting branch of the soil
hydraulic properties without the need to account forthe macropore flow that might occur during downward(gravity) infiltration. So far, several laboratory experiments on small soil samples have been conducted using
either a flux condition (Hudson et al., 1996) or a pressure
condition (Šimůnek et al., 2001; Young et al., 2002) atthe bottom boundary.
In this work, the hydraulic properties of a silty clay

The accuracy of numerical water flow models for the vadose zone
depends on the estimation of the soil hydraulic properties. In thisstudy, the hydraulic parameters for a silty clay soil in a large lysimeter
were determined through inverse modeling of a fallow period withupward water flow from a shallow groundwater table. Parameter
uniqueness was studied by simulating a hypothetical soil with known
hydraulic properties under comparable conditions. Sensitivity analysis
showed that the pressure head h(z,t ), the volumetric water content
␪(z,t ), and the cumulative bottom flux Q(t ) were least sensitive tothe residual volumetric water content ␪r and the pore-connectivity
parameter ␭ in the van Genuchten–Mualem (VGM) model. Parameter response surfaces showed that least squares fitting with ␪(z,t ) data
is more likely to result in a unique hydraulic parameter set thanleast squares fitting with h(z,t ) or Q(t ) data. With only ␪(z,t ) in
the objective function, the least squares minimization algorithm wascapable of finding the correct soil hydraulic parameters, provided that␪r and ␭ were fixed and that multiple initial parameter estimates wereused. The protocol that was developed for the hypothetical soil was
subsequently applied to the actual groundwater table lysimeter. The
soil hydraulic parameters for the lysimeter for two (x,y ) locations
were determined using ␪(z,t ) data as measured by capacitance sensors.
The variability in the optimized inverse of the air-entry value ␣ andthe saturated hydraulic conductivity Ks in the VGM model was relatively high because of the high parameter correlation between theseparameters. The optimized soil hydraulic properties can be used to
study capillary rise from the groundwater table.

Environmental impact and irrigation water management studies often rely on numerical models to
predict water flow and solute transport in the vadosezone. The predictions of these models are sensitive tothe input of the soil hydraulic properties. The soil hydraulic properties are usually described using empiricalor semitheoretical water retention and hydraulic conductivity functions (e.g., Gardner, 1958; Brooks andCorey, 1964; Mualem, 1976; van Genuchten, 1980; Kosugi, 1999).
Traditionally, hydraulic properties have been measured in the laboratory using soil samples taken fromthe field. Applying the soil hydraulic properties derived
from these small-scale experiments to the larger field
T.J. Kelleners, J. Šimůnek, and T.H. Skaggs, USDA-ARS, George
E. Brown, Jr. Salinity Lab., 450 W. Big Springs Road, Riverside, CA
92507; T.J. Kelleners, R.W.O. Soppe, and J.E. Ayars, USDA-ARS,
Water Management Research Lab., 9611 S. Riverbend Ave., Parlier,CA 93648. T.J. Kelleners, presently at Dep. of Plants, Soils, andBiometeorology, Utah State Univ., Logan, UT 84322; R.W.O. Soppe,
presently at Alterra-ILRI, P.O. Box 47, 6700 AA Wageningen, The
Netherlands; J. Šimůnek, presently at Dep. of Environmental Sciences,
Univ. of California, Riverside, CA 92521. Received 9 Aug. 2004.
*Corresponding author (tkelleners@cc.usu.edu).
Published in Vadose Zone Journal 4:558–572 (2005).
Original Research Paper
doi:10.2136/vzj2004.0118© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: CV, coefficient of variability; EC, electrical conductivity; LM, Levenberg–Marquardt; NRMSE, normalized root mean squareerror; VGM, van Genuchten–Mualem.

558

559

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

soil in a weighing lysimeter are determined by inverse
modeling of a summer fallow period where capillary risefrom the groundwater table replenishes the depleted rootzone. The lysimeter is 4 m long, 2 m wide, and 3 m
deep, and as such, constitutes an intermediate scale between the sample scale and the field scale. The Hydrus1D model (Šimůnek et al., 1998a), which includes aparameter estimation routine, is used to solve the flow
problem. The objectives of the study were (i) to investigate the issue of parameter uniqueness for the upward
infiltration problem by testing the inverse procedure ona hypothetical soil with known soil hydraulic propertiesand (ii) to determine the soil hydraulic properties ofthe silty clay soil in the lysimeter.

THEORY
Governing Flow Equation
The water flow calculations are conducted with theHydrus-1D model, which numerically solves the Richards equation:

冢

⳵h
⳵
⳵␪
K
⫽
⫹K
⳵t
⳵z
⳵z

冣

[1]

where ␪ is the volumetric water content (L3 L⫺3), t is
the time (T), z is the vertical coordinate (L) (positive
upward), h is the pressure head (L), and K is the hydraulic conductivity (L T⫺1).
Initial and boundary conditions need to be specified
to solve Eq. [1]. The initial condition is specified interms of pressure head or water content:h(z,t) ⫽ hi(z)

t⫽0

[2]

␪(z,t) ⫽ ␪i(z)

t⫽0

[3]

The soil surface boundary is described by an atmospheric boundary condition that switches between a prescribed flux condition and a prescribed head condition,
depending on the prevailing transient pressure head at
the soil surface. Whenever the surface pressure head isin the range hA ⬍ h(0,t) ⬍ hS, the atmospheric boundary
is a prescribed flux condition:兩⫺K ⳵h⳵z ⫺ K 兩 ⫽ R

z⫽0

[4]

where R is the potential rate of infiltration or evaporation (L T⫺1) under the prevailing atmospheric conditions, and hA (L) and hS (L) are minimum and maximum
allowed values of the surface pressure head, respectively. If the surface pressure reaches hA or hS, the surface boundary switches to a prescribed pressure head
condition: h(0,t) ⫽ hA or h(0,t) ⫽ hS. In this study, we
use hA ⫽ ⫺105 cm (hypothetical soil), hA ⫽ ⫺106 cm(lysimeter soil), and hS ⫽ 0 cm.
A pressure head boundary condition is specified at
the bottom of the lysimeter:h(z,t) ⫽ h0(t)z ⫽ ⫺L[5]

where h0 is the prescribed value of the pressure head(L) and L is the depth (L) of the lysimeter.

Soil Hydraulic Properties
The soil hydraulic properties are described with the
VGM model (van Genuchten, 1980; Mualem, 1976):1h⬍0␪ ⫺ ␪r⫽ (1 ⫹ |␣h|n )m[6]
Se ⫽␪s ⫺ ␪ r1hⱖ0冦K(Se) ⫽ KsSe␭ 共1 ⫺ 关1 ⫺ S e1/m兴 兲

m 2

[7]

m ⫽ 1 ⫺ 1/nn⬎1[8]
where Se is the effective saturation, ␪r is the residual
volumetric water content, ␪s is the saturated volumetricwater content, ␣ is the inverse of the air-entry value(L⫺1), n is a pore-size distribution index, Ks is the saturated hydraulic conductivity (L T⫺1), and ␭ is a poreconnectivity parameter.

Inverse Procedure
The objective function ⌽ that is minimized during
the parameter estimation process is (e.g., Šimůnek etal., 1998a):mq

nqk oqjk

k⫽1

j⫽1 i⫽1⌽(b;q) ⫽ 兺 wk 兺 兺 [q*k (zj,ti) ⫺ qk (zj,ti,b)]2[9]

where b is the vector of parameters to be optimized(e.g., b ⫽ [␪r, ␪s, ␣, n, Ks, ␭]), q*k is the measured value
for the kth measurement set at depth zj and at time ti,
qk is the corresponding predicted value for parameter
vector b, mq is the total number of measurement sets,
nqk is the number of depths z for the kth measurement
set, and oqjk is the number of observation times t for
depth z and measurement set k.
The different measurement sets are weighted using
weighting coefficients wk:
1[10]
wk ⫽ nqk
␴2k 兺 oqjk
j⫽1

where ␴2k is the variance of the data in the kth measurement set. The use of wk minimizes unwanted differencesin weighting among data types caused by differences inmagnitude and the number of data points.
Minimization of the objective function ⌽ is accomplished by using the Levenberg–Marquardt (LM) nonlinear minimization method (Marquardt, 1963; Šimůnek
et al., 1998a). The LM method uses the steepest descentmethod when the objective function is far from its minimum and switches to the inverse Hessian method as theminimum is approached (Šimůnek and Hopmans, 2002).
The LM method is a local gradient-type search algorithm, as opposed to global search algorithms that search
the entire parameter space. Local search algorithms are
generally sensitive to the initial parameter estimates.
Consequently, different initial estimates need to be examined.

HYPOTHETICAL SOIL
Model Setup
Inverse estimation procedures for the lysimeter were
investigated by considering a hypothetical soil with knownReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

560

VADOSE ZONE J., VOL. 4, AUGUST 2005

soil hydraulic properties. Upward flow from a shallow
groundwater table was simulated for conditions closely
resembling the summer fallow period in the lysimeter. A
homogeneous silty clay soil was assumed with a constantgroundwater table at 1 m below the soil surface. The
vertical profile was divided into 200 elements, each witha thickness of 0.5 cm. The parameters describing thehydraulic properties of the soil were obtained from the
Rosetta pedotransfer program of Schaap et al. (2001).
Soil texture data (2.6% sand, 43.9% silt, and 53.5%
clay) and dry bulk density data (1.38 g cm⫺3) were used
as input for Rosetta. These soil data were obtained
during a field experiment at the location where thelysimeter soil was taken (see Kelleners et al., 2004). The
data were averages for the B1 horizon of this soil (20–75
cm below soil surface). The Rosetta-predicted soil hydraulic parameter values were: ␪r ⫽ 0.101, ␪s ⫽ 0.492,
␣ ⫽ 0.015 cm⫺1, n ⫽ 1.321, Ks ⫽ 3.47 cm d⫺1, and ␭ ⫽⫺1.055. Note that the physical meaning of ␭ as a tortuosity factor is lost by allowing ␭ ⬍ 0. With ␭ ⫽ ⫺1.055 theVGM model should be interpreted as a semitheoreticalfitting equation (Schaap and Leij, 2000).
A 100-d fallow period was simulated with an initiallydry soil profile. The initial pressure head values wereassumed to increase linearly over five depth intervals
with hi(z ⫽ 0 cm) ⫽ hA ⫽ ⫺100 000 cm, hi(⫺10) ⫽⫺30 000 cm, hi(⫺50) ⫽ ⫺15 000 cm, hi(⫺70) ⫽ ⫺1000
cm, hi(⫺90) ⫽ ⫺30 cm, and hi(⫺100) ⫽ 0 cm. This initialpressure head distribution was intended to represent asoil profile that is depleted by root water uptake from a
preceding crop. The top 10 cm of the soil was consideredfurther depleted by continued evaporation from the soil
surface after irrigation was halted toward the end of thegrowing season. An atmospheric boundary condition
was specified at the top boundary with a constant potential evaporation rate of 0.5 cm d⫺1 and zero rainfall andirrigation. At the bottom boundary a groundwater tablecondition was specified (h0 ⫽ 0 cm).
The bottom boundary during the inverse analysis was
specified as a “time-variable boundary condition” inHydrus, albeit with constant zero pressure head. This
prevented negative pressure heads at the bottom boundary when the LM algorithm tested ␪s values that werehigher than the true ␪s value of 0.492. This problem
only occurred because ␪i was used to specify the initialcondition during the inverse analysis (with ␪i, each Hydrus run starts with the conversion of ␪i into hi using
the current estimate of the soil hydraulic parameters).
The problem is a consequence of the incompatibility
between ␪i(z) and h0(t) at t ⫽ 0 and z ⫽ ⫺L. Specifying
the bottom boundary as a time-variable boundary condition may still result in some artificial redistribution near
the bottom during the first few time steps, but the effectof the redistribution on the calculated water flow is negligible.

Data Generation
Simulated pressure heads and volumetric water contents at the end of each day at depths of 5, 15, 25, 35,
45, 55, 65, 75, 85, and 95 cm below the soil surface were

Fig. 1. Simulated volumetric water content for the 100-d fallow periodwith hypothetical soil hydraulic properties. The numbers in the
figure refer to the depth below the soil surface.

recorded for use in the inverse analysis. This resulted
in 1000 data points for ␪(z,t) and h(z,t) each. As an
example, the resulting volumetric water contents as afunction of time and depth are shown in Fig. 1. Notethat during the 100-d calculation period the upward
moving wetting front passes all measurement depths, thus
providing a maximum amount of information for the inverse analysis. In addition, the cumulative upward flux
through the bottom of the soil profile (22.1 cm) was subdivided into 1-d intervals, resulting in 100 Q(t) data points.
It should be pointed out that the setup for the hypothetical soil is relatively favorable for the inverse process
as compared with an actual lysimeter soil. First, we
assumed a homogeneous soil profile between 0 and 100
cm depth. In real-life situations the soil profile will beheterogeneous and may even encompass completely different soil layers. Second, we did not consider type I
model errors (e.g., Neuman, 2003). Flow in a real soil
will not adhere strictly to the Richards equation. For
example, preferential flow and vapor transport willcause deviations between field reality and the model.
Also, the hydraulic properties of a real soil might notalways fit the VGM model. Third, we did not consider
measurement errors during most of the inverse calculations. In practice, instrument errors and observation errors
will result in a certain degree of scatter in the data, aswill be shown during the stability analysis. Fourth, thelarge range of pressure head values in the initial soil
profile is beneficial for the parameter identifiability during the inverse process. In the field, the pressure headgradients may be smaller. Therefore, the hypotheticalsoil constitutes a “best case” scenario with respect toinverse analysis.

Sensitivity Analysis
A sensitivity analysis was performed to study the relationship between the measurement data and the model
parameters for the upward flow problem. The higher thesensitivity of a model parameter to the data, the higher
the chance that the parameter is identifiable during theinverse process. Sensitivity coefficients for the hypothetical soil were calculated from (Šimůnek et al., 1998b):

561

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 2. Sensitivity of the pressure head in the hypothetical soil to a
1% change in the soil hydraulic parameters in the van Genuchten–
Mualem model at depths of 5, 15, and 35 cm below the soil surface.

兩skl ⫽ 0.01bl

兩 兩

⳵qk
q (b ⫹ ⌬bel) ⫺ qk (b)⫽ 0.01bl k
⳵bl
1.01bl ⫺ bl

兩 [11]⫽ 冨qk (b ⫹ ⌬bel) ⫺ qk (b)冨where skl is the change in variable qk corresponding to
a 1% change in parameter bl, el is the lth unit vector,and ⌬b ⫽ 0.01 b. The 0.01bl term in Eq. [11] is included
to allow comparison between different parameters, independent of their invoked unit or absolute value.
Sensitivities skl of the pressure head h(z,t), the volumetric water content ␪(z,t), and the cumulative bottom
flux Q(t) were calculated for all six soil hydraulic parameters. All three measurement sets showed a decrease insensitivity according to n ⬎ ␣ ⬎ ␪s ⬎ Ks ⬎ ␭ ⬎ ␪r, withn being the most influential parameter and ␪r being theleast influential parameter. As an example, Fig. 2 shows
the sensitivity of the pressure head at the 5-, 15-, and35-cm depths to the six hydraulic parameters. Maximum
sensitivity is observed at times when the upward moving
wetting front passes a certain depth and large changes
in the pressure head occur. Sensitivities are highest forthe shallowest depth; this is due primarily to the nonlinear nature of the retention curve (small changes in thewater content of relatively dry soil result in large changes
in the pressure head). By comparison, the differences
are less pronounced for the sensitivity of the water contents at different depths (data not shown). Figure 2 also
shows that the duration of elevated sensitivity levels islongest for the shallowest depth. This is attributed to themore diffuse moisture front at this depth, a consequenceof the increased distance from the groundwater table.
The sensitivity of the cumulative bottom flux to thehydraulic parameters is shown in Fig. 3. Clearly, thesensitivity increases as the simulation progresses in time.
However, most of the increase in sensitivity occurs during the first 10 to 30 d. Figures 2 and 3 show that n
is the most influential parameter while ␣, ␪s, and Ks
constitute a middle group. The flow variables are least
influenced by ␭ and ␪r, which describe the dry part of
the hydraulic conductivity curve and the dry part of thewater retention curve, respectively. However, interpretation of the sensitivities for the different parameters

Fig. 3. Sensitivity of the cumulative upward bottom flux in the hypothetical soil to a 1% change in the soil hydraulic parameters in thevan Genuchten–Mualem model.

is not straightforward. The sensitivity value should be
evaluated against the nature of the parameter. For example, a 1% change in the value of Ks, which can change
several orders of magnitude, is less significant than a
1% change in the value of n, which changes only between
1.0 and 3.0 for most natural soils (Šimůnek and van
Genuchten, 1996). Keeping the above in mind, the sensitivities seem to indicate that the upward flow problem
is more suitable to identify n, ␣, ␪s, and Ks, and lesssuitable to identify ␭ and ␪r. A flow problem involving
a drying soil will probably be more suitable for theidentification of ␭ and ␪r.

Response Surfaces
The uniqueness of the inverse problem was investigated by calculating two-dimensional response surfacesof the objective function ⌽ as a function of pairs of soil
hydraulic parameters (e.g., Toorman et al., 1992). Each
response surface was created by varying the two selectedparameters around their true value using 50 discrete
points, while keeping the other parameters constant.
This resulted in 2500 simulated ⌽ values for each response surface. The six soil hydraulic parameters could
be paired in 15 different ways, yielding 15 response
surfaces for each case considered. Objective functions
were calculated for ␪(z,t), h(z,t), Q(t), and for all possible combinations of these measurement sets. In the current work we only present a selection of the results.
We need to stress that two-dimensional response surfaces provide only cross sections of the full six-dimensional parameter space. As such, these two-dimensionalsurfaces do not provide full proof about the uniquenessof the inverse problem. Nevertheless, the response surfaces are useful to study the behavior of the objectivefunction in the parameter space. The inverse parameter
estimation technique is expected to be unsuccessful if
response surfaces do not display a clearly defined global
minimum in the two-dimensional parameter planes (Šimůnek et al., 1998b).
First we studied parameter uniqueness with only ␪(z,t),h(z,t), or Q(t) in the objective function. Contour plotsof the response surfaces for ␣–n, ␣–Ks, and n–Ks areReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

562

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 4. Contour lines of the objective function ⌽ for the parameter combinations ␣–n, ␣–Ks, and n–Ks using either ␪(z,t ), h(z,t ), or Q(t ) datain the objective function.

shown in Fig. 4. Single minima are evident for ⌽(␣,n;␪)(Fig. 4a) and ⌽(n,Ks;␪) (Fig. 4g) but not for ⌽(␣,Ks;␪)(Fig. 4d). In contrast, identifiable minima are absent in
all response surfaces for h(z,t) and Q(t). From Fig. 4 itappears that ␪(z,t) data contain more useful information
for uniquely identifying ␣, n, and Ks than h(z,t) or Q(t)data. Similarly, response surfaces also showed that ␪rand ␪s can only be identified with ␪(z,t) data in theobjective function (comparison not shown). With solelyh(z,t) and Q(t) data in the objective function, identi-

fying water retention parameters is not possible because
there is no information about the volumetric water content of the soil. Note that some of the local minima in
Fig. 4 are artifacts resulting from the selected discretization of 50 by 50 points. The lower the resolution, themore artificial local minima will appear, especially nearparameter combinations where ⌽ changes only gradually, such as near the global minimum.
Parameter uniqueness with ␪(z,t) in the objectivefunction was investigated further by studying nine more

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

563

Fig. 5. Contour lines of the objective function ⌽ for nine parameter combinations using ␪(z,t ) data in the objective function.

response surfaces (Fig. 5). Response surfaces for ␪r and␪s combined with ␣, n, and Ks all exhibit clear identifiable
minima (Fig. 5a, 5b, 5d, 5e, 5g, and 5h). The response
surfaces for ␭ combined with ␣, n, and Ks are less conclusive because the shape of the contour plots is moreelongated and local minima appear to be present (Fig. 5c,
5f, and 5i). Whether ␪(z,t) alone will suffice to uniquely
define all soil hydraulic parameters in the inverse process can be doubted. The absence of well-defined minima in some of the two-dimensional response surfaces
and the divergent shape of some of the contour lines

in the surface plots indicate that ␪(z,t) alone may not
be sufficient.
The benefits of combining different measurement setsin the objective function are explored in Fig. 6. Response
surfaces for ␣–n, ␣–Ks, and n–Ks are given for ␪(z,t)
and Q(t) combined. Comparison with Fig. 4 shows thatthe contour lines of the response surfaces in Fig. 6 are
more convergent, indicating an increased potential fora unique solution during the inverse process. However,⌽(␣,Ks;␪,Q) still lacks a single minimum. Apparently,
different combinations of ␣ and Ks may result in approx-

564

VADOSE ZONE J., VOL. 4, AUGUST 2005

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

imately the same value of ⌽. Combining ␪(z,t) withh(z,t) and combining ␪(z,t) with both h(z,t) and Q(t) in
the objective function gave more or less similar response
surfaces as those shown in Fig. 6. In contrast, combining
h(z,t) and Q(t) data in the objective function resultedin ␣–n and ␣–Ks response surfaces that showed a longelongated region with low ⌽ values, indicating the absence of a minimum (results not shown).

Inverse Solutions
Different inverse procedures were tested for their ability
to retrieve the true soil hydraulic parameters. We used␪(z,t) and Q(t) data, both separately and combined, forfitting. Pressure head data were not considered for threereasons. First, pressure head data are not available forthe groundwater table lysimeter. Second, the measurement range of most tensiometers does not go below h ⫽⫺800 cm (Young and Sisson, 2002), making it impossible
to obtain a good data set in a depleted soil profile. Thisis especially true if the soil is fine textured like the silty
clay soil considered here. Third, the response surfaces
show that soil water content data are more valuable
than pressure head data for the current inverse problem.
Model grid and boundary conditions for the inverse
analysis were the same as for the forward simulation.
However, the initial conditions were specified in termsof water content ␪i(z) instead of pressure head hi(z),consistent with h(z,t) data not being available for thelysimeter. The LM minimization method was repeated
50 times for each case, using 50 different initial estimatesof the soil hydraulic parameters. Each initial parameter
set was generated by random selection from predetermined parameter intervals ([0–0.13] for ␪r, [0.3–0.6] for␪s, [0.0005–0.05 cm⫺1] for ␣, [1.05–2.0] for n, [0–25 cmd⫺1] for Ks, and [⫺5–10] for ␭). These intervals werealso used as bounds during the parameter optimizations.
The choice of 50 realizations was a compromise between
statistical considerations (more realizations give morereliable means) and the computational effort. We omitted LM minimizations that experienced numerical errors during the Hydrus model runs or that did not converge within 20 iterations.
The results of the parameter estimations for the hypothetical soil are summarized in Table 1 [␪(z,t) data],
Table 2 [Q(t) data], and Table 3 [␪(z,t) and Q(t) data].
Each table shows four different cases. In the first caseall six hydraulic parameters are optimized simultaneously. In the other three cases ␪r, ␭, and both ␪r and ␭
are fixed to their true values while the remaining parameters are optimized. Fixing parameters simplifies the minimization process by reducing the degrees of freedom
in the inverse process. We chose to fix ␪r and ␭ becauseof the relatively low sensitivity of ␪(z,t) and Q(t) tothese two parameters.
The coefficient of variability (CV, %) in Tables 1, 2,
and 3 relates to the variability in the optimized parameters:
Fig. 6. Contour lines of the objective function ⌽ for the parameter
combinations ␣–n, ␣–Ks, and n–Ks using ␪(z,t ) and Q(t ) data inthe objective function.

CV ⫽

100
冨bl 冨

冪 兺(bN ⫺⫺ 1b )
l

l

2

[12]

565

www.vadosezonejournal.org

Table 1. Results of the hydraulic parameter estimation for the hypothetical soil with volumetric water content ␪(z,t ) data in theobjective function.†
␪r

␪s

␣

n

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm⫺1
Mean
CV, %
NRMSE, %
␪r
Mean
CV, %
NRMSE, %
␭
Mean
CV, %
NRMSE, %␪r and ␭
MeanCV, %NRMSE, %True values

Ks

␭

SSQ␪

Success rate

cm d⫺1

0.118
23.2
32.4

0.497
1.8
2.1

0.023
40.3
80.9

1.409
13.5
15.9

15.00
72.6
463.1

2.771
156.8
554.9

4.00 ⫻ 100

0/20

0.101
–
–

0.491
2.0
2.0

0.019
44.0
63.0

1.320
9.5
9.5

15.58
70.7
477.5

4.968
94.0
732.5

5.95 ⫻ 100

0/23

0.107
35.0
37.4

0.492
0.3
0.3

0.014
11.6
12.5

1.402
8.8
11.2

2.99
55.6
49.9

⫺1.055––1.25 ⫻ 10⫺1

0/22

0.101
–
–
0.101

0.493
0.8
0.8
0.492

0.016
31.0
34.3
0.015

1.317
1.7
1.8
1.321

4.64
103.0
141.9
3.47

⫺1.055––⫺1.0555.07 ⫻ 10⫺3

17/35

–

–

† Number of attempted Levenberg–Marquardt minimizations is 50.

31.0%) and Ks (CV ⫽ 103.0%) with ␪r and ␭ fixed. The
uncertainty in the ␣ and Ks values is reflected in thehigh parameter correlation of between 0.95 and 0.98
in the converged minimizations (numbers not shown).
Note that ⫺1 stands for perfect negative correlation, 0
for no correlation, and ⫹1 for perfect positive correlation. The difficulty of identifying ␣ and Ks separately
was anticipated because the response surface ⌽(␣,Ks;␪)(Fig. 4d) lacked a clear minimum. Table 1 also showsthat failure to fix ␭ during the inverse process will leadto severe errors in the optimized value of Ks (NRMSE
of 463.1 and 477.5%).
None of the LM minimizations with solely Q(t) datain the objective function converged to the true soil hydraulic parameters (Table 2). The absence of a single
minimum in the response surfaces for ␣–n, ␣–Ks, andn–Ks (Fig. 4c, 4f, and 4i) already suggested that Q(t)data alone were insufficient to identify uniquely thehydraulic parameters. It is interesting to note that theabsence of water content data in the objective function
did not result in large errors in the optimized ␪s values(NRMSE 2.3–15.2%). It appears that the use of watercontent data for the initial condition helped in narrowing the range of possible values for this parameter.

where bl is the optimized value of the lth soil hydraulicparameter, bl is the average optimized value of the lthsoil hydraulic parameter, and N is the number of converged LM minimizations. The normalized root meansquare error (NRMSE, %) in Tables 1, 2, and 3 is ameasure of the variability of the optimized parameters
around the true parameter value:
NRMSE ⫽

100
冨b*l 冨冪兺 N ⫺ 1(bl ⫺ b*l )2[13]

where bl* is the true value of the lth soil hydraulic parameter. Both CV and NRMSE allow comparisons among
different parameters, irrespective of their invoked unit
or absolute value. The success rate in the last columnof Tables 1, 2, and 3 refers to the number of minimizations for which all optimized parameter values were
within 5% of their true values (number before the slash)as a fraction of the total number of converged minimizations N (number behind the slash).
Table 1 shows that fixing ␪r and ␭ is essential if only
␪(z,t) data are used in the objective function. Even then,
only 17 of 35 converged minimizations arrive at thecorrect soil hydraulic parameters. Also, considerable
variability remains in the optimized values of ␣ (CV ⫽Table 2. Results of the hydraulic parameter estimation for the hypothetical soil with cumulative bottom flux Q(t ) data in the objective function.†
␪r

␪s

␣

n

cm⫺1
Mean
CV, %
NRMSE, %
␪r
Mean
CV, %
NRMSE, %
␭
Mean
CV, %
NRMSE, %␪r and ␭
MeanCV, %NRMSE, %True values

Ks

␭

SSQQ

Success rate

cm d⫺1

0.087
50.4
45.7

0.492
15.2
15.2

0.024
58.8
113.9

1.281
11.7
11.8

11.05
55.1
288.5

1.046
354.7
408.9

2.40 ⫻ 103

0/11

0.101
–
–

0.485
6.7
6.8

0.018
50.5
66.8

1.196
5.3
11.1

11.83
70.6
350.9

1.465
187.1
362.8

7.79 ⫻ 102

0/9

0.059
70.2
59.5

0.477
5.4
6.1

0.017
17.6
25.7

1.214
7.6
11.0

11.76
66.5
335.1

⫺1.055––3.28 ⫻ 100

0/14

0.101
–
–
0.101

0.484
1.7
2.3
0.492

0.017
9.8
15.7
0.015

1.259
7.1
8.4
1.321

6.68
42.8
127.2
3.47

⫺1.055––⫺1.0556.90 ⫻ 10⫺2

0/11

–

–

† Number of attempted Levenberg–Marquardt minimizations is 50.

566

VADOSE ZONE J., VOL. 4, AUGUST 2005

Table 3. Results of the hydraulic parameter estimation for the hypothetical soil with volumetric water content ␪(z,t ) and cumulative
bottom flux Q(t ) data in the objective function.†
␪r

␪s

␣

n

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm⫺1
Mean
CV, %
NRMSE, %
␪r
Mean
CV, %
NRMSE, %
␭
Mean
CV, %
NRMSE, %␪r and ␭
MeanCV, %NRMSE, %True values

Ks

␭

SSQ␪

SSQQ

Success rate

cm d⫺1

0.107
33.7
36.2

0.521
6.4
9.1

0.021
31.4
59.3

1.353
12.9
13.5

11.08
81.8
346.0

1.104
333.2
408.1

2.70 ⫻ 100

1.78 ⫻ 102

0/15

0.101
–
–

0.507
6.1
7.1

0.018
47.3
60.6

1.393
14.9
16.7

8.9
90.5
282.0

1.146
279.6
372.2

1.67 ⫻ 100

1.05 ⫻ 102

0/17

0.095
33.6
32.2

0.492
0.4
0.4

0.015
11.3
11.3

1.341
6.8
7.1

3.58
42.1
43.7

⫺1.055––

1.56 ⫻ 10⫺2

1.04 ⫻ 100

1/22

0.101
–
–
0.101

0.492
0.0
0.0
0.492

0.015
0.9
0.9
0.015

1.322
0.2
0.2
1.321

3.43
2.3
2.5
3.47

⫺1.055––⫺1.0554.81 ⫻ 10⫺56.99 ⫻ 10⫺2

17/20

–

–

–

† Number of attempted Levenberg–Marquardt minimizations is 50.

The error in the optimized ␪r value was larger (NRMSE
45.7–59.5%) because the lowest ␪i value of 0.138 at z ⫽0 was higher than the true ␪r value of 0.101. Estimationof ␪r therefore required extrapolation beyond the measurement range. With pressure head instead of watercontent as the initial condition in the inverse analysis it
would have been impossible to approximate even ␪sbecause only the difference ␪s ⫺ ␪r could have been
optimized with only Q(t) data in the objective function.
Combining ␪(z,t) and Q(t) data in the objective function for the upward flow problem (Table 3) reduces thedeviations between the optimized and true values of Ks
and ␭ for all cases (judging from the NRMSE values).
Comparison between Tables 1 and 3 also shows that
the SSQ␪ is lowered by combining ␪(z,t) and Q(t) in
the objective function, signaling a better fit between themeasured and simulated water contents. However, ␪r
and ␭ still need to be fixed during the inverse process
to obtain the true soil hydraulic parameter values (17of 20 converged minimizations are successful). With ␪rand ␭ fixed, the variability in the optimized values of ␣and Ks is now small (NRMSE ⬍ 3%), despite the factthat the parameter correlation between these parameters remains high (≈0.97; numbers not shown).
The above results agree with findings of Šimůnek andvan Genuchten (1996), who, after studying downward
infiltration from a tension disc infiltrometer, concluded
that cumulative infiltration alone will not provide a
unique solution for the inverse problem. In subsequentstudies these authors suggested augmenting the cumulative infiltration data with measured final water contents(Šimůnek and van Genuchten, 1997) or with ␪(z,t) data
(Šimůnek et al., 1999) to obtain unique solutions. Šimůnek et al. (1999) also noted that augmenting cumulativeflux data with h(z,t) data did not improve the results.
It appears that the upward flux problem (this study)
and the downward infiltration study (Šimůnek and coworkers) pose similar challenges for inverse parameter
estimation. This is not a complete surprise since bothflow problems involve infiltration into a dry soil froma pressure head boundary condition.

Stability of the Inverse Solutions
The stability of the soil hydraulic parameter estimatesis examined by altering the setup of the inverse process
in four ways. The computational effort is limited to thecase where only ␪(z,t) data are used for fitting, andwhere ␪r and ␭ are being fixed (success rate of 17/35
in Table 1). This case is most relevant for the studied
groundwater table lysimeter, with no useful Q(t) orh(z,t) data available. First, the adequacy of using 50
replicates is demonstrated by increasing the number ofLM minimizations from 50 to 100. Second, the effect ofboth underestimating and overestimating the ␪r value
during the inverse analysis is quantified by fixing ␪r to
0.075 and 0.125, respectively (true ␪r ⫽ 0.101). Third,
we assessed the consequences of using ␭ ⫽ 0.5, a number
suggested by Mualem (1976) for most natural soils, instead of ␭ ⫽ ⫺1.055. Fourth, we tested the impact ofmeasurement errors in the water content data by adding
normally distributed random errors with zero mean and
0.0025 cm3 cm⫺3 standard deviation to the ␪(z,t) data.
The results of the stability analysis are summarizedin Table 4. Increasing the number of LM minimizations
from 50 to 100 did not result in significant changes inthe soil hydraulic parameter estimates (compare Tables
1 and 4). The biggest changes occurred in the mean Ks
value (increases from 4.64 to 4.98 cm d⫺1) and in the
NRMSE value for Ks (increases from 141.9 to 152.7%).
These increases are insignificant for a parameter like
Ks, which is often found to vary one or two orders ofmagnitude, even in homogeneous materials. This confirms that 50 replicates suffice to obtain reliable statistics
for the inverse solutions.
Underestimating or overestimating ␪r does not vitiate
the parameter estimates. The alterations in ␪r are offset
by small changes in the optimized mean n and Ks values.
The CV values for all parameters actually decrease,indicating that the optimizations have become less sensitive to the initial parameter estimates. Also, the SSQ␪values of 1.78 ⫻ 10⫺3 and 2.67 ⫻ 10⫺3 are lower thanthe SSQ␪ value of 5.07 ⫻ 10⫺3 found in Table 1. The
effect of fixing ␭ to 0.5 is somewhat detrimental. The

567

www.vadosezonejournal.org

Table 4. Stability of the hydraulic parameter estimates for the hypothetical soil with volumetric water content ␪(z,t ) data in theobjective function.†
␪r

␪s

␣

n

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm⫺1
Mean
CV, %
NRMSE, %

0.101
–
–

0.493
0.9
1.0

Mean
CV, %
NRMSE, %

0.075
–
–

0.493
0.5
0.6

Mean
CV, %
NRMSE, %

0.125
–
–

0.491
0.1
0.2

Mean
CV, %
NRMSE, %

0.101
–
–

0.499
2.6
3.1

Mean
CV, %
NRMSE, %True values

0.101
–
–
0.101

0.493
0.9
0.9
0.492

Ks

100 LM minimizations
1.315
4.98
1.9
101.8
1.9
152.7
␪r underestimated
0.016
1.280
5.12
21.3
1.1
71.2
24.4
3.3
115.7
␪r overestimated
0.014
1.373
2.58
1.8
0.2
3.8
6.2
4.0
26.3
␭ fixed to 0.5
0.020
1.470
8.84
44.4
8.7
94.1
67.1
15.0
286.6
Random errors in water content data
0.016
1.318
4.68
32.8
1.9
106.8
35.9
1.9
148.3
0.015
1.321
3.47
0.017
33.1
38.1

␭

SSQ␪

Success rate

⫺1.055––

6.38 ⫻ 10⫺3

35/66

⫺1.055––2.67 ⫻ 10⫺3

0/39

⫺1.055––

1.78 ⫻ 10⫺3

0/32

0.5
–
–

1.40 ⫻ 10⫺1

0/31

⫺1.055––⫺1.0551.25 ⫻ 10⫺2

5/37

–

–

cm d⫺1† Number of attempted Levenberg–Marquardt minimizations is 50, except for the case with 100 minimizations.

mean values of all four optimized parameters change
to offset the alteration in the ␭ value. The CV valuesfor ␪s, ␣, and n all increase, as well as the SSQ␪ value(to 1.40 ⫻ 10⫺1). However, the changes in the mean
parameter values remain relatively small and are unlikely to affect the calculated upward water flow in asignificant manner. For example, the calculated cumulative upward bottom flux changes only from 22.1 cm(using the true parameter values) to 22.9 cm with ␭ ⫽ 0.5.
The incorporation of random measurement errors in
the water content data does not significantly change theoutcome of the inverse solutions. The mean parametervalues remain about the same, the CV and NRMSE
values increase slightly, and the SSQ␪ values increase
to 1.25 ⫻ 10⫺2. Furthermore, the success rate is still 5/37,
despite the erroneous water contents. Note that a standard deviation for the water content measurement errorof 0.0025 cm3 cm⫺3 is representative of certain electromagnetic techniques like time domain reflectometry (e.g.,
Heimovaara and Bouten, 1990; Lambot et al., 2002). The
use of other, less consistent sampling techniques might
increase the standard deviation and thereby the uncertainty in the estimated soil hydraulic parameters. The
occurrence of persistent instrument errors might increase the uncertainty even more.

Implications
The response surfaces already indicated that the upward flow problem might suffer from nonuniquenessproblems, even when ␪(z,t) and Q(t) data are combined.
Different combinations of parameter values may leadto equally small values of the objective function. Thisis especially true for the ␣ and Ks parameters. The LM
minimization method will not always find the exact
global minimum under these circumstances. It is therefore not surprising that the inverse solutions yield a collection of parameter sets that perform about equally well.

We would like to stress that the failure to consistently
recover the true parameters during the inverse solutions
is not due to the LM algorithm itself. There is simply
no distinct global minimum in the objective function.
The use of a global search algorithm (e.g., Abbaspour
et al., 1997; Lambot et al., 2002), which examines thecomplete parameter space, would probably not help
under these circumstances. A global search algorithm
would only help if the LM minimization was frequently
ending up in a local minimum that is distinctly different
from the global minimum. Inspection of the individual
minimization results showed that this is seldom the case.
The nonuniqueness problem implies that we cannot
expect to find a unique set of soil hydraulic parameters
for the groundwater table lysimeter. We can only expect
to find a collection of parameter sets that describe thedata equally well. Analysis of the CV values will provide
important information about the reliability of the calculated mean parameter values. The stability analysis indicated that 50 LM minimization runs will suffice forthis purpose.
The failure to consistently find the true soil hydraulicparameters for the upward flow problem introduces uncertainty into the modeled soil water fluxes for the
groundwater table lysimeter. Calculated state variables
like water content and pressure head will be a functionof the chosen parameter set, which may be nonunique.
However, the use of nonunique parameter sets may
still yield valuable information on spatial and temporal
trends in the state variables and in the water balancecomponents.
MATERIALS AND METHODS
Experimental Setup
The groundwater table lysimeter (2 m width, 4 m long, and3 m deep) was installed in 1995 at the USDA-ARS facility inReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

568

VADOSE ZONE J., VOL. 4, AUGUST 2005

Parlier, CA. The top 1.7 m of the lysimeter consisted of an
undisturbed soil monolith from the west side of the San Joaquin Valley, CA. The soil was a saline silty clay (fine, smectitic,
thermic Sodic Haploxerert). The bottom 1.3 m of the lysimeterconsisted of disturbed soil from the same location, hand
packed to a dry bulk density of 1.3 to 1.35 g cm⫺3 (Schneider
et al., 1996). A truck scale measured changes in soil watercontent with time. A mariotte bottle was used to maintain the
groundwater table at about 1.0 m below soil surface. Overthe years, the lysimeter was planted with cotton (Gossypium
hirsutum L.), safflower (Carthamus tinctorius L.), and alfalfa(Medicago sativa L.), and was irrigated with good quality water(electrical conductivity, EC ⫽ 0.4 dS m⫺1) by subsurface drip,
surface drip, and sprinkler systems. The EC of the groundwater in the lysimeter was about 14 dS m⫺1.
All components of the lysimeter water balance were measured directly except evapotranspiration, which was calculatedby difference on an hourly basis. The depth of the groundwatertable was measured with an observation well in the center ofthe lysimeter using a pressure transducer. The distribution ofsoil water content with depth was monitored at two locationsin the lysimeter with multisensor EnviroSCAN capacitanceprobes (Sentek Pty Ltd., Kent Town, South Australia).1 The
capacitance probes each held 11 sensors at 10-cm intervals
starting at 5.5 cm (Location L1) and 7.0 cm (Location L2)
below the soil surface. The frequency readings of the capacitance sensors were converted into volumetric water contents
by the procedure of Kelleners et al. (2004). In this procedure,
the frequency was related to the soil permittivity by an equation that described the electromagnetic properties of the sensor–plastic access tube–soil system. The permittivity was thenrelated to the volumetric water content using the empiricalmodel of Malicki et al. (1996).

Model Setup
A summer fallow period (10 July 2001–26 Sept. 2001) was
selected for the inverse analysis. During this period the root
zone, which was depleted of water by the preceding safflower
crop, was gradually replenished by capillary rise from thegroundwater table. Daily rainfall (practically none) and reference evapotranspiration for the simulation period were takenfrom a California Irrigation Management Information Systemweather station located about 500 m from the lysimeter. The
reference evapotranspiration was converted to potential evaporation for bare soil using the dual crop coefficient procedureof Allen et al. (1998). For two 4-d periods (3–6 Aug. 2001 and
14–17 Sept. 2001), the lysimeter was covered with a plastic
sheet during sprinkler irrigation of the surrounding (fallow)
field. Potential evaporation during these times was set to zero.
The depth of the groundwater table as measured in the observation well varied between 88 and 100 cm below the soil
surface during the calculation period.
The safflower crop preceding the fallow period was drip
irrigated and planted in rows. This resulted in a horizontallyand vertically heterogeneous soil moisture pattern. The heterogeneous pattern persisted after the crop was harvested andthe irrigation was halted. During the inverse analysis of the
fallow period only one-dimensional vertical flow was assumed
for a single (x,y ) location. Thus, the water balance fluxes
measured with the lysimeter could not be used in the inverse
analysis because they represent an areal average for the complete lysimeter. In contrast, the capacitance probes, which repre1
The mention of trade or manufacturer names is made for information only and does not imply an endorsement, recommendation, or
exclusion by the USDA-ARS.

sent single locations, could be used, provided that the horizontal
soil water fluxes at these locations were negligible comparedwith the vertical (upward) fluxes. Alternatively, the lysimeter
could be described with a complete three-dimensional model.
This was not attempted because of the lack of spatial data (e.g.,
initial water contents for only two [x,y] locations) and becauseof the large computational requirements.
Thus, only water content data derived from the capacitancesensors were used in the inverse analysis. However, not all
capacitance data were included in the objective function. The
top two L1 sensors (at the 5.5- and 15.5-cm depths) and the
top three L2 sensors (at depths of 7 to 27 cm) were excluded
because the upward moving wetting front never reached thesesensors. In fact, these sensors recorded a decrease in the watercontent during the 79-d calculation period as a result of continued evaporation to the atmosphere. Elimination of the “drying” sensors enabled a pure estimate of the soil hydraulicproperties during wetting. These wetting soil hydraulic properties might be different from the drying soil hydraulic properties
because of hysteresis (e.g., Dane and Wierenga, 1975). The
problem of “drying” sensors did not occur for the hypotheticalsoil because of the extremely dry initial conditions in thatcase. The bottom L1 sensor (at the 105.5-cm depth) and thebottom L2 sensor (at the 107-cm depth) were excluded because they were always below the groundwater table andtherefore added little information for the parameter estimation process. The water content data used for fitting wererestricted to eight L1 sensors (25.5–95.5 cm depth) and seven
L2 sensors (37–97 cm depth).
Only the top 100 cm of the lysimeter was modeled with
Hydrus-1D. The soil profile was divided into 109 elementsthat varied in size from 0.1 cm at the soil surface to 1.0 cm
at depth. An atmospheric boundary condition was used at thetop, and a variable pressure head boundary was used at thebottom. The capacitance probe data from 9 July 2001 wereused to specify the initial ␪i(z ) condition. In the objectivefunction, simulated water contents were averaged over 6-cm
depth intervals to be consistent with the vertical zone of influence of the capacitance sensors (e.g., Paltineanu and Starr,1997). The value of hA was decreased from ⫺105 to ⫺106 cm
after trial calculations showed that maintaining hA ⫽ ⫺105
sometimes caused unrealistically high water contents in the
topsoil.
Note that a homogeneous soil profile was assumed in themodel while the lysimeter soil actually contained three horizons (silty clay Ap, 0–20 cm; silty clay B1, 20–75 cm; and clay
B2, 75–120 cm). This means that the optimized soil hydraulicparameters constitute a compromise between the hydraulicproperties of the B1 horizon and the B2 horizon (water content
data from shallow depths were not included in the objectivefunction). The hydraulic properties should therefore be viewed
as effective properties for the subsoil. In practice, the parameters will mainly represent the B1 horizon, since the hydraulicparameters are most sensitive to the water content changes
at shallow depths where the soil is initially relatively dry.
Simultaneous optimization of the soil hydraulic properties forthe individual soil horizons was not attempted because of
the nonuniqueness problems discussed above, and because
previous studies indicated that this requires not only ␪(z,t )data, but also h(z,t ) and Q(t ) data (e.g., Jacques et al., 2002;
Ritter et al., 2003).

RESULTS AND DISCUSSION
Inverse Solutions
The earlier results showed that the parameters ␪r and␭ need to be fixed to obtain unique soil hydraulic param-

569

www.vadosezonejournal.org

Table 5. Results of the soil hydraulic parameter estimation for the water content monitoring locations L1 and L2 in the groundwatertable lysimeter.†
␪r

␪s

␣

n

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

cm⫺1
L1
Mean
CV, %
L1
Mean
CV, %
L2
Mean
CV, %
L2
Mean
CV, %

Ks

␭

SSQ␪

Numberconverged

cm d⫺1

0.078
–

0.483
0.0

0.002
29.1

1.320
1.1

0.10
34.3

0.5
–

2.98 ⫻ 10⫺1

6

0.078
–

0.483
0.0

0.002
8.1

1.215
0.2

0.15
6.9

⫺1.055–2.99 ⫻ 10⫺1

2

0.078
–

0.487
0.0

0.004
13.0

1.241
0.6

0.26
17.4

0.5
–

1.31 ⫻ 10⫺1

11

0.078
–

0.487
0.0

0.004
7.5

1.178
0.2

0.32
9.0

⫺1.055–1.31 ⫻ 10⫺1

6‡

† Number of attempted Levenberg–Marquardt minimizations is 50.
‡ One solution was rejected because the minimization converged in to a local minimum (SSQ␪ of 1.90 as compared with an average SSQ␪ of 1.31 ⫻ 10⫺1).

eter sets. Even then, not all optimizations for the hypothetical soil arrived at the true parameter values; with␪(z,t) data in the objective function, only 17 out of 35
converged optimizations resulted in parameter sets thatwere ⬍5% from the true parameter values. The problem
now is to estimate ␪r and ␭ for the lysimeter soil. The
value of ␪r was fixed to 0.078, which is an estimate of thehygroscopic water content based on thermogravimetricmeasurements on undisturbed samples from a soil similar to that of the lysimeter (Kelleners et al., 2004). The
stability analysis for the hypothetical soil showed thatthe estimation of ␭ is more critical, so two different
values were selected. The value of ␭ was fixed to 0.5,
as suggested by Mualem (1976), and to ⫺1.055, as predicted by Rosetta. Each of the LM minimizations was
repeated 50 times with 50 initial parameter estimates
selected randomly from the same predetermined parameter intervals as for the hypothetical soil. These parameter intervals were also used as bounds during the minimization process.
The results of the parameter estimation for L1 andL2 in the lysimeter soil are shown in Table 5. The optimized ␪s and n values of 0.483 to 0.487 and 1.178 to
1.320, respectively, fall within the expected range for asilty clay soil. The low CV values for ␪s and n indicatethat these parameters were reliably identified. The optimized values for ␣ and Ks were 0.002 to 0.004 cm⫺1and 0.1 to 0.32 cm d⫺1, respectively. These values are
relatively low but not unrealistic. The elevated CV values for ␣ and Ks show that there is more uncertainty
in these parameters. Also, there is strong correlation
between ␣ and Ks (correlation coefficient ⬎ 0.99). Low
CV values for ␪s and n, high CV values for ␣ and Ks,
and high parameter correlation between ␣ and Ks werealso found for the hypothetical soil when ␪(z,t) data
were used for fitting (Table 1, ␪r and ␭ fixed).
The differences between the optimized soil hydraulicparameters for L1 and L2 are small. This is encouraging
because it indicates that the assumption of strictly vertical flow is not unreasonable. The water contents at thestart of the calculation period were different for L1 andL2. If significant horizontal fluxes would have occurredin the lysimeter, these different initial conditions would
probably have translated into different optimized soil

hydraulic parameters. Table 5 also shows that the valuefor ␭ has only a limited influence on the value of theother parameters and on the value of the objective function SSQ␪.
It is important to note that the converged LM minimizations generally end up at approximately the samevalues for ␣ and Ks, despite the high parameter correlation between these parameters. The sole exception is
one solution for L2 with ␭ ⫽ ⫺1.055 that yielded Ks ⫽17.52 cm d⫺1 and ␣ ⫽ 0.022 cm⫺1. However, the SSQ␪value of 1.90 for this solution was clearly out of linewith all other solutions, and was therefore not includedin Table 5. Apparently, the LM minimization ended up
in a local minimum in this case. It appears that highparameter correlation does not prevent the LM algorithm from finding the approximate solution, provided
multiple initial parameter estimates are used. However,parameter correlation slows down the minimization process in the sense that more iterations are required toarrive at a solution. This is probably the reason why
several of the attempted minimizations did not converge
within 20 iterations, both for the hypothetical soil andthe lysimeter soil.

Soil Hydraulic Properties
The soil hydraulic parameter estimates for ␭ ⫽ 0.5
were selected for further analysis. We preferred the ␭ ⫽0.5 estimates over the ␭ ⫽ ⫺1.055 estimates because of
the higher number of converged minimizations (Table5). The calculated soil hydraulic properties for L1 andL2 are depicted in Fig. 7, together with the Rosetta
prediction based on soil texture and dry bulk density
data. Two water retention curves and two hydraulicconductivity curves are calculated for each location, representing the upper and lower bounds according to thehydraulic parameter estimates.
The optimized (Hydrus) water retention curves point
toward a finer textured soil than the predicted (Rosetta)
curve. Also, the optimized hydraulic conductivity isroughly one order of magnitude lower than the predicted hydraulic conductivity for ⫺100 ⬍ h ⬍ 0. The
optimized soil hydraulic properties should only be used
for the conditions under which they were determined.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

570

VADOSE ZONE J., VOL. 4, AUGUST 2005

Fig. 7. Optimized (a) soil water retention and (b) hydraulic conductivity functions for the groundwater table lysimeter for Locations 1
and 2 compared with the Rosetta predicted soil hydraulic functions
for a silty clay soil.

Thus, the parameters can be used to study capillary risefrom the groundwater table. Flow problems involving
the drying of soil or flow problems involving downward
infiltration should not be studied using the present parameters. The first category demands a separate studyof the soil hydraulic properties during drying while thesecond category requires the inclusion of preferential
flow phenomena through macropores.

Measured vs. Calculated Water Content
Measured vs. calculated water contents with time forall depths included in the objective function are shownin Fig. 8. The effect of uncertainty on the soil hydraulicparameter estimates is assessed by showing both thelowest and highest calculated water contents at each
depth. The speed of the upward moving wetting front
is described reasonable well at both locations and does
not differ much for the range of parameter estimates.
However, the rate at which the water content increases
at individual depths is generally overestimated. Thismay be related to the measurement volume of the capacitance sensor, which was approximated by a vertical zoneof influence of 6 cm. Increasing or decreasing the zoneof influence in Hydrus did not result in significant improvements (results not shown). Since the measurement
volume of the sensor in the presence of a sharp wetting

Fig. 8. Measured (symbols) and calculated (lines) volumetric watercontent for the 79-d fallow period in the groundwater table lysimeter for (a) Location 1 and (b) Location 2. The numbers in the
figures refer to the depth below the soil surface of the center ofthe sensors.

front remains unknown, no further conclusions are warranted.
Figure 9 shows the measured and calculated watercontents with depth for the initial and final day of thecalculation period. Again, both the lowest and highest
calculated water contents are included and found to bepractically the same. Measured and calculated watercontents compare reasonably well, except in the topsoilat L2, where the measured values show a clear decreasein water content with time while the calculated values
hardly changed. This discrepancy could be due to several
factors. First, the “wetting” soil hydraulic properties
may be unsuitable to describe the drying process inthe topsoil due to the effect of hysteresis. Second, thehydraulic properties of the topsoil may differ from thehydraulic properties of the subsoil. Third, vapor flow,which is not included in Hydrus, could be the main
driving force behind the drying of the topsoil. Finally,
formation of cracks may result in a loss of contact between the soil and the plastic access tube of the capacitance probe, resulting in an underestimation of the water content. The true reason for the discrepancy in the
topsoil of L2 is probably due to a combination of thesefactors.

571

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

in Hydrus resulted in an overestimation of the increasein water content for the 100-cm profile. The 27-mmincrease in the water content of the entire lysimeteris considerably smaller than the capacitance-measured
increase in water content at L1 (112 mm) and L2 (50mm). Also, the cumulative upward bottom flux for thelysimeter (58 mm) was lower than the Hydrus-calculated
mean upward flux (125 and 93 mm for L1 and L2, respectively). This indicates that both L1 and L2 represented
relatively dry locations in the lysimeter at the start ofthe calculation period. This is probably due to the factthat both capacitance probes were located inside a crop
row where root water uptake from the preceding safflower crop was highest. If horizontal soil water fluxes
occurred in the lysimeter, they would have probably
added water to locations L1 and L2. This may haveresulted in an overestimation of the hydraulic conductivity in this study.

CONCLUSIONS

Fig. 9. Initial and final measured (symbols) and calculated (lines)
volumetric water content with depth for the 79-d fallow period in
the groundwater lysimeter for (a) Location 1 and (b) Location 2.

Water Balance
The cumulative water balances for L1 and L2 and forthe entire lysimeter are summarized in Table 6. Ideally,
the Hydrus-simulated change in water content shouldbe the same as the capacitance-measured change in water content at each location. This is not the case for L2,
where the underestimation of the drying of the topsoil
Table 6. Cumulative water balance terms for the summer fallow
period for the groundwater table lysimeter. Results for Hydrus
are averages using all converged soil hydraulic parameter estimates. Values in parentheses denote coefficient of variability(CV).
UpwardChange inPrecipitation Runoff Evaporation bottom flux water content

HydrusCap. probe

0

0 (0)

mm
L1
0 (4.6)

HydrusCap. probe
Lysimeter

0

0 (0)

L2
3 (5.9)93 (0.3)

0

0

31

58

125 (0.4)

125 (0.4)
112
91 (0.3)
50
27

Results for the hypothetical soil showed that the flow
variables ␪(z,t), h(z,t), and Q(t) for the upward flow
problem were most sensitive to the n parameter in the
VGM model and least sensitive to the ␪r and ␭ parameters. The response surfaces showed that the inclusionof ␪(z,t) data in the objective function is essential forparameter identifiability. With only h(z,t) or Q(t) datain the objective function, there is insufficient information about the water retention ␪(h) of the soil. The
inverse analysis showed that the ␪r and ␭ parameters
needed to be fixed during parameter estimation, evenif ␪(z,t) and Q(t) data were combined in the objectivefunction. Fixing the parameters that were least sensitive
to the flow variables decreased the degrees of freedom
in the parameter optimization in the most efficient manner. With only ␪(z,t) in the objective function, the LM
algorithm was capable of finding the correct soil hydraulicparameters, provided that ␪r and ␭ were fixed, and provided that multiple initial parameter estimates were used.
The hydraulic properties of the silty clay lysimeter
soil were determined for two (x,y) locations using ␪(z,t)
data and multiple initial parameter estimates. The ␪rparameter was fixed to 0.078, and the ␭ parameter was
fixed to ⫺1.055 or 0.5. The results of the parameter
estimation showed that consistent sets of soil hydraulicparameters could be determined for both locations. Differences in the hydraulic parameters for the two locations were limited. The variability in the optimized ␣and Ks values was higher than the variability in the ␪sand n values because of the high parameter correlation
between ␣ and Ks. The optimized parameters were notvery sensitive to the value of ␭. Measured and calculated
␪(z,t) values agreed reasonably well, although the model
generally overestimated the rate at which the water content increased at individual depths. The optimized soil
hydraulic properties can be used to study capillary risefrom the groundwater table.
REFERENCES
Abbaspour, K.C., M.Th. van Genuchten, R. Schulin, and E. Schläppi.
1997. A sequential uncertainty domain inverse procedure for esti-Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

572

VADOSE ZONE J., VOL. 4, AUGUST 2005

mating subsurface flow and transport parameters. Water Resour.
Res. 33:1879–1892.
Allen, R.G., L.S. Pereira, D. Raes, and M. Smith. 1998. Crop evapotranspiration, guidelines for computing crop water requirements.
FAO Irrig. and Drain. Pap. 56. FAO, Rome.
Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous
media. Hydrol. Pap. 3. Colorado State Univ., Fort Collins.
Dane, J.H., and P.J. Wierenga. 1975. Effect of hysteresis on the prediction of infiltration, redistribution and drainage of water in a layeredsoil. J. Hydrol. (Amsterdam) 25:229–242.
Gardner, W.R. 1958. Some steady-state solutions of the unsaturated
moisture flow equation with application to evaporation from awater table. Soil Sci. 85:228–232.
Heimovaara, T.J., and W. Bouten. 1990. A computer-controlled 36-channel time domain reflectometry system for monitoring soil watercontents. Water Resour. Res. 26:2311–2316.
Hopmans, J.W., and J. Šimůnek. 1997. Review of inverse estimationof soil hydraulic properties. p. 643–659. In M.Th. van Genuchten
et al. (ed.) Characterization and measurement of the hydraulicproperties of unsaturated porous media. Part 1. University of California, Riverside.
Hopmans, J.W., J. Šimůnek, N. Romano, and W. Durner. 2002. Inverse
modeling of transient water flow. p. 963–1008. In J.H. Dane andG.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser.
5. SSSA, Madison, WI.
Hudson, D.B., P.J. Wierenga, and R.G. Hills. 1996. Unsaturated hydraulic properties from upward flow into soil cores. Soil Sci. Soc.
Am. J. 60:388–396.
Jacques, D., J. Šimůnek, A. Timmerman, and J. Feyen. 2002. Calibration of Richards’ and convection-dispersion equations to field-scalewater flow and solute transport under rainfall conditions. J. Hydrol.
(Amsterdam) 259:15–31.
Kelleners, T.J., R.W.O. Soppe, J.E. Ayars, and T.H. Skaggs. 2004.
Calibration of capacitance probe sensors in a saline silty clay soil.
Soil Sci. Soc. Am. J. 68:770–778.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1987. Parameter
estimation for unsaturated flow and transport models, a review. J.
Hydrol. (Amsterdam) 91:255–293.
Kosugi, K. 1999. General model for unsaturated hydraulic conductivity
for soils with lognormal pore-size distribution. Soil Sci. Soc. Am.
J. 63:270–277.
Lambot, S., M. Javaux, F. Hupet, and M. Vanclooster. 2002. A global
multilevel coordinate search procedure for estimating the unsaturated soil hydraulic properties. Water Resour. Res. 38(11):1224.
doi:10.1029/2001WR001224.
Malicki, M.A., R. Plagge, and C.H. Roth. 1996. Improving the calibration of dielectric TDR soil moisture determination taking into
account the solid soil. Eur. J. Soil Sci. 47:357–366.
Marquardt, D.W. 1963. An algorithm for least-squares estimation of
nonlinear parameters. SIAM J. Appl. Math. 11:431–441.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
Neuman, S.P. 2003. Maximum likelihood Bayesian averaging of uncer-

tain model predictions. Stochastic Environ. Res. Risk Assess. 17:
291–305.
Paltineanu, I.C., and J.L. Starr. 1997. Real-time soil water dynamics
using multisensor capacitance probes: Laboratory calibration. Soil
Sci. Soc. Am. J. 61:1576–1585.
Ritter, A., F. Hupet, R. Muñoz-Carpena, S. Lambot, and M. Vanclooster. 2003. Using inverse methods for estimating soil hydraulic properties from field data as an alternative to direct methods. Agric.
Water Manage. 59:77–96.
Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated
hydraulic conductivity with the Mualem-van Genuchten model.
Soil Sci. Soc. Am. J. 64:843–851.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. Rosetta: A
computer program for estimating soil hydraulic parameters withhierarchical pedotransfer functions. J. Hydrol. (Amsterdam) 251:
163–176.
Schneider, A.D., J.E. Ayars, and C.J. Phene. 1996. Combining monolithic and repacked soil tanks for lysimeters from high water tablesites. Appl. Eng. Agric. 12:649–654.
Šimůnek, J., and J.W. Hopmans. 2002. Parameter optimization andnonlinear fitting. p. 139–157. In J.H. Dane, and G.C. Topp (ed.)
Methods of soil analysis. Part 4. SSSA, Madison, WI.
Šimůnek, J., K. Huang, M. Šejna, and M.Th. van Genuchten. 1998a.
The Hydrus-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variablysaturated media. Version 1.0. IGWMC-TPS-70. International
Ground Water Modeling Center, Colorado School of Mines,
Golden, CO.
Šimůnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated
soil hydraulic properties from tension disc infiltrometer data bynumerical inversion. Water Resour. Res. 32:2683–2696.
Šimůnek, J., and M.Th. van Genuchten. 1997. Estimating unsaturated
soil hydraulic properties from multiple tension disc infiltrometerdata. Soil Sci. 162:383–398.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1998b. Parameter estimation analysis of the evaporation method for determining
soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1999. Estimating unsaturated soil hydraulic properties from laboratory tension
disc infiltrometer experiments. Water Resour. Res. 35:2965–2979.
Šimůnek, J., O. Wendroth, N. Wypler, and M.Th. van Genuchten.
2001. Non-equilibrium water flow characterized by means of upward infiltration experiments. Eur. J. Soil Sci. 52:13–24.
Toorman, A.F., P.J. Wierenga, and R.G. Hills. 1992. Parameter estimation of hydraulic properties from one-step outflow data. Water
Resour. Res. 28:3021–3028.
van Genuchten, M.Th. 1980. A closed-form equation for predicting
the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
Young, M.H., A. Karagunduz, J. Šimůnek, and K.D. Pennell. 2002.
A modified upward infiltration method for characterizing soil hydraulic properties. Soil Sci. Soc. Am. J. 66:57–64.
Young, M.H., and J.B. Sisson. 2002. Tensiometry. p. 575–608. In J.H.
Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA,
Madison, WI.