3.6.2. Inverse Methods
JAN W. HOPMANS, University of California, Davis, CaliforniaJIÌÍ ŠIMæNEK, George E. Brown Jr. Salinity Laboratory, Riverside, California
NUNZIO ROMANO, University of Naples “Federico II”, Naples, Italy
WOLFGANG DURNER, Braunschweig Technical University, Braunschweig, Germany

3.6.2.1 Introduction
An adequate hydrological description of water flow and contaminant transport in
the vadose zone relies heavily on soil water retention and unsaturated hydraulic conductivity data of the considered spatial domain. The importance of an accurate soil
hydraulic description of the vadose zone, including the root zone, is increasingly
recognized in the fields of environmental engineering, soil physics and groundwater hydrology. With the current focus on the entire vadose zone, with increasing applications of watershed and land–atmosphere models, the spatial and temporalscales of interest have shifted to larger dimensions. This trend in increasing larger
spatial scales brings along with it the presence of increasing soil heterogeneity.
Hence, methodologies need to be available that allow for a rapid and accurate soil
hydraulic characterization, including its spatial variability.
Currently, many laboratory and field methods exist to determine the highly
nonlinear soil hydraulic functions in the vadose zone, represented by soil water retention and unsaturated hydraulic conductivity curves. Most methods require either
static or steady-state flow conditions to satisfy the assumptions of the corresponding analytical solution, which can make measurements time consuming. Excellent reviews of these types of direct methods are presented by Dirksen (1991), Reeve
and Carter (1991), and in Sections 3.3, 3.4, and 3.5 of this book. In contrast, the inverse modeling approach presented here estimates soil hydraulic properties fromtransient experiments, giving much more flexibility in experimental boundary conditions than required for steady-state methods. As an additional advantage, inverse
modeling allows the simultaneous estimation of both the soil water retention andunsaturated hydraulic conductivity function from a single transient experiment. In
other ways, inverse modeling of transient water flow is not much different than methods applied to steady flow. In either case, inversion of the governing equation is required to estimate the unsaturated hydraulic conductivity function from experimental
data. Whereas the steady-state methods invert Darcy’s equation, transient methods
invert Richards’ equation. Inversion of Richards’ equation requires repeated numerical simulation of the governing transient flow problem. Successful applicationof the inverse modeling technique improves both speed and accuracy, as there isno specific need to attain steady-state flow.
963

964

CHAPTER 3

By definition, inverse modeling is a general mathematical method to determine unknown causes on the basis of observation of their effects, as opposed to modeling of direct problems whose solution involves finding effects on the basis of a
description of their causes. Inverse modeling is widely accepted in engineering andphysical sciences for system characterization. For example, in Section 6.6, inverse
modeling is applied to estimate solute transport parameters. In this section, we focus
on the application of inverse modeling towards parameter estimation of soil hydraulic functions. Inverse modeling was first applied to the pressure plate outflow
method, where an initially saturated soil core was subjected to a series of step increases in air pressure with the drainage or outflow measured after each pressure
step increase (Gardner, 1956). Assuming a constant soil water diffusivity within pressure steps, the analytical solution yielded the soil water diffusivity as a function ofsoil water content. Doering (1965), using Gardner’s (1962) solution, simplified theoutflow method by proposing a one-step experiment, so that considerable timesaving
was achieved without loss in accuracy. Additional modifications were introduced
by Passioura (1976), Gupta et al. (1974), and Valiantzas and Kerkides (1990). This
last paper extended the outflow method to the simultaneous determination of thesoil water retention and unsaturated hydraulic conductivity functions using the
Brooks and Corey (1966) formulation of the soil hydraulic properties. Among thefirst to suggest the application of inverse modeling to estimate soil hydraulic parameters were Whisler and Watson (1968), who suggested the use of drainage volume and the soil water retention curve to estimate the unsaturated conductivity curveof a draining soil by matching observed and simulated drainage flows. Significant
progress has been made regarding the nonuniqueness of the optimized parameters
and the type of flow variable(s) to be included in the objective function. Zachman
et al. (1981) used simulated drainage data to show that the best parameter optimization results were obtained if cumulative drainage flow was matched with predicted outflow. Moreover, in a subsequent paper, Zachman et al. (1982) showed theimportance of selecting the correct parametric form for the soil hydraulic functions.
Kool and Parker (1988) and Russo et al. (1991) applied this approach to one-dimensional unsaturated flow problems. Their goal was to estimate parameter values for unsaturated porous media in laboratory soil columns. Kool and Parker(1988) discussed the advantage of including tensiometric data in the inverse approach, measured simultaneously with drainage rates for a hypothetical infiltration
and redistribution experiment.
The inverse modeling approach assumes a priori that the applied processmodel and the selected hydraulic relationships are an exact description of the soil’s
physical behavior, and therefore assumes that the model error is negligible. Thisimplies that deviations between simulation and observation are caused only by randomly distributed inaccuracies of measurements. Russo (1988) studied the influence of the parametric form of the soil hydraulic functions on the one-step outflow
optimization. For that purpose he investigated the Mualem–van Genuchten (vanGenuchten, 1980), Brooks and Corey (1966), and Gardner–Russo (Russo, 1988)
soil hydraulic models, using data from Kool et al. (1985) and Parker et al. (1985).
In a subsequent paper, Russo et al. (1991) pointed out that the larger number of parameters in the van Genuchten model might enhance the likelihood of nonuniquenessand instability in the inverse solution. Also, Chen et al. (1999) tested seven differ-

THE SOIL SOLUTION PHASE

965

ent soil hydraulic models in their ability to fit multistep outflow experimental data,and concluded that only four of these models were able to describe the outflow data
successfully.
Although laboratory experiments have the advantage of being quick andprecise, they often lead to soil hydraulic properties that are not representative of the
field. Among the benefits of inverse methods is the fact that they are equally applicable to field experiments, even under nontrivial boundary conditions. The first
application of the inverse problem to field data was conducted by Dane and Hruska(1983), in which the parameters of van Genuchten’s (1980) soil hydraulic functions
were optimized from transient drainage data. An in-depth discussion of the application of the inverse parameter estimation method as applied to field experiments
is presented by Kool et al. (1987) and Kool and Parker (1988).
The objective of this section is to introduce the reader to the theory and applications of inverse modeling for parameter estimation of soil hydraulic functions,using various proven experimental techniques. New developments are continuously
made, so this discussion must be regarded as “work in progress”. In any case, inverse modeling in general has been shown to be an exciting new tool that allows
for soil hydraulic characterization using a wide spectrum of transient laboratory andfield experiments.
In the subsequent sections we summarize the recent developments in various
laboratory and field techniques to indirectly determine soil hydraulic properties using
the inverse method. First, we introduce the reader to the basic theoretical concepts
underlying all presented applications of inverse modeling for hydraulic parameter
estimation. For a better understanding of this theoretical background, we recommend to first read the general review of the application and theory of parameter optimization in Section 1.7. Subsequently, we discuss specific established experimental
methods for hydraulic parameter estimation. The experiments by themselves are reviewed elsewhere (e.g., Sections 3.3, 3.5, and 3.6.1); however, the type and number of measurements may need to be adapted to better serve the inverse methodology. Accordingly, issues that relate to the analysis and interpretation of the parameter
estimation results (both simulation and optimization modeling) will be specificallyaddressed. Finally, we will introduce a step-by-step example of the inverse method,
analyzing the influence of experimental conditions on parameter optimization results.
3.6.2.2 Theory of Flow and Optimization
The inverse method includes three interrelated functional parts (Fig. 3.6.2–1):(i) a controlled transient flow experiment for which boundary and initial conditions
are prescribed and various flow variables are measured, such as cumulative infiltration and/or drainage cumulative and/or matric head and/or water content; (ii) a
numerical flow model simulating the transient flow regime of this experiment, using
initial estimates of the parametric soil hydraulic functions; and (iii) an optimization algorithm, which estimates the unknown parameters through minimization of
the difference between observed and simulated flow variables (residuals) definedin an objective function (φ) through an iterative solution of the transient flow equation. The quality of the final solution of the parameter estimation problem is de-

966

CHAPTER 3

Fig. 3.6.2–1. Flow chart of the inverse method illustrating the integration of measurement, modeling,
and optimization (modified from Hopmans & Simçnek, 1999).

pendent on each of these three individual components as well as their integration
within a computational framework. The three components are interfaced throughdata files that include experimental, numerical flow model, and parameter optimization results. Parameters of the soil hydraulic functions are updated iteratively
in the optimization routine, thereby continuously reducing the residuals until a predetermined convergence criterion (e.g., reduction in φ value between two consecutive iterations) has been achieved. As any other soil hydraulic characterization
method, the magnitude of the residuals in the objective function will depend on theproper selection of the soil hydraulic model; that is, an inappropriate model will
increase deviations between measurements and simulations. In this section we describe the water flow simulation and parameter optimization components, andelaborate on the consequences of implied assumptions.
3.6.2.2.a Water Flow Modeling
One-dimensional vertical variably saturated water flow in soils is usually described using Richards’ equation:∂θ/∂t = (∂/∂z)[K(hm)(∂hm/∂z) + K(hm)][3.6.2–1]

where θ denotes the volumetric water content (L3 L−3), hm is the soil water matrichead (L), K is the unsaturated hydraulic conductivity (L T−1), z is a vertical coordinate (L), positive upwards, and t is time (T). Equation [3.6.2–1] can be solved nu-

THE SOIL SOLUTION PHASE

967

merically using either the standard approximation of the hm-based formulation byreplacing ∂θ/∂t by C(hm)∂hm/∂t, where C(hm) = dθ/dhm denotes the water capacity
function (L−1), or the mixed formulation represented by Eq. [3.6.2–1] (Celia et al.,1990). Although Eq. [3.6.2–1] is generally accepted, its derivation for unsaturated
macroscopic water flow, including the formulation of the effective hydraulic properties, [i.e., the soil water retention function, θ(hm), and unsaturated hydraulic conductivity function, K(hm)], includes many assumptions. Most importantly these are:
1. Flow is isothermal, thereby neglecting possible temperature effects on thesoil hydraulic properties and temperature-induced water flow.
2. The influence of the air phase on water flow is negligible.
3. The relation between water content and water energy status for a monotonic drainage or wetting process is unique, that is, independent of the waterflow rate.
These assumptions are not necessarily fulfilled for transient experiments employed
in inverse modeling. Therefore, if deemed appropriate, other flow models may besubstituted for Eq. [3.6.2–1]. For example, this was demonstrated by Schultze etal. (1999) and Chen et al. (1999), who used two-fluid phase flow models to simulate unsaturated water flow. To further elaborate on the consequences of these assumptions for simulated water flow is beyond the scope of this study, and we refer
the reader to Bear (1972), Nielsen et al. (1986), Quintard and Whitaker (1999), andvan Genuchten and Sudicky (1999) for an in-depth analysis.
Some of the presented experimental methods require numerical solution of
the two-dimensional flow equation, which for radially symmetric Darcian flow isassumed to be described by:∂θ/∂t = (1/r)(∂/∂r)[rK(∂hm/∂r)] + (∂/∂z)[K(∂hm/∂z)] + (∂K/∂z)[3.6.2–2]

where r is the radial coordinate (L). The applications of Eq. [3.6.2–2], presented
below, assume that the soil is uniform and isotropic, and that a multidimensional
approach of the problem is required to allow for spatially variable boundary conditions. Otherwise, the same assumptions as for Eq. [3.6.2–1] are applicable.
The soil water retention curve and unsaturated hydraulic conductivity function can be described by various expressions, as presented in Section 3.3.4. However, for convenience, here we use the soil water retention function of van Genuchten(1980)
Se(hm) = (θ − θr)/(θs − θr) = 1/(1 + |αhm|n)m[3.6.2–3]

which after substitution in the capillary model of Mualem (1976) yields the following hydraulic conductivity model
K(hm) = KsSel[1 − (1 − Se1/m)m]2[3.6.2–4]

In Eq. [3.6.2–3] and [3.6.2–4], θr and θs denote the residual and saturated volumetricwater contents (L3 L−3), respectively; Se is the effective saturation, Ks is the saturated hydraulic conductivity (L T−1), l is a pore connectivity coefficient, and α (L−1),

968

CHAPTER 3n, and m = 1 − 1/n are empirical coefficients. Most hydraulic models assume thatthe predictive K(θ) is coupled with soil water retention relations through the use ofcommon parameters, as in Eq. [3.6.2–4]. However, the general inverse modelingapproach does not require the coupling of the two soil hydraulic functions, and recent studies suggest that improved parameter optimization results are obtained bytheir decoupling (Vogel et al., 1999; Durner et al., 1999a). Moreover, it has beenshown that increased flexibility in the representation of the soil water retention curve,
for example, by using a bimodal pore-size distribution model (Durner, 1994; Zurmühl & Durner, 1998) or spline functions (Kulkarni et al., 1998), may improve parameter optimization results. Typically, the number of parameters needed to describe
the soil hydraulic functions varies between four and seven. Some parameters canbe measured independently (e.g., θs), assigned an assumed value (e.g., the tortuosity parameter, l), or related, such as n and m in the van Genuchten (1980) expression.
3.6.2.2.b Parameter Optimization
Objective Function—How Can the Optimal Parameter Combination
Be Determined? The following is a summary of the various analyses that are routinely used in the theory and application of inverse modeling. For a more in-depth
treatise, we refer the reader to Hopmans and Šimçnek (1999) and Section 1.7.
Desired hydraulic parameters are determined by systematically minimizing
the differences between observed and simulated state variables. The total of thesedifferences is expressed by an objective function, φ, which may be defined as
j=my

i=nj

β,y) = Σ
β)]2
φ(β
vj Σ
wi,j[yj*(z,ti) − yj(z,ti,βj=1
i=1[3.6.2–5]

where the right-hand side represents the residuals between the measured (yj*) andcorresponding model-predicted (yj) space–time variables using the soil hydraulicparameters of the optimized parameter vector, β. The first summation sign sums
the residual for all measurement types my, whereas the variable nj in the second summation denotes the number of measurements for a certain measurement type j. Typically, in water flow studies, yj* may represent water flux density, cumulative water
flow, soil water matric head, or soil water content values. Assuming that the measurement errors within a measurement type are independent and uncorrelated, aweighted least-squares problem represents the maximum likelihood estimator (Section 1.7). Weighting factor values for vj can be selected such that data types areweighted equally using a normalization procedure or such that they are equal to thereciprocal of the measurement variance of measurement type j (Section 1.7); additional weighting (wi,j) can be assigned to individual data (Hollenbeck et al.,2000). In addition to the transient measurements, the objective function can also
include independently measured soil water retention or unsaturated hydraulic conductivity data points.
Among the various nonlinear optimization techniques available to minimize
the objective function (Section 1.7), the Levenberg–Marquardt method is the mostTHE SOIL SOLUTION PHASE

969

widely used. This method combines the Newton method with the steepest descend
method, providing confidence intervals for the optimized parameters.
Ill-Posedness—Can an Optimum Solution Be Found? A solution of theinverse problem is attained by minimization of Eq. [3.6.2–5]. To determine whether
the inverse problem is at all solvable, it must be “correctly posed”. An incorrectly
posed, or ill-posed, inverse solution causes nonuniqueness or divergent results. It
is generally characterized by nonuniqueness, and/or nonidentifiability, and/or instability of the identified parameters. Instability stems from the fact that small errors in the measured variable may result in large changes of the optimized parameters. Nonuniqueness occurs when a given response leads to more than one set of
optimized parameters, β. Generally, the degree of nonuniqueness will increase withan increase in measurement errors of the measured data in φ. Nonidentifiability occurs if more than a single parameter set leads to the same model response. If a parameter set is nonidentifiable, the problem is certainly nonunique and therefore illposed. Identifiability can be increased, and hence nonuniqueness reduced, bydecreasing the number of parameters to be optimized. Nonuniqueness can also be
caused by a lack of sensitivity of the flow variables to certain parameter combinations. Accordingly, nonuniqueness depends on the type of measured data, appliedvalues of the weighting factors in Eq. [3.6.2–5], and the suitability of the boundary conditions for the specific flow experiment. Moreover, sensitivity is influenced
by the type and number of optimized parameters, and by model and input measurement errors. An experiment must be designed such that direct information isavailable for the least sensitive parameters, thereby eliminating them from the parameter set, and by providing well-constrained initial estimates. It is generally recommended to test for nonuniqueness by solving the inverse problem repeatedly using
different initial parameter estimates.
Response Surface Analysis—How to Visualize Uniqueness of a Solution?
The behavior of the inverse problem can be evaluated by plotting the value of theobjective function (φ), against pairs of optimized parameters to obtain response surfaces. Each response surface is obtained by solving the flow equation, with the appropriate boundary and initial conditions, for many possible combinations of a selected pair of parameter values within a predetermined range, while keeping the otherparameters constant. This type of analysis should also be done before the experiment is conducted to investigate the well-posedness of the inverse problem. Since
the influence of only two parameters on φ is shown in a response surface, each surface represents only a cross section of the full parameter space. The behavior of φacross the complete parameter space requires calculating response surfaces for all
possible pairs of parameter combinations. Since only two parameters appear within
a single response surface analysis, the behavior of φ in the various parameter planes
is only an indication of the uniqueness of the solution in the full parameter space.
For example, local minima may be present, but not appear in any of the cross-sectional planes (Šimçnek & van Genuchten, 1996). On the other hand, response surfaces can reveal the occurrence of local minima, the presence of a well-defined
global minimum, and can assist in evaluating parameter sensitivity and correlation.
Since the shape of the response surfaces will depend on the measured variables and970

CHAPTER 3

the applied experimental boundary conditions, response surface analysis is essential in the evaluation of optimal experimental designs.
Parameter Uncertainty—How Precise Are the Parameter Estimates? The
precision of the estimated parameters can be assessed by uncertainty analysis,
which yields confidence intervals for the estimated parameters at the φ minimum.
This analysis assumes that uncertainty is caused by measurement errors only, thatis, that the model error is zero, and that the inverse solution has converged to theglobal minimum. Under these conditions, the residuals in φ are independent and normally distributed. Uncertainty analysis further assumes that any variable, yj, in Eq.
[3.6.2–5] is a linear function of the optimized parameters within their computed confidence interval. Although restrictive and only approximately valid for nonlinear
problems, the uncertainty analysis provides a means to compare confidence intervals among parameters; thereby indicating which parameters should be measured
or estimated independently. For the stated assumptions, parameter standard deviation values can be determined and confidence intervals can be estimated. However, these confidence interval estimates are not necessarily correct if parameters
are correlated. In addition, correlated parameters can cause a slow convergence rate,
and will increase nonuniqueness and parameter uncertainty. Hollenbeck et al.
(2000) pointed out that testing of the hydraulic model’s adequacy is required first
before parameter uncertainty analysis is conducted. Clearly, uncertainty analysis
of optimized hydraulic parameters is not meaningful if the tested model fails theadequacy test, for example, because of experimental or modeling errors. For that
purpose, a correlation analysis between measured and optimized flow variables only
is not sufficient, but should include a residual analysis, obtained from a plot of observed and simulated data vs. the independent variable (space and/or time). Modeladequacy can be tested statistically by comparing the average residual with a predefined (expected) measurement error, as estimated from experience or sensor
specifications. In Eq. [3.6.2–5], inadequate hydraulic models yield φ values that arelarger than one, if vj = 1/njσj2 (σj2 is measurement error variance) and wi,j = 1.
Sensitivity Analysis—Is the Experimental Design Suitable? Parameter sensitivity is determined by the derivative of the objective function with respect to aparticular parameter. A flow experiment designed for inverse modeling should include measurements that are most sensitive to changes in the optimized parameters, thereby eliminating ill-posedness. The sensitivity analysis of the optimization
problem depends on type, location, and frequency of measurements; thereby providing information on the optimal experimental design. Moreover, a higher sensitivity will result in quicker convergence of the parameter estimation problem. Sensitivity coefficients can be computed a priori from hypothetical experiments, using
the Jacobian, as defined in Eq. [1.7–15] and [1.7–17] in Section 1.7. Large sensitivity values indicate a well-defined minimum with precise parameter estimates. Sensitivity and parameter uncertainty are inversely related; that is, highly sensitive measurement variables yield small parameter uncertainties for a given measurement error.
Since the Jacobian matrix is determined for the optimum parameter set only, thesensitivity coefficients characterize the behavior of φ at its minimum, and do notdescribe problem sensitivity anywhere else in the parameter space. Sensitivity coefficients can be calculated as a function of time (Inoue et al., 1998), and can beTHE SOIL SOLUTION PHASE

971

used to determine optimum sensor location within the spatial domain of the experiment. Moreover, as was reported by Zurmühl and Durner (1998), a similar shapeof plots with parameter sensitivity coefficient vs. time indicates correlated parameters.
In summary, it must be pointed out that the success of inverse modeling depends on the suitability and quality of each of the three earlier mentioned components, namely, (i) the experimental design, that is, the choice of boundary conditions, the location and time resolution of the measurement sensors, and the degreeof accuracy of the experimental data; (ii) the suitability of the transient flow modeland hydraulic functions; and (iii) the robustness of the optimization algorithm asdetermined by its convergence to a global minimum. If any of these three components is unsatisfactory, inverse modeling may diverge (fail to find a global minimum), or result in an ill-posed problem yielding inaccurate parameters with highuncertainties. Further requirements for well-posed problems will be illustrated in
the example in Section 3.6.2–8.
3.6.2.3 Multistep Outflow Method
3.6.2.3.a Introduction
Kool et al. (1985) were among the first to apply the inverse approach by numerical solution of the Richards equation for a one-step outflow experiment. They
concluded that uniqueness problems are minimized if the experiment is designed
to cover a wide water content range. Moreover, they also determined that initial parameter estimates must be reasonably close to their true values and that outflow
measurement errors must be small. Parker et al. (1985) subsequently experimentally applied the one-step method to four different soils of different texture and concluded that θ(hm) and K(θ) can be optimized simultaneously by using cumulative
outflow as a function of time. It was also found that the optimized soil hydraulicfunctions could be extrapolated to a water content range beyond that achieved withthe single pressure step by including in the objective function an independently
measured point on the soil water retention curve.
The need for independently measured soil water retention data in the one-stepoutflow optimization procedure was demonstrated by van Dam et al. (1992). In theirstudy, measured unsaturated hydraulic conductivity data were compared with optimized hydraulic functions using the one-step outflow procedure with and without measured θ(hm) data. They concluded that optimization using measured outflow data alone is inadequate and that additional retention data are necessary toobtain accurate K(θ) functions. Similarly, Toorman et al. (1992) concluded thatuniqueness problems in the transient one-step outflow experiment were minimized
if matric head data were included in the objective function.
To circumvent the need for additional soil water pressure measurements inthe outflow experiment, van Dam et al. (1994) conducted outflow experiments inwhich the pneumatic pressure was increased in several smaller steps. Their work,
using a loam soil, showed that the outflow data of a multistep experiment contained
sufficient information for unique estimates of the soil hydraulic functions. The experimental work by Eching and Hopmans (1993a, b) and Eching et al. (1994)
showed how the multistep method, when combined with automated matric head

972

CHAPTER 3

measurements during drainage of the soil core, resulted in unique parameter values for the optimized soil hydraulic functions for four different textured soils. In
their analysis, an excellent match between optimized and independently measured
soil water retention and unsaturated hydraulic conductivity data was found. The multistep outflow method was also recommended by Durner et al. (1999b) after they
compared it with the classical one-step method, although they also showed that theperformance of the one-step method depends largely on the pressure step size relative to the shape of the soil water retention function. Their work further demonstrated that a reliable estimation of the unsaturated hydraulic conductivity function
is dependent on the accuracy of the retention model, if the coupled approach of Eq.
[3.6.2–4] is used. They concluded that combination of the multistep method with
outflow and tensiometric data in the objective function would yield accurate estimations of both the retention and the hydraulic conductivity parameters for a wide
range of soil textures.
3.6.2.3.b Experimental Procedures
The experimental apparatus (Fig. 3.6.2–2) is based on the pressure cellmethod, presented in Section 3.3.2, and can be adapted (Eching et al., 1993b) toinclude one or two miniature tensiometers connected to pressure transducers(Wildenschild et al., 2001). The porous membrane can be either a ceramic plate or
porous nylon (MSI Inc., Westborough, MA) with a sufficiently high air-entry value.
The advantage of the porous nylon is its low resistance, thereby minimizing pressure differences across the membrane. However, the thin membrane accelerates thepossible diffusion of soil gas across the membrane, thereby causing accumulationof gas beneath the membrane. Gas accumulation, caused by exsolution of air fromthe water phase can be eliminated by using pressurized N2 gas instead of air. Continuous and unattended monitoring of drainage flow rates can be easily accomplishedby installation of a pressure transducer in the bottom of a water-receiving burette.
Details of various experimental setups can be found in Eching and Hopmans(1993a, b) and Zurmühl (1998).
After assembling the soil core between the end plates and insertion of the tensiometer(s), the soil is saturated with a CaCl2 solution (0.0001–0.01 M, depending
on soil type and Na content). If needed, AgNO3 solution, mercuric chloride and/or
thymol (see Section 3.3.2.1.d) should be added to the saturating solution to prevent
microbial clogging of the porous membrane and/or soil pores. After saturation, thepressure transducers can be calibrated from known pressures in the soil water andthe burette. Subsequently, the soil is slightly unsaturated by exceeding its air-entry
value across to height of the sample to assure the presence of a continuous gas phaseat the start of the multistep experiment. As was demonstrated by Hopmans et al.
(1992), this is needed to obtain a correct solution of the Richards equation, by ensuring the availability of air to displace the water phase during the drainage experiment.
After hydraulic equilibrium has been established, the first pressure increment
is applied and the soil matric head and drainage volume are continuously monitored
and logged using the pressure transducers. The volumetric water content at the conclusion of the outflow experiment is determined by oven-drying the soil core. Sat-

THE SOIL SOLUTION PHASE

973

Fig. 3.6.2–2. Schematic of the apparatus for outflow experiments (modified from Wildenschild et al.,2001).

urated and initial water content values, in addition to intermediate volumetric watercontent values, are back calculated using this oven-dry value and the measured cumulative outflow volumes between pressure steps. The multistep outflow techniquewill be improved by including the following considerations:
1. The choice of the number and magnitude of the pressure increments is notknown a priori. In the limiting case, one may select a single pressure increment (onestep outflow experiment). However, it has been clearly demonstrated that selectionof multiple steps increases the sensitivity of the measurements with regard to theparameters to be optimized. We recommend using available information on soil
water retention to select the number and size of the pressure increments, thereby
making sure that the general shape of the retention curve is preserved. This meansthat pressure differences should be small near the air-entry point, where the retention curve is steep, but can be large for the dry range, where the retention curve is974

CHAPTER 3

Tables 3.6.2–1. Thickness (d), saturated hydraulic conductivity (Ks), resistance (R), and air-entry pressure (ha,e) for a ceramic plate, sintered glass plates (Robu, Hattingen, Germany), and a nylon membrane (Osmonics Laboratory Products; www. osmolabstore.com).

Ceramic plate
Sintered glass plate P5
Sintered glass plate P4
Sintered glass plate P3
Nylon membrane

d

Ks

R

ha,e

cm

cm h−1

h

kPa

0.74
0.70
0.70
0.70
0.01

0.0498
0.1
2.1
8.5
0.025

14.86
7
0.33
0.08
0.4

100
45
15
7
170

nearly flat. Ideally, one selects the pressure increments such that about equal
drainage volumes are obtained for each increment.
2. Questions are often asked about the duration of pressure increments. There
are various issues that must be considered. First, it has been shown (e.g., Schultze
et al., 1999) that one of the most important advantages of the inverse method is thatexcellent results can be obtained without waiting for hydraulic equilibrium between
pressure steps. Moreover, true hydraulic equilibrium is often not attained in the drysoil moisture range, because of the low unsaturated hydraulic conductivity valueof the soil core, especially near the porous membrane. As the duration of the outflow experiment increases, there is increasing danger of air accumulation under the
porous membrane by gaseous diffusion and clogging of the membrane and soil bymicrobial growth. If no tensiometric measurements are taken, however, one might
prefer increasing the gas pressure only after there is no measurable outflow, whichmakes the method equivalent to the one presented in Section 3.3.2.
3. To increase the sensitivity of the parameter optimization procedure, it isimportant to minimize the influence of the porous membrane on drainage rate. For
example, the resistance of the porous membrane must be low relative to that of theunsaturated soil. This might not be the case when using thick porous ceramic
plates for near-saturated, coarse-textured soils. Therefore, it is recommended to use
a thin nylon porous membrane with a large air-entry value (Table 3.6.2–1). The low
resistance nylon membrane creates only minor differences in water pressure acrossthe membrane, and is now routinely used in outflow experiments. Moreover, thelow resistance nylon membrane eliminates the need for a separate calculation of thematric head at the soil–membrane interface for direct conductivity calculations.
However, care must be taken to maintain hydraulic contact between the membrane
and the soil sample. Especially for sandy soils, highly conductive sintered glass
plates with limited air-entry pressures are a viable alternative (Table 3.6.2–1). If theresistance of the porous plate or nylon membrane is not negligible, it must be included in the numerical simulation of the experiment, effectively requiring simulation for a two-layered domain. Under these circumstances, the porous membranein the flow domain must be replaced by an equivalent porous medium with the same
total resistance. In any case, it is good practice to measure the plate resistance or
its saturated hydraulic conductivity separately, before and after the experiment, toassess whether clogging has occurred due to microbial growth or migration of soil
particles.
4. A well-posed inverse problem is achieved if the matric head is measured
in the soil core during the outflow experiment. A single tensiometer is best placed

THE SOIL SOLUTION PHASE

975

near the center of the soil core. However, insertion of a second tensiometer close
to the bottom outflow boundary is recommended, since it provides valuable information about the validity of the pressure continuity assumption across the porous
membrane.
5. In the dry range, it is often found that, although outflow has essentiallyceased, the soil water matric head is not in equilibrium with the water beneath theporous membrane. That is, at low water content values, the corresponding hydraulicconductivity values are too low for drainage to occur, thereby preventing the soil
from attaining hydraulic equilibrium. This is especially the case for coarse-textured
soils and has wide implications for soil water retention measurements in general(see Section 3.3.1).
6. As static equilibrium conditions are not required in the inverse analysis,
there is additional freedom in the selection of the physical dimensions of the soilcore. However, column length affects parameter sensitivity in two ways. It is expected that the sensitivity for the hydraulic conductivity parameters increases withincreasing core length because of the resulting increasing flow rates. Alternatively,
the sensitivity of the soil water retention parameters increases as the soil core
length decreases. The practical core length ranges from 5 to 20 cm. The core diameter has no effect on parameter sensitivity; however, collection of larger-diameter cores reduces among-core variability in undisturbed soils.
7. The general recommendation to start a transient outflow experiment at nearsaturation by exceeding the soil’s air-entry value first will reduce the sensitivity ofthe saturated hydraulic conductivity on the inverse solution. Moreover, an initialdrainage experiment on the same core is required to define the soil’s air-entry value.
The near-saturation requirement is not needed if the initial pressure step, startingfrom saturation, is slightly larger than the air-entry value of the soil.

3.6.2.3.c Simulations and Optimization
Initial and boundary conditions applicable to the multistep experiment arehm(z,t) = hm,i(z)t = 0, 0 < z < L

q(z,t) = 0

t > 0, z = L

hm(z,t) = h(z,t) − hat > 0, z = 0[3.6.2–6]

where hm,i is the initial matric head, q denotes the flux density (L T−1), z = 0 is thebottom of the porous membrane or plate, z = L is the top of the soil core, h(0,t) is
the water pressure head at the bottom of the porous membrane, and ha is either thepneumatic gas pressure applied to the top of the soil core (z = L) or the suction applied beneath the porous membrane (z = 0). For example, if h(0,t) = 1 cm (heightof water above the bottom of the porous membrane), and ha = 80 cm (applied to the
top of the soil core), then hm = −79 cm. Similarly, for the same water level in theburette, if an 80-cm suction head is applied to the water in the burette, hm = −79
cm. The objective function, Eq. [3.6.2–5], includes the matric head measurements
inside the soil core (j = 1) and cumulative drainage volume (j = 2) vs. time. It is rec-

976

CHAPTER 3

ommended to select a relatively small time interval between measurements (input
to data logger), and to smooth and eliminate data later, rather than selecting the number of required measurements (nj) a priori. Additional comments that provide guidance towards successful application of the multistep outflow method are:
1. Although the limitations with respect to the experiment or flow modelingare few, the inverse approach relies on the availability of a universally applicablenonlinear optimization algorithm. Problems with the parameter optimization technique generally are associated with the difficulty of defining an objective functionthat will yield unique and convergent solutions. Since ill-posedness of the inverse
problem can be caused by correlation between the parameters to be optimized,
uniqueness of the optimized parameters is generally increased by reducing the number of free parameters. For example, if parameter values can be measured independently, their values should be fixed or be included in the objective function as
Bayesian estimates (Section 1.7). To further reduce the number of parameters to
be optimized, we often make use of the relationship between m and n (m = 1 − 1/n)
and couple the retention curve model with the conductivity curve model (vanGenuchten, 1980). Furthermore, the tortuosity parameter of the conductivity function, l, is frequently set to a fixed value (l = 0.5 according to Mualem, 1976; l = −1
according to Schaap & Leij, 2000).
2. Because the Richards equation includes the soil water retention function,
θ(hm), only by way of its derivative, ∂θ/∂t = C(hm)∂hm/∂t, where C(hm) = dθ/dhm,
the residual and saturated water contents (θr and θs) are perfectly correlated. Consequently, only one of these two parameters may be optimized. In practice, θs, isindependently measured from oven-drying at the conclusion of the outflow experiment. The parameters θr and θs can be estimated simultaneously only if some additional water content-related information is included in the optimization (e.g., theinitial condition is given in terms of water content; the initial or final water volume
in the sample is included into the objective function).
3. It must be stressed that uniqueness of the coupled retention and conductivity solution will depend on the selection of the hydraulic models and their ability to accurately represent the true soil hydraulic properties. For example, Zurmühl
and Durner (1998) and Durner et al. (1999a) demonstrated that optimization witheight or more fitting parameters was stable and unique if bimodal hydraulic models were used for soils with a bimodal pore-size distributions. Furthermore, a nonsuitable retention function can lead to a meaningless estimate of the conductivity
parameter, l (Durner et al., 1999b).
4. Extrapolation beyond the range of measurement is associated with a highlevel of uncertainty. To increase the range of validity of the optimized hydraulic functions, it is recommended to include independently measured soil water retention
and/or conductivity data in the objective function. For example, with the experimental range limited by the largest applied pressure, soil water retention points atlower matric head values, as measured with the pressure plate extractor (Section
3.3.2.4) or the evaporation method (Section 3.6.2.4), may be included to augment
the water content range of the optimized hydraulic functions.
5. Outflow methods provide few dynamic data at near saturation. This implies that the sensitivity of the method to the shape of the conductivity function nearTHE SOIL SOLUTION PHASE

977

saturation, and in particular the sensitivity for estimating the saturated conductivity parameter, Ks, is low. Nevertheless, it is frequently recommended not to fix Ks
at the value of an independently measured saturated conductivity (Section 3.4) but
to treat it as an empirical fitting parameter because this often improves the description of the overall conductivity function, particularly for undisturbed soils.
However, disagreement between fitted and measured Ks may indicate that theshape of the assumed parametric form of the conductivity function near saturation
might not be correct (Durner, 1994).
6. As was demonstrated by Eching and Hopmans (1993a), the multistep outflow experiment can be extended to optimize soil hydraulic parameters for soil wetting by reversing the time sequence of the pressure increments starting with an initially dry soil. Hence, combined analysis of a drying experiment followed by awetting experiment allows simultaneous optimization of hysteresis (Schultze et al.,
1996; Zurmühl, 1998). However, as reported by Durner et al. (1999b), poor agreements between optimization results and independent measurements might occur if
an incorrect hysteresis model is used.
7. Recently, the choice of appropriate weighting factors in the objective function, Eq. [3.6.2–5], has been put into question, especially with regard to the calculations of the uncertainty of the optimized parameters. A poor choice of the measurement error can result in either too narrow or extremely large confidence intervals(Hollenbeck & Jensen, 1998a), thereby leading to incorrect interpretations of the
optimization results. Using a parameter sensitivity analysis, Vrugt et al. (2001b)
showed that uniqueness of a multistep outflow experiment without matric headmeasurements can be significantly improved by including in the objective function
only the outflow measurements immediately following the applied pressure increments, when flow rates are highest, in addition to the total outflow at the end of eachpressure step.
8. The ability of estimating soil hydraulic functions using inverse procedures
has raised the question whether hydraulic properties might be influenced by the flowrate or boundary conditions. The effect of the flow dynamics can be caused by water
entrapment or by discontinuity of water-filled pores for large matric head gradients(see Section 3.3.1). Research is ongoing (Wildenschild et al., 2001; Mortensen etal., 1998) to investigate these dynamic effects. Additional complications were reported by Hollenbeck and Jensen (1998b), who addressed the difficulty of experimental reproducibility, and Schultze et al. (1999), who cautioned that outflow experiments might include air-phase effects that should be incorporated by two-phaseflow modeling instead of describing water flow by the traditional Richards equation (Eq. [3.6.2–1]) only.
9. Some have expressed concerns about the step-wise changes of the boundary condition, which may cause a flow behavior that does not occur in nature (vanDam et al., 1994). As an alternative to the multistep method, a continuous outflow
method was tested using a gradual change of the pressure boundary condition.
Durner et al. (1999b) discussed the comparison of the multistep outflow methodwith this continuous method. They concluded that the two methodologies areequally suitable for identifying the water retention parameters. However, the small
flux rates reduced the sensitivity of the hydraulic conductivity, thereby making thecontinuous outflow method less suitable.

978

CHAPTER 3

10. Finally, measured drainage rate and matric head data can also be used toestimate the unsaturated hydraulic conductivity data directly. The concepts of thedirect hydraulic conductivity estimation were discussed by Eching et al. (1994) and
Liu et al. (1998). The procedure requires estimation of the matric head at thesoil–membrane interface. Since outflow rates are relatively high at the beginningof each pressure step, it is preferred to use data from the time periods immediately
following an increase in applied air pressure. If so required, the matric head at thesoil–membrane interface is estimated from the saturated hydraulic conductivity andthickness of the saturated porous membrane in combination with known measured
values of the water pressure at the bottom of the membrane and the measured
drainage rate. The effective permeability of the soil can be estimated subsequently
from the Darcy equation solving for K(Se), after substituting the average drainage
rate and the assumed hm gradient in the soil core using the measured matric headvalues in the center of the core and at the soil–membrane interface (Wildenschild
et al., 2001).
3.6.2.4 Evaporation Method
3.6.2.4.a Introduction
The parameter estimation technique has also been successfully applied to thelaboratory evaporation technique (Section 3.6.1.1.c). The evaporation method was
first introduced by Gardner and Miklich (1962), who imposed a series of constantevaporation rates to one side of a soil sample after first allowing the sample to attain hydraulic equilibrium before the next evaporation rate was applied. They
measured the matric head response of two tensiometers. Becher (1971) simplified
the evaporation method by using a constant evaporation rate. Several other modifications of the evaporation method, with simultaneous measurements of evaporation rate and matric head values at different heights in the sample, have since beendeveloped (Wind, 1968; Boels et al., 1978; Schindler, 1980; Tamari et al., 1993;
Wendroth et al., 1993; Halbertsma & Veerman, 1994). Experimental data obtainedwith the evaporation method can be analyzed using either a simple approach introduced by Schindler (1980), the classical Wind analysis or its modifications(e.g., Wind, 1968; Wendroth et al., 1993; Halbertsma & Veerman, 1994), or increasingly more often from numerical inversion (Ciollaro & Romano, 1995; Santini et al., 1995; Šimçnek et al., 1998c,1999c).
In their review, Feddes et al. (1988) obtained reasonable agreement between
hydraulic conductivities determined from parameter estimation by the inverse procedure and using Wind’s method (Wind, 1968). As part of a study of the spatial variability of soil hydraulic properties, Ciollaro and Romano (1995) successfully used
the inverse procedure to estimate parameter values, based on evaporation experiments, to determine the soil hydraulic parameters of a large number of samples. Santini et al. (1995) also used parameter estimation by the inverse procedure in connection with evaporation experiments. Their results compared favorably withindependently measured retention and saturated hydraulic conductivity data.
Šimçnek et al. (1998c) obtained excellent correspondence between retention curves
and hydraulic conductivity functions obtained with parameter estimation by the inverse technique and the modified Wind’s method (Wendroth et al., 1993). They also

THE SOIL SOLUTION PHASE

979

showed that, contrary to Wind’s method, which requires matric head measurements
at several locations, comparable results could be obtained with tensiometer readings at just one location. Romano and Santini (1999) showed that the inversion results compared satisfactorily with soil hydraulic data as measured independently
with the instantaneous profile method (Section 3.6.1.2.a). However, they concluded that in some cases unsaturated hydraulic conductivity functions with parameters independent of the soil water retention functions are needed for accurate
hydraulic characterization.

3.6.2.4.b Experimental Procedures
Detailed experimental procedures are outlined in Šimçnek et al. (1998c) andRomano and Santini (1999) and will only be summarized here. Initially, saturated
10-cm-high soil cores are placed on a ceramic plate to establish a static initial matric head distribution. Two to five miniature tensiometers are placed horizontallyat different vertical locations in the soil core. After tensiometer equilibration, as evidenced from the tensiometer readings, the soil core is placed on an impermeable
plate for the evaporation experiment. Total soil water storage changes are determinedfrom weight measurements of the soil core. A strain-gauge load cell placed under
the plate bearing the soil sample is used for soil sample weight measurements, while
the matric head is monitored at the various soil depths connecting the tensiometers
to pressure transducers. Calibrating each transducer across the working pressure
range should be performed before and after the test. An evaporation rate can be induced using either natural laboratory conditions, or can be accelerated using a fan
to blow air across the sample surface. A two-rate experiment can also be used inorder to induce initially sufficiently large matric head gradients (Wendroth et al.,1993). Once the hm gradient is about 1.5 to 2.5 m m−1, evaporation is allowed tocontinue without the fan. The experiment is terminated after the matric head values become too low for reliable functioning of the tensiometers.
Additional comments regarding the experimental setup are:
1. Extrapolation beyond the measurement range is associated with a high levelof uncertainty. Inclusion of independently measured information beyond the measurement range, that is, additional soil water retention data or a residual water content value, could greatly decrease this uncertainty.
2. Using the inverse modeling approach, it has been shown that matric headreadings from a single tensiometer in combination with a final soil water storage
measurement may be adequate to guarantee precise estimation of the soil hydrauliccharacteristics within the range of measurements (Šimçnek et al., 1999c). Maximum sensitivities are generally attained by placing this single tensiometer near theevaporating soil surface where the largest soil water content changes occur. However, tensiometer placement very close to the soil surface is discouraged, as the soil
surface remains dry for most of the latter part of the evaporation experiment, so thatits information is limited. Moreover, to guarantee sufficient time measurements ofchanges in matric head, especially when hydraulic head gradients become large or
as backup information if a tensiometer fails, at least two tensiometer locations arerecommended.

980

CHAPTER 3

3.6.2.4.c Simulation and Optimization
The governing flow equation for one-dimensional isothermal Darcian flow
is given by Eq. [3.6.2–1]. Boundary and initial conditions applicable to an evaporation experiment are as follows:hm(z,t) = hm,i(z)t = 0, 0 < z ≤ L

q(z,t) = qevap(t)t > 0, z = L

q(z,t) = 0

t > 0, z = 0[3.6.2–7]

where qevap(t) is the time-variable evaporation rate (L T−1) imposed at the soil surface, and all other variables are as defined in Eq. [3.6.2–6].
The objective function to be minimized includes all matric head measurements
and a soil water storage value. This single water storage value provides a water content reference from which the saturated water content can be estimated using theevaporation rate data. The inclusion of more values of water storage in the objective function does not improve the optimization process, since the evaporation rate
(upper boundary condition) is always enforced. A sensitivity analysis can be used
to determine the optimum tensiometer location for the parameter estimation procedure.
Additional comments include:
1. Parameter sensitivity is not affected by soil core height. However, duration of the evaporation experiment will be less for shorter soil cores. Šimçnek etal. (1999c) analyzed the sensitivity of the parameters of the van Genuchten relationship for a hypothetical two-rate evaporation experiment with tensiometer locations at depths of 1, 3, 5, 7, and 9 cm within a 10-cm-high sample. Their results
showed that (i) the sensitivity increased as the experiment progressed in time andthe soil core became drier, (ii) the sensitivity was highest for tensiometer 1, closest to the evaporating soil surface, and (iii) the matric head was most sensitive to nand θs, whereas the sensitivities of the parameters α, θr, and Ks were, by comparison, much smaller.
2. While a two-rate evaporation experiment has important advantages over a
one-rate experiment when using the modified Wind method, the two-stage approach
did not show an advantage in the parameter estimation approach by the inverse
method, except that it will speed up the experiment. In addition, it was concluded
by Romano and Santini (1999), who used a parameter uncertainty analysis of hypothetical numerical experiments for coarse-textured soils, that smaller parameter
confidence intervals were obtained when evaporative fluxes increased.
3. When numerically solving the evaporation process, certain combinationsof soil hydraulic parameters, selected by the parameter optimization procedure, maycause physically impossible water delivery rates towards the soil surface to satisfy
the required evaporation flux. This may especially occur at the later stages of theevaporation experiment, when water contents near the soil surface are low. To eliminate convergence problems in such situations, the evaporation rate boundary condition must be replaced by a minimum allowable matric head value (about −100 toTHE SOIL SOLUTION PHASE

981−150 m) at the expense of disagreement between simulated and measured evaporation fluxes. Since this will usually result in larger residuals, selected hydraulic parameter values must be adjusted in the subsequent iterations so that imposed evaporation fluxes can be enforced.
4. In principle, the evaporation and the multistep outflow methods represent
similar flow processes. In practice, however, there are some important differences.
First, the multistep outflow method imposes abrupt step-wise changes of the pressure boundary condition, thereby inducing high flow rates, particularly under wet
conditions. In contrast, the evaporation method imposes a smooth change of theboundary condition, which is more typical of a natural drying process. Second, as
the soil drains, the outflow rates in the multistep method decrease due to the hydraulic head gradient approaching zero, thereby reducing the sensitivity of the hydraulic parameters. In the evaporation method, however, the hydraulic head gradients increase as the soil dries, thereby gaining parameter sensitivity within theoperable range of tensiometric measurements (usually hm > −800 cm). Consequently,
the evaporation method is preferable for hydraulic characterization in the intermediate water content range, whereas the multistep outflow method is recommended
for estimation of soil hydraulic functions in the wet range. As stated above, the rangeof validity can be extended for both methods by including independently measured
data in the objective function.
3.6.2.5 Tension Disc Infiltrometer3.6.2.5.a Introduction
Tension disc infiltrometers have recently become very popular devices for insitu measurements of the near-saturated soil hydraulic properties (Perroux & White,
1988; Ankeny et al., 1991; Reynolds & Elrick, 1991; Logsdon et al., 1993). Tension infiltration data have been used primarily for evaluating saturated and unsaturated hydraulic conductivities, and to quantify the effects of macropores and preferential flow paths on infiltration. Tension infiltration data are generally used toevaluate the saturated hydraulic conductivity, Ks, and the sorptivity parameter in
Gardner’s (1958) exponential model of the unsaturated hydraulic conductivity
using Wooding’s (1968) analytical solution (see Section 3.5.4). Adequate parameter estimation requires either infiltration measurements using two different disc
diameters (Smettem & Clothier, 1989), or infiltration measurements using a single
disc diameter but multiple tensions (e.g., Ankeny et al., 1991).
Although early-time infiltration data from the tension disc infiltrometer canbe used to estimate the sorptivity (White & Sully, 1987) and the matrix flux potential,
only steady-state infiltration rates are usually used for Wooding-type analyses.
Šimçnek and van Genuchten (1996, 1997) suggested using the entire cumulative
infiltration curve in combination with parameter estimation to estimate additionalsoil hydraulic parameters. From an analysis of numerically generated data for one
supply tension experiment they concluded that the cumulative infiltration curve byitself does not contain enough information to provide a unique inverse solution.
Hence, additional information about the flow process, such as the water content
and/or matric head measured at one or more locations in the soil profile is needed
to successfully obtain unique inverse solutions for the soil hydraulic functions.

982

CHAPTER 3

3.6.2.5.b Experimental Procedures
Details of the experimental procedure are outlined in Section 3.5.4. Usually,
the soil surface is covered with a thin layer of sand to ensure hydraulic contact between the disc and the underlying soil. The sand should have a sufficiently large
air-entry value to remain saturated at all applied tensions. It should also have a Ks
value that is large enough to prevent extensive flow impedance effects (Reynolds& Zebchuk, 1996). The disc is connected with a tension-controlled water supply
system, and transient infiltration rates can be accurately determined from pressuretransducers connected to a data logging system (Ankeny et al., 1988). Water is supplied initially at the greatest suction and is decreased consecutively to a lower suction each time steady state has been attained. If desired, time domain reflectometry (TDR) or tensiometers can be installed below the supply disc to monitor watercontent and/or matric head changes (Wang et al, 1998; Šimçnek et al., 1999d). Alternatively, soil samples can be taken before and after an experiment (near and below
the supply disc) and be used for initial and final volumetric water content determinations.
Šimçnek and van Genuchten (1997) studied infiltration at several consecutive supply tensions. They considered several scenarios with different levels of information and concluded that the most practical experiment for inverse estimationof the hydraulic parameters is the cumulative infiltration curve measured at several
consecutive tensions, augmented with the initial and final water content values directly below the disc.
3.6.2.5.c Simulation and Optimization
Simulation of flow from the tension disc infiltrometer requires use of the radially symmetric two-dimensional flow equation (Eq. [3.6.2–2]) The following initial and boundary conditions applyθ(r,z,t) = θit = 0, 0 < r < ∞, 0 < z < ∞hm(r,z,t) = hm,o(t)t > 0, 0 < r < ro, z = 0

q(r,z,t) = 0

t > 0, r > ro, z = 0hm(r,z,t) = hm,it > 0, r → ∞, z → ∞[3.6.2–8]

where θi and hm,i denote constant initial water content and matric head (at zero timeand far away from the infiltration source), respectively, ro is the infiltration disc radius (L), and hm,o is the time-variable water supply pressure at the infiltrometer membrane (L). Flow symmetry is assumed around the r = 0 axis, which is therefore also
a no-radial-flux boundary.
The objective function defined in Eq. [3.6.2–5] includes cumulative infiltration data for the successively supplied water pressures, and may include additional
information such as matric head and/or water content data below the supply discat various (r,z) positions as a function of time. In addition, the objective function
could include a Bayesian estimate (Section 1.7) of the unsaturated hydraulic con-THE SOIL SOLUTION PHASE

983

ductivity computed from Wooding’s (1968) analysis. We provide the following additional comments:
1. The inverse method was tested by Simçnek et al. (1998a) using data collected as part of the soil hydrology program of the HAPEX-Sahel regional-scale
experiment (Cuenca et al., 1997). A tension disc diameter of 25 cm with supply tensions of 11.5, 9, 6, 3, 1, and 0.1 cm was used. Agreement between the measured
and optimized cumulative infiltration curves was very good. The inverse method
and Wooding’s (1968) analysis gave almost identical unsaturated hydraulic conductivities for matric head values in the interval between −2 and −10.25 cm. However, the hydraulic conductivity in the highest-matric head range (between −1 and−0.1cm) was overestimated by a factor of two using Wooding’s analysis.
2. The Šimçnek et al. (1998a) results also suggest that tension disc experiments may provide adequate information to estimate both the unsaturated hydraulic conductivity and the soil water retention properties, and that there is no needfor additional tensiometer and TDR measurements to better define the inverse
problem. Parameters of the soil water retention curve can be closely coupled withthose of the unsaturated hydraulic conductivity function, as is the case with the van
Genuchten–Mualem model. The fitting of the cumulative infiltration curves was improved by allowing the pore-connectivity parameter to be optimized as well, rather
than assuming a constant value of l = 0.5.
3. The public domain program DISC (Šimçnek & van Genuchten, 2000) foranalyzing tension disc infiltrometer data to estimate parameter values by the inverse
procedure is available from the U.S. Salinity Laboratory, Riverside, CA(http://www.ussl.ars.usda.gov) and from Soil Measurement Systems, Tucson, AZ.
3.6.2.6 Field Drainage
3.6.2.6.a Introduction
Certainly the most popular experiment to determine the soil water retention
and hydraulic conductivity functions in the field (in situ) has been the instantaneous
profile method (Hillel et al., 1972; Vachaud et al., 1978; Section 3.6.1.2.a). The first
application of a parameter optimization technique to instantaneous profile field data
was carried out by Dane and Hruska (1983). These authors questioned the uniqueness of the solution and concluded that the sensitivity of the optimized parameters
depended on the prescribed boundary conditions. Kool et al. (1987) performed asuccessful inversion using lysimeter data and a slightly different approach than Daneand Hruska (1983). Kool and Parker (1988) investigated the numerical inversionof hypothetical in situ infiltration and redistribution flow events and showed the advantages of including simultaneous measurements of pressure heads and water contents in the optimization procedure.
Santini and Romano (1992) provided guidelines for an optimal design of the
field drainage experiment and explored the feasibility of simplifying the experimental procedures while maintaining acceptable parameter uncertainty. For example, they proposed reducing the number of θ and hm measurements within a vertical soil profile. To overcome limitations regarding the lower boundary condition,
Romano (1993) investigated the effectiveness of defining a fixed hydraulic head gra-

984

CHAPTER 3

dient and the possibility of including it as an additional optimization parameter value.
However, this approach may not apply to layered soils. Zijlstra and Dane (1996)
applied the parameter optimization technique to layered soils, and concluded thatthe inverse problem may become ill-posed because of the increased number of optimized parameters. Their soil profiles consisted of one, two, or three distinct horizons, and they used water content measurements as a function of time and depth,
and matric head values as a function of time at the bottom boundary, to minimize
the objective function. The ill-posedness of the inverse problem depended on thesize of the input data set.
3.6.2.6.b Experimental Procedures
Parameter estimation of unsaturated soil hydraulic properties from an in situ
transient drainage experiment involves soil water content and matric head measurements at multiple times and soil depths. The field setup is basically the same as
for the classic instantaneous profile method (Section 3.6.1.2.a), but some simplifications are possible and will be discussed.
First, a profile description is needed to determine the major horizon boundaries. The experiment should be conducted on a leveled, vegetation-free plot withan area not smaller than 12 m2. The plot should be instrumented in the center withmeasurement sensors located at depths depending on the soil profile description.
The experiment is initiated by ponding water over the entire plot area until matrichead values are approximately zero or do not change with time at any of the depth
locations. After the water supply is stopped and the ponded water has infiltrated,
the plot is covered with a plastic sheet and insulated to ensure a zero-flux top boundary condition and to reduce temperature variations. Subsequently, the soil is allowed
to drain by gravity. The transient drainage process is monitored by simultaneous
measurements of soil water content and matric head until no changes with timeoccur. Tensiometers are usually installed at 15-cm depth increments with at least
one tensiometer in each soil horizon. Volumetric water contents can be conveniently
monitored with a calibrated neutron probe or a TDR system for the same depth increments.
As expected, the type of variables measured and their measurement locations,timing, and frequency can have a significant influence on the well-posedness andaccuracy of the parameter estimation method. If the soil domain can be consideredhomogeneous, the unsaturated hydraulic parameters can be reasonably identified
without much loss in parameter certainty by monitoring the change with time of
matric head at a single depth and water content at two other depths close to the soil
surface. It is most convenient to impose a constant total hydraulic head gradient (e.g.,
unit-gradient) at the lower boundary. The reduction in sampling depth, however,should not be at the expense of the measurement frequency. Data must be collected
for the entire drainage experiment, especially since parameter sensitivities haveshown to increase with drainage time. Generally, an accurate, simultaneous estimation of the hydraulic parameters requires between 6 and 10 measurement times
during the drainage experiment.
For a hypothetical single-layer soil profile, 90 cm in depth, Santini and Romano (1992) discussed the parameter estimation procedure, using perturbed drainage

THE SOIL SOLUTION PHASE

985

data to simulate random measurement errors. Successful results were obtained frommeasurements of a single tensiometer at the 30-cm depth in combination with watercontent values at depths of 30 and 60 cm. Other issues are discussed in the following
comments:
1. The definition of the selected lower boundary condition of the draining soil
requires careful consideration. To correctly select the appropriate boundary condition, accurate and numerous measurements of hm at the lower boundary areneeded during the course of the drainage experiment. Taking these measurementsfor proper definition of the lower boundary condition undoubtedly increases field
operations and makes the experiment more time-consuming, but it is required forthe accurate solution of Eq. [3.6.2–1]. In this regard, the assumption of a time-invariant hydraulic gradient at the bottom of the flow domain significantly reduces
the number of measurements needed, and can be especially worthwhile if many fieldmeasurements are needed, such as in spatial variability studies (Romano, 1993).
However, this assumption will require field verification. Alternatively, if appropriate, the ∂H/∂z value at the lower boundary can be considered as an additional unknown constant parameter to be estimated by the inverse procedure together withthe unsaturated soil hydraulic parameters.
2. Soil water content measurements should be taken frequently near the soil
surface where their sensitivity to the hydraulic parameters is the greatest. The neutron probe is not recommended for near surface measurements, since its large
measurement volume may extend to above the soil surface, especially under dry soil
conditions. In contrast, TDR is especially suited for near surface measurements(Topp & Davis, 1985). Time domain reflectometry has the added advantage that it
can provide automatic and real-time monitoring of soil water content (Baker & Allmaras, 1990). The depth location of the soil matric head sensors is less critical, as
it has been shown that the sensitivity of the matric head to the hydraulic parameters is largely independent of soil depth in field drainage experiments.
3. The drainage method requires a large plot (≥12 m2) to ensure that lateralwater movement at the plot boundary will not influence the water regime in the plot
center. In that regard, soil horizon characterization is important, since the presenceof impeding horizons may induce lateral water movement and affect the required
plot size. Moreover, an impeding layer can prevent unsaturated conditions below
the impeding soil layer, thereby limiting the valid water content range of the optimized hydraulic parameters. Also, multilayered soil profiles can create perchedwater tables. Lateral water movement can be reduced by constructing an impermeable boundary around the plot perimeter to a depth slightly greater than the lower
boundary of the soil domain.
3.6.2.6.c Simulation and Optimization
Equation [3.6.2–1] is solved for the following initial and boundary conditions
hm(z) = hm,it = 0, 0 < z ≤ L

q(z,t) = 0

t > 0, z = L[3.6.2–9]

986

CHAPTER 3

where z = L denotes the soil surface, with L equal to the depth of the measured soil
domain. Several types of boundary conditions can be considered at the bottom of
the soil profile (z = 0). Soil water content or pressure head values may be prescribed
as a function of time, or a constant total hydraulic head gradient, ∂H/∂z, can be used.
Representing the lower boundary condition by a time-invariant unit hydraulic gradient can accurately describe observed drainage field studies (Libardi et al., 1980;
Ahuja et al., 1988; McCord, 1991). Unknown hydraulic parameters are estimatedby minimizing Eq. [3.6.2–5] using all hm and θ measurements. Additional comments
include:1. A limitation of the drainage experiment approach is that the experimenttakes a long time (Baker et al., 1974) and that the experimental water content range
is rather small. The presence of a shallow water table or slow drainage largely reduces the water content range. Even though the method can be extended to include
surface evaporation, success of the inverse method is quite sensitive to the upperboundary condition and would therefore require accurate evaporation rate measurements.
2. Parameter optimization using drainage data can be especially useful for hydraulic characterization of field soils with little horizon differentiation, so that the
average hydraulic behavior of the entire soil profile can be characterized. The inverse technique can be further exploited to provide effective soil hydraulic parameters for simulation of hydrological responses of large-scale areas (Kabat et al.,1997).
3.6.2.7 Additional Applications
Many more experimental techniques have demonstrated the potential benefits of inverse modeling to estimate soil hydraulic functions. Briefly, we summarize the upward infiltration method (Hudson et al., 1996), the sorptivity method (Section 3.5.3), the cone penetrometer method (Gribb, 1996), and the multistep extraction
method (Inoue et al., 1998).
Both outflow and evaporation experiments represent water extraction
processes and provide parameter estimates for the draining branches of the soil watercontent and hydraulic conductivity relationships. Parameters for the wetting
branches of these soil hydraulic relationships are typically obtained from infiltration processes. The upward infiltration experiment suggested by Hudson et al. (1996)
and the sorptivity method represent typical infiltration applications in the laboratory. For the upward infiltration experiment, packed cores are placed in a flow cell,similar to a Tempe cell (Section 3.3.2). However, the imposed flux boundary condition does not require the presence of a porous membrane at the bottom end of theflow cell. Instead, nylon fabric between the bottom end cap and the soil core allows
for unrestricted water flux into the soil core. The top end plate of the cell includes
a porous stainless-steel plate to prevent soil swelling, minimize surface evaporation, and allow escape of displaced air through the top of the cell. Moreover, tensiometers are placed horizontally at different vertical positions in the soil core. A
TDR probe can be installed horizontally in the center of the soil core at the samevertical position as one of the tensiometers. However, Hudson et al. (1996) showedthat little is gained by including water content measurements in the objective func-

THE SOIL SOLUTION PHASE

987

tion. During the upward infiltration experiment, a constant water flux is imposedat the bottom of the soil sample and the matric head distribution within the soil sample is measured with one or more tensiometers. The latter should be installed only
when the visible wetting front approaches each measurement location to prevent
their malfunctioning at low initial matric head values. Upward infiltration rates canbe controlled by a syringe pump to facilitate slow and constant movement of the
one-dimensional wetting front through the soil core.
In the sorptivity method of Bruce and Klute (Klute & Dirksen, 1986), wateris allowed to infiltrate under tension into a horizontal soil column. Using a Boltzmann transformation, the one-dimensional water content form of the flow equation
can be solved analytically for this particular boundary condition, so that the soil
water diffusivity can be computed as a function of water content. The analytical solution requires measurement of the soil water content distribution along the soil column for a specific time and supply head. In Šimçnek et al. (2000), the experimental results of Nielsen et al. (1962) are compared with the inverse solution for watersupply heads of −2, −50, and −100 cm. Excellent agreement was found between
the inverse solution and the analytical solution.
A modified cone penetrometer is instrumented with a porous filter close to
the penetrometer tip and two tensiometer rings above the filter. The device ispushed into the soil to the desired depth, and a constant positive water head is applied to the filter for a short time period. While the volume of water imbibed by thesoil is monitored, so is the response of both tensiometers to the advancing wetting
front. Subsequently, after the water supply is cut off, the tensiometer measurements
record the soil drying during water redistribution. The inverse method was used toestimate the hydraulic parameters of both the wetting and drying branches of thesoil hydraulic characteristics by simultaneously analyzing the infiltration and redistribution data (Kodešová et al., 1999; Šimçnek et al., 1999a).
A multistep extraction device consists of a ceramic soil solution sampler inserted into a wetted soil and subjected to a series of vacuum extraction pressures.
The cumulative amount of soil solution extracted, as well as the response of various tensiometers near the extractor are included in the objective function to estimate the soil hydraulic functions. The method was tested both in the laboratory andin the field (Inoue et al., 1998). The cone penetrometer and the extraction method
require solution of a two-dimensional radial flow equation (Eq. [3.6.2–2]), and are
typically applied in the field, potentially to large depths.
3.6.2.8 Example
We demonstrate the application of inverse modeling by comparing one-stepand multistep outflow experiments for estimation of the hydraulic properties of a
Columbia fine sandy loam (Wildenschild et al., 2001). We will illustrate the procedure in 12 consecutive steps. Practical comments are added where applicable toassist in solving the inverse problem. Although this example is specific to the outflow method, the same steps are required for all inverse methods.
Step 1: Problem Definition. This example determines the drainage curve ofa Columbia fine sandy loam using inverse modeling. It is assumed that the underlying flow process is described by the Richards equation, (Eq. [3.6.2–1]), and that

988

CHAPTER 3

the soil hydraulic properties are adequately represented by the coupled van
Genuchten–Mualem model (Eq. [3.6.2–3] and [3.6.2–4]). Note that although these
are common assumptions, they may not be appropriate in any specific case and may
even cause unacceptable model errors.
The analysis was conducted for a disturbed soil sample that was packed to apredetermined bulk density after sieving the soil through a 0.5-mm sieve. Consequently, the hydraulic behavior is determined mostly by its texture. If the focus of
the analysis were on the estimation of hydraulic properties near saturation, which
are largely controlled by soil structure, an undisturbed soil sample should be used.
In that case, however, the measurement of additional replicates, using larger soilcores, is recommended to account for the inherent field soil spatial variability. For
the purpose of this example, inverse modeling results are compared with the classical one-step and multistep experiments.
Step 2: Selection of Measurement Types. As discussed in Section 3.6.2.3,
outflow methods require a record of cumulative outflow as a function of time. Additionally, tensiometric data from the draining soil core will improve the well-posedness of the inverse solution. Although for this purpose measurements from a single tensiometer would be adequate, the presented example uses data from two
tensiometers. The second tensiometer was included to allow simultaneous estimation of the unsaturated hydraulic conductivity. Moreover, the extra data may “save”
the experiment if the other tensiometer fails.
Step 3: Experimental Setup. The basic requirements for outflow methods
are described in Section 3.6.2.3. We selected the experimental setup of Wildenschild
et al. (2001). A diagram of the flow cell (3.5 cm high and 7.62-cm i.d.) with gas
and water flow controls is shown in Fig. 3.6.2–2. All connections consisted of quickdisconnect fittings (Cole-Parmer, Delrin, 1/4” NPT, 06359-72; Cole-Parmer Instrument Co., Vernon Hills, IL; www.coleparmer.com) so that the cell could be periodically detached for weighing to determine the water content at various times
during drainage. Two tensiometers were inserted 1.1 cm (Tensiometer 1) and 2.4
cm (Tensiometer 2) from the top of the soil core, respectively. The tensiometer ports
were offset laterally to minimize flow disturbance caused by the presence of the tensiometers. The tensiometers were made from 0.72-cm o.d., 0.1-MPa (1-bar), highflow tensiometer cups (model 652X03-BIM3; Soilmoisture Corp., Santa Barbara,CA), glued to ~6-mm (1/4-inch) o.d. acrylic tubing. The tensiometers extended approximately 2 cm into the sample. Two 0.1-MPa (1-bar) transducers (model
136PC15G2 Honeywell, Minneapolis, MN) were used to monitor the matric head
during the outflow experiments. Two additional ports were added on opposite sidesof the cell to allow flushing the sample with CO2 prior to wetting, so that complete
saturation of the sample is ensured at the start of the experiment.
The outlet was connected to a burette for measurement of outflow as a function of time. The burette was mounted such that water drained at atmospheric pressure and was loosely covered with plastic to prevent evaporation. A 70-cm waterpressure (1-psi) transducer (model 136PC01G2; Honeywell) was attached to the bottom of the burette to measure the cumulative drainage volume. The upper boundary condition was controlled using regulated pressurized N2. The N2 was bubbled
through a distilled water reservoir before entering the pressure cell to minimize evap-

THE SOIL SOLUTION PHASE

989

oration losses from the soil core. Two layers of 1.2-µm, 0.1-mm-thick nylon filters(MSI Magna nylon disc filters, Osmonic Laboratory Products; www.osmolabstore.com) were used as a porous membrane at the bottom of the sample. We combined two nylon filters to minimize puncturing, thereby ensuring a bubbling pressure of at least 700 cm during the outflow experiments. With a saturated hydraulicconductivity of approximately 0.025 cm h−1 (Table 3.6.2–1), the hydraulic resistance of the thin nylon membrane was low compared with other commonly used but
thicker porous membranes with comparable bubbling pressure. Water pressure
differences across our porous membrane are thus minimized during drainage of thesoil core. All outflow experiments were conducted on the same soil core. The measured average saturated water content was 0.445 cm3 cm−3. Time intervals for cumulative outflow and matric head measurements were approximately 10 s, which
was sufficient to ensure that the rapid outflow after a gas pressure change wasrecorded with adequate temporal resolution.
Step 4: Selection of Initial and Boundary Conditions. The soil sample wasinitially fully saturated and then partially drained before the start of the outflow experiment so that the matric head at the bottom of the soil core (hbottom) was −2 cm,
corresponding with the following initial condition:hm,i(z) = −2 − z

t = 0, 0 < z < L[3.6.2–10]

where z = L = 3.5 cm is the height of the soil surface above a datum set at z = 0 =base of sample. Since the same sample was used for both the one-step and multistep outflow experiments, the soil sample was resaturated prior to each drainage
experiment, using the same procedure every time, including flushing with CO2 toensure complete resaturation. This procedure ensured approximately identical initial saturation values between drainage experiments.
Full saturation at the onset of the experiment may be desired to provide reproducible initial conditions. However, for undisturbed soils, maximum saturation
is usually <90% of the total pore space. Therefore, it is recommended that hydraulicmeasurements of undisturbed field soils be conducted without CO2 flushing, so thatthe estimated soil hydraulic functions reflect the soil’s natural behavior.
The boundary conditions of the one-step experiment are selected so that the
air-entry value of the membrane is not exceeded. Moreover, applied pressure steps
are soil specific and depend on soil texture. For our sandy loam soil, the pressure
step applied was 500 cm, providing soil hydraulic data for a wide soil moisture range.
Larger pressure steps may be required to drain a sufficient amount of water from afiner-textured soil. Accordingly, the boundary conditions for the one-step outflow
experiment wereq(L,t) = 0 cm d−1

t>0

hm = −500 cm

0 < t < te, z = 0 cm

[3.6.2–11]

where te denotes the duration of the one-step outflow experiment, q(0,te) ≈ 0. The
lower boundary condition was experimentally realized by increasing the applied N2

990

CHAPTER 3

pressure at the top of the soil core to ha = +500 cm of water pressure, as expressedby Eq. [3.6.2–6].
For the multistep experiment, four pressure steps were applied in sequence,
which covered the same pressure range as the one-step experiment. For each successive pressure step, time was allowed for the sample to equilibrate, as determinedfrom a near zero drainage rate. The corresponding boundary conditions wereq(L,t) = 0 cm d−1

t>0

hm = −125 cm

0 > t > t1, z = 0 cm

hm = −250 cm

t1 > t > t2, z = 0 cm

hm = −375 cm

t2 > t > t3, z = 0 cm

hm = −500 cm

t3 > t > t4, z = 0 cm

[3.6.2–12]

Step 5: Experimental Data. The cumulative outflow (Q) and matric head(hm) measurements are presented in Fig. 3.6.2–3a and 3.6.2–3b, for the one-step(O) and multistep (M) experiments, respectively. Although many more data points
were collected than are shown in Fig. 3.6.2–3, we recommend not exceeding 100
data points per measurement type. The selected data must be representative of thewhole set and include those measurements that describe the flow dynamics of theoutflow experiment. Specifically, measurements prior to a pressure step increase
must be included, and data should generally describe continuously increasing (Q)or decreasing hm values. Quality control of the raw data includes removal of spurious data caused by failure of the monitoring and data logging equipment, and comparison of hm values with their expected equilibrium values at the end of each pressure step. Moreover, measurement uncertainty can be inferred from the time seriesof measurements.
Step 6: Definition of the Objective Function. The definition of the objective function is one of the most important steps in the inverse procedure. For the
purpose of this example, we will compare optimizations with (+) and without (−)
tensiometric data. As discussed in Section 3.6.2–2, selected weighting factors
should be equal to the inverse of the expected measurement variances. However,in practice, this information may not be available. Moreover, for cases that include
more than one measurement type, using measurement errors as weighting factors
often resulted in local minima. It is, therefore, recommended to set all individual
weighing factors (wi,j) in Eq. [3.6.2–5] equal to one and to compute weighing factors for the measurement types, vj, that are reciprocal to their average magnitude
multiplied with the number of data, nj. For Measurement Type 1 (cumulative outflow), 2 (Tensiometer 1) and 3 (Tensiometer 2), the corresponding weights in Eq.
[3.6.2–5] are then given by:–)
v1 = 0.5/(n1Q–
v2 = 0.25/(n2h m,1)–v3 = 0.25/(n3h m,2)

– denotes the time-averaged cumulative outflow value, and h– and h– repwhere Qm,1m,2
resent the time-averaged tensiometric data of the top (Tensiometer 1) and bottom

THE SOIL SOLUTION PHASE

991

Fig. 3.6.2–3a. Measured (symbols) and optimized cumulative outflow and matric head curves for the
one-step outflow experiment. In the notation of O4+ and O4−, the O stands for one-step outflow, thedigit refers to the number of parameters to be optimized, and the + or − indicates the use or lack ofuse of tensiometric data, respectively. The 1 and 2 in the graph for the matric head data refer to the
top and bottom tensiometer, respectively.

992

CHAPTER 3

Fig. 3.6.2–3b. Measured (symbols) and optimized cumulative outflow and matric head curves for themultistep outflow experiment. In the notation of M4+ and M4−, the M stands for multistep outflow,
the digit refers to the number of parameters to be optimized, and the + or − indicates the use or lackof use of tensiometric data, respectively. The 1 and 2 in the graph for the matric head data refer to the
top and bottom tensiometer, respectively.

THE SOIL SOLUTION PHASE

993

tensiometer (Tensiometer 2), respectively. In this example, the value of the objective function was calculated with vj = 1/(njσj2) and wi,j = 1, where σj2 denotes thevariance of all fitted data of measurement type j (Clausnitzer & Hopmans, 1995).
Step 7: Selection of Fitting Parameters. As stated in Section 3.6.1–3, the“classical” strategy in obtaining the van Genuchten–Mualem parameters is (i) to
optimize α, n, θr, and Ks and to constrain the parameter m by the relation m = 1 −1/n; (ii) to determine θs from an independent measurement; and (iii) to set the tortuosity parameter l to a fixed value of 0.5. However, since it has been increasingly
demonstrated that improved fitting of the unsaturated hydraulic conductivity function is achieved with including l as a fitting parameter, we compare optimizations
using l as a fitting parameter and as a fixed parameter set to 0.5 (Mualem, 1976).
The saturated water content was fixed to its average value of 0.455 cm3 cm−3.
Step 8: Inverse Simulations. Inverse simulations were conducted with the
HYDRUS-1D code (Šimçnek et al., 1998b). Inverse simulation requires the same
information as the forward problem, that is, time and geometric information, initial conditions, and boundary conditions, plus initial parameter estimates, positionof observation points, and measurement times with corresponding data and weighing factors. Additionally, some convergence parameters need to be defined. It is,however, recommended to apply the default values as provided with HYDRUS-1D,
which were determined from experience.
The soil sample was discretized using a variable grid spacing with 50 nodes,
using finer spacings at the bottom of the soil column. The porous nylon membrane
was not considered, since its resistance to water flow was negligible. If a ceramicporous plate is used, it must be explicitly represented in the model, as its hydraulicresistance will underestimate the soil’s hydraulic conductivity near saturation. The
optimizations were repeated three times with different initial estimates for the parameters to be optimized. Conveniently, initial estimates were taken from the database that is included in HYDRUS-1D.
In total, we completed eight sets of inverse simulations (Table 3.6.2–2). The
optimization results of the one-step (O) experiment will be compared with the multistep (M) experimental results for cases with tensiometric data (O4+, O5+, M4+,
M5+) and without tensiometric data (O4−, O5−, M4−, M5−), using either four fitting parameters (l is fixed; cases denoted by 4) or five fitting parameters (l is an additional fitting parameter; cases denoted by 5).
For all multistep simulations, the inverse solution converged to similar optimization results, indicating that solutions were unique without local minima. Thiswas different for the one-step simulations without tensiometric measurements(O4− and O5−), which converged towards various local minima depending on theinitial parameter values. Although impractical, improved optimization results forthese cases were obtained using initial parameter estimates equal to their optimized
values of the one-step experiments with hm measurements.
Step 9: Comparison of Observed and Simulated Data. Residual Analysis.
After evaluation of the uniqueness of an inverse solution, the next logical step is tocompare simulated results with the corresponding observations. Figures 3.6.2–3a
and 3.6.2–3b show this comparison for the optimizations of the one-step (O4+ and994

CHAPTER 3

Table 3.6.2–2. Optimization results for one-step (O) and multistep (M) outflow experiments to estimate
either four (O4+, O4−, M4+, and M4−) or five (O5+, O5−, M5+, and M5−) parameters with (+) orwithout (−) the help of tensiometric data. Values for the objective function and the average residuals
were calculated according to Eq. [3.6.2–5] and [3.6.2–13], respectively. The values in parentheses in
the last five columns refer to standard errors.
Case

Data

φ

AR, Q/hm

R2

θr

α

n

cm−1

cm
O5+

Q, hm

0.01120.0172/6.54 0.99

O4+

Q, hm

0.03850.0135/23.7 0.99

O5−Q

0.002620.00917

0.99

O4−

Q

0.008810.01230.99

M5+Q, hm 0.00960 0.0143/5.00 0.99

M4+

Q, hm

0.0249M5−Q

0.00625

0.0126

0.99

M4−

Q

0.00835

0.01670.99

0.0273/8.56 0.99

0.191
0.00724(0.00156) (0.00021)
0.02450.0101(0.0238) (0.00115)
0.09480.00520(0.0189) (0.00942)
0.00000.00970(0.00018) (0.00394)
0.194
0.00921(0.00270) (0.00013)
0.117
0.0106(0.01074) (0.00039)
0.192
0.00920(0.00475) (0.00023)
0.185
0.00897(0.00451) (0.00026)

Ksl

cm h−1
4.05(0.170)1.57(0.0787)
2.29(2.28)
1.52(0.780)3.02(0.0110)1.86(0.0755)
3.06(0.208)3.11(0.258)0.126(0.0081)1.08(0.288)0.219(0.485)1.08(0.780)0.278(0.0224)
1.09(0.09917)
0.358(0.0488)0.534(0.580)−1.23(0.0167)
0.5(fixed)−0.186(1.05)
0.5(fixed)−0.849(0.0539)
0.5(fixed)−0.318(0.176)
0.5(fixed)O4−) and the multistep (M4+ and M4−) simulations, respectively. The respective
optimization results are presented in Table 3.6.2–2. The comparison shows that thegeneral dynamics of the measurements is matched by all simulations. This is a necessary, but not sufficient condition to evaluate the accuracy of the optimized hydraulic parameters. Upon closer inspection, we notice the presence of small systematic deviations between measured and predicted values, indicating that errors
are autocorrelated. Moreover, it is apparent that not all four models fit the measured data equally well.
Regarding the systematic deviations, outflow data for the multistep experiment show that during the first pressure step the simulated outflow reaches equilibrium much faster than the observed data (Fig. 3.6.2–3b). This is not uncommon
and indicates that the process model (Richards’ equation) is unable to accurately
describe the true flow behavior in the near-saturated water content region. This effect is most likely caused by the presence of a discontinuous air phase at matric headvalues near the soil’s air-entry value (Schultze et al., 1999; Wildenschild et al., 2001).
Improved fitting of the measured data can be achieved by either applying the firstpressure step after the soil is first slightly desaturated or by solving the multiphaseflow problem. Other deviations of the simulated multistep outflow curves can be
related to the assumption of a fixed l parameter in the unsaturated conductivity model(Durner et al., 1999b).
When comparing cumulative outflow (Q) residuals (Column 4 in Table
3.6.2–2) between the one-step and the multistep method, it appears that the onestep method gives better results (i.e., the fit is better). However, it must be intuitively
clear that as the number and complexity of different measurement types increase,the difficulty of fitting all combined data will increase as well, especially if the op-

THE SOIL SOLUTION PHASE

995

timization problem contains model and/or measurement errors. Because of this, caseO5− (one-step method without tensiometric measurements and five fitting parameters) produces the lowest φ value of all eight cases (Table 3.6.2–2, third row). Inspection of Table 3.6.2–2 also clearly shows that the fitting of the cumulative outflow data is much improved if the connectivity parameter, l, in the conductivity
expression is allowed to vary: the φ values for the five-parameter cases (l fitted) arealways much lower than for the four-parameter cases, where l is fixed at 0.5.
The average residual (Column 4 in Table 3.6.2–2) includes the average deviation between the measured and fitted Q or hm values. These values can be compared with the measurement error of each measurement type for model adequacy
testing. Conservative estimates of measurement errors for the tensiometer and burette measurements are 1 and 0.011 cm, respectively. The burette measurement error
corresponds with 0.5 mL of drainage for a soil core radius of 3.81 cm (sample areaof 45.6 cm2). By inspection of Table 3.6.2–2, it is clear that the time-averaged residual (AR) of the matric head (measurement type j = 2 and 3), defined as:
3nj

AR(hm) = j=2
Σ (1/nj) Σ [h*m,j−1(ti) − hm,j−1(ti,ββ)]2
i=1[3.6.2–13]

is larger than the conservative measurement error in all cases (for notation, see Eq.
[3.6.2–5]). Specifically, the lowest average residual is 5 cm (case M5+, Table
3.6.2–2), whereas the assumed measurement error was only 1 cm. When considering the average residual for cumulative outflow (measurement type j = 1 only;
replace hm by Q in Eq. [3.6.2–13]), the AR of almost all cases is about equal to orslightly larger than the assumed measurement error of 0.011 cm. This is specificallyso for the one-step experiments without hm measurements (cases O4− and O5−),
since these are the easiest to fit.
These considerations, in addition to the observations of the residual analysis, which showed that the deviations between measurements and observations are
not randomly distributed, allow us to assess the adequacy of the flow model. We
conclude that the combined measurement and model error for the outflow experiment using the described experimental conditions is larger than the precision of themeasurement devices. The larger errors are caused by assuming that (i) the Richards’
type water flow equation, in combination with its numerical solution and the selected hydraulic model, is the correct physical model for simulating drainage of the
one-step and multistep experiments and/or (ii) the experimental conditions are
such that they exactly conform to the assumptions of the flow model. Clearly, thedivergence between observed and measured outflow during the first two steps ofthe multistep method indicate that these assumptions are not fully met (Fig.
3.6.2–3b). For example, although analytical models have proven useful in integratingknowledge of soil hydraulic properties for numerical models and other applications,they rarely fit the measured retention data exactly, thereby introducing AR values
larger than the measurement error. For a more thorough discussion of the influenceof model errors on the interpretation of inverse modeling results see Durner et al.
(1999a, b).
Correlation of Simulated and Observed Data. The fifth column of Table
3.6.2–2 lists the R2 values that quantify the correlation between measured and sim-

996

CHAPTER 3

ulated flow variables (for an optimum parameter set β). As all correlation values
are very high, their magnitude is of little value for model adequacy testing.
Summarizing the results, we conclude that the residual errors are smallenough to accept the assumed flow and hydraulic models for soil hydrological applications. Accordingly, we proceed further.
Step 10: Analysis of Parameter Values and Hydraulic Functions. Parameter Uncertainty. Columns 6 through 10 of Table 3.6.2–2 list the optimized hydraulic parameters. The values in parentheses denote the standard error of estimation as determined from the 95% confidence intervals around the optimized value.
The standard errors of all fitted parameters are computed from (Cii)1/2, where Cii isthe diagonal element of the parameter covariance matrix C (see Section 1.7, Eq.
[1.7–23]). It was already concluded from the residual analysis that these standard
error values are not very useful statistically, since the requirement of independent
residuals is violated. However, their relative magnitude may be useful when comparing parameters and hydraulic models.
Parameter Correlation. Table 3.6.2–3 shows correlation coefficients between
the optimized parameters for all eight cases. Values larger than ± 0.90 are bolded.
Obviously, high correlation between parameters unnecessarily increases the number of optimized parameters. High correlation causes underestimation of parameter uncertainty, slows down convergence rate, and increases nonuniqueness. It is expected that the number of highly correlated parameters increases as the number offitted parameters increases. As a result, the available information in the objectivefunction is reduced. Table 3.6.2–3 confirms that the highest number of correlations
occur for cases O4− and O5−, that is, the one-step outflow experiment without matric head measurements. As discussed in Section 3.6.2–3, the information content
of these experiments is not sufficient to obtain unique parameters. The data show
large correlations between n and θr for the multistep experiments, a result that iscommon when fitting retention data to the van Genuchten relationship if information for the dry range is missing.
Hydraulic Functions. The optimized soil water retention and unsaturated hydraulic conductivity functions for all four cases using four fitting parameters arepresented in Fig. 3.6.2–4. The optimized functions are accurate only for the experimental water content range of the outflow experiment (i.e., hm > −500 cm), andcare must be exercised in their extrapolation to drier soil conditions.
Independent soil water retention data (thick solid line) were obtained by a syringe pump procedure (Wildenschild et al., 1997) in a separate experiment bywhich the soil sample was drained at a constant, low flow rate of 0.5 mL h−1, simulating quasistatic conditions at any time during the drainage experiment. The conditions of the steady-state experiment are considered optimal; therefore, the optimized retention curves are compared with the retention data from the syringe
experiment using cumulative drainage with corresponding tensiometric data. Unfortunately, it was impossible to estimate the hydraulic conductivity from the syringe pump experiments because of errors resulting from the extremely small hydraulic head gradients. Consequently, the optimized curves are compared withindependently estimated unsaturated hydraulic conductivity data reported in Wildenschild et al. (2001).

THE SOIL SOLUTION PHASE

997

Table 3.6.2–3. Correlation matrices for the van Genuchten parameters for one-step (O) and multistep(M) outflow experiments to estimate either four (O4+, O4−, M4+, and M4−) or five (O5+, O5−, M5+,and M5−) parameters with (+) or without (−) the help of tensiometric data.
Parameter

Case

θr

θr

O5+
O4+
O5−
O4−
O5+
O4+
O5−
O4−
O5+
O4+
O5−
O4−
O5+
O4+
O5−
O4−
O5+
O4+
O5−
O4−
M5+
M4+
M5−
M4−
M5+
M4+
M5−
M4−
M5+
M4+
M5−
M4−
M5+
M4+
M5−
M4−
M5+
M4+
M5−
M4−

1

α

n

Ks

l

θr

α

n

Ks

l

−0.162−0.363−0.282−0.683
0.581
0.864
0.346
−0.662−0.104−0.352−0.323−0.661
0.083-0.161-1−0.499−0.481−0.547−0.650
0.913
0.950
0.887
0.912
−0.368−0.267−0.347−0.006−0.330-−0.235
--

α

n

Ks

l

1

−0.761−0.774−0.998−0.999
0.673
0.992
0.998
0.998
0.037
-0.981--1

0.715−0.765−0.998−0.992−0.063-0.971--1

0.715-−0.966--

1

1−0.737−0.715−0.804−0.8400.255
0.621−0.0560.158
0.111-−0.172--1−0.408−0.412−0.080−0.1700.322-0.077--1

0.811-0.751--

1

When comparing the functions, the results are surprisingly close. Note, however, that the optimization results of the one-step experiment without tensiometricmeasurements (O4−) were obtained by using initial, optimized values from the optimization with hm measurements (O4+) as initial parameter values. We also findthat the tensiometric measurements in the multistep experiments gave results similar to multistep experiments without tensiometric data, especially for cases with
an optimized l. This may be caused partly by our experimental strategy to wait for
static equilibrium before increasing the applied gas pressure. Generally speaking,it shows that for the conditions of our experiment, the use of outflow data alone froma multistep experiment may be sufficient to obtain the correct hydraulic properties

998

CHAPTER 3

Fig. 3.6.2–4. Comparison of optimized and independently measured retention and unsaturated hydraulicconductivity data for the four parameter optimizations (l = 0.5). For an explanation of the notationof O4+, O4−, M4+, and M4−, see the figure captions of Fig. 3.6.2–3a,b.

by inverse simulation. The unsaturated hydraulic conductivity comparison in Fig.
3.6.2–4 also supports the general conclusion that multistep experiments are preferredover one-step experiments.
Step 11: Response Surface Analysis. The evaluation of response surfaces
is usually done a priori to design the experimental conditions for inverse modelingusing forward simulations. We like to illustrate these analysis in this final step to
support the concepts introduced in Section 3.6.2–2.
Response Surface Analysis. Response surface analysis can be used to investigate the posedness of the optimization problem. The behavior of the objective function within a multiparameter space can only be visualized by limiting the numberof variable parameters to two. Consequently, the value of the objective function, φ,
can be shown for combinations of two parameters, while keeping all other param-

THE SOIL SOLUTION PHASE

999

eters at their “true” (i.e., optimized) values. It thus follows that response surfaces
can only be calculated posteriori if the true optimum is known as determined fromforward modeling results.
The case with four adjustable parameters requires analysis of six possible parameter pairs. Figure 3.6.2–5a shows response surfaces for (n,α), (Ks,α), (θr,α),(θr,Ks), (Ks,n), and (θr,n) for the one-step experiment without hm measurements (caseO4− in Table 3.6.2–3). Along a response surface, the optimum two-parameter
combination is determined by a φ valley or minimum enclosed by contour lines. The
shape of the valley indicates the rate of convergence, degree of parameter correlation, parameter sensitivity, and presence of local minima. When comparing response
surfaces, the parameter sensitivity decreases as the φ interval increases. The ideal
response surface shows a narrow minimum area with circular shape, indicating nocorrelation between fitting parameters. Response surfaces that are parallel to oneof the axes are insensitive to that respective parameter, indicating high parameter
uncertainty and a large confidence interval. Unfavorable response surfaces consistof long narrow valleys with an approximate angle of 45°, indicating a high correlation between parameters, while L-shaped valleys indicate slow convergence.
The response surfaces calculated from real data, contrary to those obtained
from numerically generated data, are much more difficult to interpret since the residuals are not necessarily normally distributed and may exhibit some systematic bias.
The response surfaces of the ill-posed one-step experiment in Fig. 3.6.2–5a are unfavorable because of (i) multiple minima within a low φ value valley (response surface θr,n), (ii) extremely narrow valleys for parameter combinations θr,n and Ks,n,and (iii) wide interval spacings. When comparing the results of Fig. 3.6.2–5a withthose of Fig. 3.6.2–5b (multistep with hm data), we note that interval spacings are
smaller and that most minima are elliptical in shape. Still, the Ks parameter appears
to be the least sensitive, which is not surprising when considering that almost alldata in the objective function are related to unsaturated flow.
Response surfaces are usually calculated to describe the behavior of the objective function in the cross sections of the two-parameter plane for the multiparameter space. However, they cannot illustrate the shape of the objective function inother parts of the parameter space. Thus, if response surfaces do not display welldefined minima, the parameter optimization problem is certainly ill-posed (Figure
3.6.2–5a, parameter pair θr,n). There is, however, no guarantee that the problem iswell-posed if well-defined minima are demonstrated for all parameter pairs. Thisis because minimization (especially for gradient-type methods; see Section 1.7) canlead to optimized parameters for any local minimum of the objective function, faraway from the global minimum. For example, although the response surfaces forthe one-step and multistep experiments look fairly similar (Fig. 3.6.2–5a and3.6.2–5b), the response surfaces for the one-step method were obtained after using
initial parameter estimates that were close to the actual global minimum. Thus, response surface analyses can show that a specific experimental setup will result in
an ill-posed optimization problem, but it cannot guarantee a well-posed optimization problem by itself.
Summary and Conclusions. Analysis of the one-step data without considering tensiometric information suffered from large uncertainties and cross corre-

1000
CHAPTER 3

Fig. 3.6.2–5a. Response surfaces for the one-step experiment without tensiometric data (O4-).

THE SOIL SOLUTION PHASE
1001

Fig. 3.6.2–5b. Response surfaces for the multistep experiment with tensiometric data (M4+).

1002

CHAPTER 3

lation among parameters, with optimized parameters strongly dependent on theirinitial estimates. Results for the multistep experiment (both with and without additional tensiometric information) and the one-step experiment with measured matric head values were fairly similar. This was especially the case when the l parameter
was allowed to be optimized as well. However, significant smaller confidence intervals were obtained when combining matric head and cumulative outflow data inthe objective function. Thus, although matric head data are not required for inverse
modeling of a multistep outflow experiment, their inclusion may lead to smaller confidence intervals of the optimized parameters. Tensiometric data are required when
conducting a one-step outflow experiment to avoid nonuniqueness problems. Finally, the analysis has shown that the inverse modeling of multistep outflow can provide accurate and reliable estimates of both the retention and the unsaturated hydraulic conductivity curve parameters for the intermediate soil water content range.
3.6.2.9 Discussion
Soil hydraulic parameter estimation by inverse modeling is a relatively complex procedure that provides a quick method for soil hydraulic characterization,yielding parameters for both the soil water retention and unsaturated hydraulic conductivity function from a single experiment. Its successful application requires suitable experimental procedures as well as advanced numerical flow codes and optimization algorithms. Numerical codes with user-friendly interfaces are becoming
available that can be used for both inverse and direct simulations (e.g., Šimçnek etal., 1998b, 1999b). However, since the method as a whole is not fully developed
yet, both experimental and numerical modeling expertise is required for successful application of the methodology and correct interpretation of the results.
When compared with other measurement methods, the inverse modeling approach renders a suite of benefits. First and foremost, it mandates the combinationof experimentation with numerical modeling. Since the optimized hydraulic functions are mostly needed as input to numerical flow and transport models for prediction purposes, it is an added advantage that the hydraulic parameters are estimated using similar numerical models as used for predictive forward modeling. An
additional benefit of the inverse procedure is their application to transient experiments, thereby providing relatively fast results. Finally, the parameter optimization
procedure computes confidence intervals of the optimized parameters, although theirinterpretation can be misleading.
Inverse problems for parameter estimation of soil hydraulic functions can be
ill-posed because of inadequate experimental design, measurement errors, andmodel errors. Analysis of such flow problems must include a search for the optimum number of flow variables required in the objective function. For example, as
the number of optimized parameters increases, increased information content of themeasurements is required, for example, by including observations of different
types, or by using time-variable boundary conditions. Sensitivity analysis canlargely optimize the need for type, number, and spatial location of the measurements.
For example, it has been shown in a variety of applications that measured transient
flow data, as induced by multiple step changes at the domain boundary, are moresensitive to the estimated parameters than using a single step. Also, an increase inTHE SOIL SOLUTION PHASE

1003

the number of parameters to be optimized will generally lead to a reduction in model
error, but will increase the parameter uncertainty. A well-posed problem requiresa priori testing for nonuniqueness using response surface analysis and parameter
sensitivity and correlation. Insensitive parameters should be measured independently, whereas highly correlated parameters will affect uniqueness of the inverse
problem, requiring independent measurement of one of the correlated parameters.
To reduce nonuniqueness and to extend the range of application beyond the experimental range of measurements, independently measured information on the soil
hydraulic functions can be included in the objective function. This prior information can reduce parameter uncertainty, but may reduce the goodness of fit between
model and data.
As with all other laboratory and field methods to estimate soil hydraulic functions, it is assumed that the functional forms used are capable of accurately describing the soil hydraulic data. It is therefore essential to perform a model adequacy
test by comparing objective function residuals with measurement errors, to identify model errors. Obviously, if the parametric models (Section 3.3.4) are not adequate for a tested soil, the resulting fitting parameters will not be valid. Especially,
if the coupled Mualem approach is used, their improper selection may compromise
the accuracy of both hydraulic functions.
When comparing the optimized hydraulic functions with the results of othermethods, one must consider differences in model assumptions and experimental
range. Laboratory measurements, although accurate, provide hydraulic information
for a relatively small soil core, detached from its surroundings. On the other hand,field experiments will generally include the continuum of soil horizons that will influence water flow and the estimated soil hydraulic functions. Moreover, as is thecase for any method, the parameter estimates are only valid for the range of the experimental conditions, and care must be exercised in their extrapolation.
This chapter must be regarded as a work in progress, since the mathematical, analytical, and experimental procedures that constitute the inverse method asa whole are still an area of intensive research. Improvements in parameter estimation methods in combination with experimental requirements and optimization algorithms continue to appear in a steady stream of publications. Nevertheless, thegeneral inverse method has demonstrated to be an excellent new tool that allows
for soil hydraulic characterization using a wide spectrum of transient laboratory andfield experiments. To date, the application of the inverse parameter estimation
method in the vadose zone has been limited to the estimation of soil hydraulic properties. This is not surprising because hydraulic parameters are required in most flow
and transport models and their direct measurements are time-consuming. We posethat inverse modeling can be used to estimate other soil properties as well, such as
solute transport, heat flow, and gaseous transport parameters. Moreover, the methodology can be applied to better understand processes, such as differentiation between
matrix and macropore contributions to water flow for two-domain simulation models (Šimnek et al., 2001), or to infer root water uptake parameters in crop growthsimulation models (Vrugt et al., 2001a). Summarizing, we conclude that parameter estimation by inverse modeling has tremendous potential in characterizing vadose flow and transport processes, while simultaneously presenting us with an additional tool to better understand their fundamental mechanisms.

1004

CHAPTER 3

3.6.2.10 References
Ahuja, L.R., B.B. Barnes, D.K. Cassel, R.R. Bruce, and D.L. Nofziger. 1988. Effect of assumed unit
gradient during drainage on the determination of unsaturated hydraulic conductivity and infiltration parameters. Soil Sci. 145:235–243.
Ankeny, M.D., M. Ahmed, T.C. Kaspar, and R. Horton.1991. Simple field method for determining unsaturated hydraulic conductivity. Soil Sci. Soc. Am. J. 55:467–470.
Baker, F.G., P.L.M. Veneman, and J. Bouma. 1974. Limitations of the instantaneous profile method forfield measurement of unsaturated hydraulic conductivity. Soil Sci. Soc. Am. J. 38:885–888.
Baker, J.M., and R.R. Allmaras. 1990. System for automating and multiplexing soil moisture measurement by time-domain reflectometry. Soil Sci. Soc. Am. J. 54:1–6.
Bear, J. 1972. Dynamics of fluids in porous media. American Elsevier, New York, NY.
Becher, H.H. 1971. Ein Verfahren zur Messung der ungesättigten Wasserleitfähigkeit. Z. Pflanzenernaehr. Bodenkd. 128:1–12.
Boels, D., J.B.H.M. van Gils, G.J. Veerman, and K.E. Wit. 1978. Theory and system of automatic determination of soil moisture characteristics and unsaturated hydraulic conductivities. Soil Sci.
126:191–199.
Brooks, R.H., and A.T. Corey. 1966. Properties of porous media affecting fluid flow. J. Irrig. Drain. Div.
Am. Soc. Civ. Eng. 92:61–88.
Celia, M.A., E.T. Bouloutas, and R.L. Zarba. 1990. A general mass-conservative numerical solution for
the unsaturated flow equation. Water Resour. Res. 26:1483–1496.
Chen, J., J.W. Hopmans, and M.E. Grismer. 1999. Parameter estimation of two-fluid capillary pressuresaturation and permeability functions. Adv. Water Resour. 22:479–493.
Ciollaro, G., and N. Romano. 1995. Spatial variability of the soil hydraulic properties of a volcanic soil.
Geoderma 65:263–282.
Clausnitzer, V., and J.W. Hopmans. 1995. Nonlinear Parameter estimation. LM-OPT. General purpose
optimization code based on the Levenberg–Marquardt algorithm. LAWR Report 100032. University of California, Davis, CA.
Cuenca, R.H., J. Brouwer, A. Chanzy, P. Droogers, S. Galle, S.R. Gaze, M. Sicot, M., H. Stricker, R.
Angulo-Jaramillo, S.A. Boyle, J. Bromly, A.G. Chebhouni, J.D. Cooper, A.J. Dixon, J.-C. Fies,
M. Gandah, J.-C. Gaudu, L. Laguerre, M. Soet, H.J. Steward, J.-P. Vandervaere, and M. Vauclin. 1997. Soil measurements during HAPEX-Sahel intensive observation period. J. Hydrol.
188–189:224–266.
Dane, J.H., and S. Hruska. 1983. In-situ determination of soil hydraulic properties during drainage. Soil
Sci. Soc. Am. J. 47:619–624.
Dirksen, C. 1991. Unsaturated hydraulic conductivity. p. 209–269. In K.A. Smith and C.E. Mullins (ed.)
Soil analysis: Physical methods. Marcel Dekker, New York, NY.
Doering, E.J. 1965. Soil water diffusivity by the one-step method. Soil Sci. 99:322–326.
Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water
Resour. Res. 30:211–233.
Durner, W., E. Priesack, H.-J. Vogel, and T. Zurmühl. 1999a. Determination of parameters for flexible
hydraulic functions for inverse modeling. p. 817–827. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. University of California, Riverside, CA.
Durner, W., E.B. Schultze and T. Zurmühl. 1999b. State-of-the-art in inverse modeling of inflow/outflow experiments. p. 661–681. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. University of California,
Riverside, CA.
Eching, S.O., and J.W. Hopmans. 1993a. Optimization of hydraulic functions from transient outflow
and soil water pressure data. Soil Sci. Soc. Am. J. 57:1167–1175.
Eching, S.O., and J.W. Hopmans. 1993b. Inverse solution of unsaturated soil hydraulic functions fromtransient outflow and soil water pressure data. Land, Air and Water Resources Paper No.
100021. Univ. of California, Davis, CA.
Eching, S.O., J.W. Hopmans, and O. Wendroth. 1994. Unsaturated hydraulic conductivity from transient multi-step outflow and soil water pressure data. Soil Sci. Soc. Am. J. 58:687–695.
Feddes, R.A., P. Kabat, P.J.T. van Bakel, J.J. B. Bronswijk, and J. Habertsma. 1988. Modeling soil water
dynamics in the unsaturated zone—State of the art. J. Hydrol. 100:69–111.
Gardner, W.R. 1956. Calculation of capillary conductivity from pressure plate outflow data. Soil Sci.
Soc. Am. Proc. 20:317–320.
Gardner, W.R. 1958. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85:228–232.

THE SOIL SOLUTION PHASE

1005

Gardner, W.R. 1962. Note on the separation and solution of diffusion type equations. Soil Sci. Soc. Am.
Proc. 26:404.
Gardner, W.R., and F.J. Miklich. 1962. Unsaturated conductivity and diffusivity measurements by a constant flux method. Soil Sci. 93:271–274.
Gribb, M.M. 1996. Parameter estimation for determining hydraulic properties of a fine sand from transient flow measurement. Water Resour. Res. 32:1965–1974.
Gupta, S.C., D.A. Farrell, and W.E. Larson. 1974. Determining effective soil water diffusivities fromone-step outflow experiments. Soil Sci. Soc. Am. Proc. 38:710–716.
Halbertsma, J.M., and G.J. Veerman. 1994. A new calculation procedure and simple set-up for the evaporation method to determine soil hydraulic functions. Report 88. DLO Winand Staring Centre,
Wageningen, The Netherlands.
Hillel, D., V.D. Krentos, and Y. Stylianou. 1972. Procedure and test of an internal drainage method formeasuring soil hydraulic characteristics in situ. Soil Sci. 114:395–400
Hollenbeck, K.J., and K.H. Jensen. 1998a. Maximum-likelihood estimation of unsaturated hydraulicparameters. J. Hydrol. 210:192–205.
Hollenbeck, K.J., and K.H. Jensen. 1998b. Experimental evidence of randomness and nonuniquenessin unsaturated outflow experiments designed for hydraulic parameter estimation. Water Resour.
Res. 34:595–602.
Hollenbeck, K.J., J. Šimçnek, and M.Th. van Genuchten. 2000. RETCML: Incorporating maximumlikelihood estimation principles in the soil hydraulic parameter estimation code RETC. Comput. Geosci. 26:319–327.
Hopmans, J.W., and J. Šimçnek. 1999. Review of inverse estimation of soil hydraulic properties. p.
643–659. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulicproperties of unsaturated porous media. University of California, Riverside, CA.
Hopmans, J.W., T. Vogel, and P.D. Koblik. 1992. X-ray tomography of soil water distribution in onestep outflow experiments. Soil Sci. Soc. Am. J. 56:355–362.
Hudson, D.B., P.J. Wierenga, and R.G. Hills. 1996. Unsaturated hydraulic properties from upward flow
into soils cores. Soil Sci. Soc. Am. J. 60:388–396.
Inoue, M., J. Šimçnek, J.W. Hopmans, and V. Clausnitzer. 1998. In situ estimation of soil hydraulic functions using a multi-step extraction technique. Water Resour. Res. 34:1035–1050.
Kabat, P., R.W.A. Hutjes, and R.A. Feddes. 1997. The scaling characteristics of soil parameters: From
plot scale heterogeneity to subgrid parameterization. J. Hydrol. 190:363–396.
Klute, A., and C. Dirksen. 1986. Conductivities and diffusivities of unsaturated soils. p. 687–734. In A.
Klute (ed.) Methods of soil analysis. Part 1. 2nd ed. Agron. Monogr. 9. ASA and SSSA, Madison, WI.
Kodešová, R., S.E. Ordway, M.M. Gribb, and J. Šimçnek. 1999. Estimation of soil hydraulic properties with the cone permeameter: Field studies. Soil Sci. 164:527–541.
Kool, J.B., and J.C. Parker. 1988. Analysis of the inverse problem for transient unsaturated flow. Water
Resour. Res. 24:817–830.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1987. Parameter estimation for unsaturated flow andtransport models—A review. J. Hydrol. 91:255–293.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties for onestep outflow experiments by parameter estimation. I. Theory and numerical studies. Soil Sci.
Soc. Am. J. 49:1348–1354.
Kulkarni, R., A.T. Watson, J.E. Nordtvedt, and A. Sylte. 1998. Two-phase flow in porous media: Property identification and model validation. AICHE J. 44:2337–2350.
Libardi, P.L., K. Reichardt, D.R. Nielsen, and J.W. Biggar. 1980. Simple field methods for estimating
soil hydraulic conductivity. Soil Sci. Soc. Am. J. 44:3–7.
Liu, Y.P., J.W. Hopmans, M.E. Grismer, and J.Y. Chen. 1998. Direct estimation of air–oil and oil–water
capillary pressure and permeability relations from multi-step outflow experiments. J. Contam.
Hydrol. 32:223–245.
Logsdon, S.D., E.L. McCoy, R.R. Allmaras, and D.R. Linden. 1993. Macropore characterization by indirect methods. Soil Sci. 155:316–324.
McCord, J.T. 1991. Application of second-type boundaries in unsaturated flow modeling. Water Resour. Res. 27:3257–3260.
Mortensen, A.P., R.J. Glass, and K. Hollenbeck. 1998. Visualization of quasi-2D unsaturated flow during dynamic outflow experiments. Am. Geophys. Union. Eos Trans. 79:368.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media.
Water Resour. Res. 12:513–522.
Nielsen, D.R., J.W. Biggar, and J.M. Davidson. 1962. Experimental consideration of diffusion analysis in unsaturated flow problems. Soil Sci. Soc. Am. Proc. 26:107–111.

1006

CHAPTER 3

Nielsen, D.R., M.Th. van Genuchten, and J.W. Biggar. 1986. Water flow and solute transport processes
in the unsaturated zone. Water Resour. Res. 22:89S–108S.
Parker, J.C., J.B. Kool, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties fromone-step outflow experiments by parameter estimation. II. Experimental studies. Soil Sci. Soc.
Am. J. 49:1354–1359.
Passioura, J.B. 1976. Determining soil water diffusivities from one-step outflow experiments. Austr. J.
Soil Res. 15:1–8.
Perroux, K.M., and I. White. 1968. Design for disc permeameters. Soil Sci. Soc. Am. J. 52:1205–1215.
Quintard, M., and S. Whitaker. 1999. Fundamentals of transport equation formulation for two-phase flow
in homogeneous and heterogeneous porous media. p. 3–57. In M.B. Parlange and J.W. Hopmans(ed.) Vadose zone hydrology: Cutting across disciplines. Oxford University Press, Oxford, UK.
Reynolds, W.D., and D.E. Elrick. 1991. Determination of hydraulic conductivity using a tension infiltrometer. Soil Sci. Soc. Am. J. 55:633–639.
Reynolds, W.D., and W.D. Zebehuk. 1996. Use of contact material in tension infiltrometer measurements. Soil Technol. 9:141–159.
Reeve, M.J., and A.D. Carter. 1991. Water release characteristic. p. 111–160. In K.A. Smith and C.E.
Mullins (ed.) Soil analysis: Physical methods. Marcel Dekker, New York, NY.
Romano, N. 1993. Use of an inverse method and geostatistics to estimate soil hydraulic conductivity
for spatial variability analysis. Geoderma 60:169–186.
Romano, N., and A. Santini. 1999. Determining soil hydraulic functions from evaporation experiments
by a parameter estimation approach: Experimental verifications and numerical studies. Water
Resour. Res. 35:3343–3359.
Russo, D. 1988. Determining soil hydraulic properties by parameter estimation: On the selection of amodel for the hydraulic properties. Water Resour. Res. 24:453–459.
Russo, D., E. Bresler, U. Shani, and J.C. Parker. 1991. Analysis of infiltration events in relation to determining soil hydraulic properties by inverse problem methodology. Water Resour. Res.
27:1361–1373.
Santini, A., and N. Romano. 1992. A field method for determining soil hydraulic properties. Proceedings of XXIII Congress of Hydraulics and Hydraulic Constructions. Florence, Italy. 31 Aug.–4Sept. Tecnoprint S.n.c., Bologna, I:B.117–B.139 (in Italian; available also in English).
Santini, A., N. Romano, G. Ciollaro, and V. Comegna. 1995. Evaluation of a laboratory inverse methodfor determining unsaturated hydraulic properties of a soil under different tillage practices. Soil
Sci. 160:340–351.
Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with theMualem–van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.
Schindler, U. 1980. Ein Schnellverfahren zur Messung der Wasserleitfähigkeit im teilgesättigten Bodenan Stechzylinderproben. Arch. Acker-Pflanzenbau u. Bodenkd. 24:1–7.
Schultze, B., O. Ippisch, B. Huwe, and W. Durner. 1999. Dynamic nonequilibrium in unsaturated water
flow. p. 877–892. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of thehydraulic properties of unsaturated porous media. University of California, Riverside, CA.
Schultze, B., T. Zurmühl, and W. Durner 1996. Untersuchung der Hysterese hydraulischer Funktionen
von Böden mittels inverser Simulation. Mitteilungen der Deutschen Bodenkundlichen
Gesellschaft 80:319–322.
Šimçnek, J., R. Angulo-Jaramillo, M.G. Schaap, J.-P. Vandervaere, and M.Th. van Genuchten. 1998a.
Using an inverse method to estimate the hydraulic properties of crusted soils from tension disc
infiltrometer data. Geoderma 86:61–81.
Šimçnek, J., J.W. Hopmans, D.R. Nielsen, and M.Th. van Genuchten. 2000. Horizontal infiltration revisited using parameter optimization. Soil Sci.165:708–717.
Šimçnek, J., R. Kodešová, M.M. Gribb, and M.Th. van Genuchten. 1999a. Estimating hysteresis in thesoil water retention function from modified cone penetrometer test. Water Resour. Res.
35:1329–1345.
Šimçnek, J., M. Šejna, and M.Th. van Genuchten. 1998b. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated
media. Version 2.0. IGWMC-TPS-70. International Ground Water Modeling Center, Colorado
School of Mines, Golden, CO.
Šimçnek, J., M. Šejna, and M.Th. van Genuchten. 1999b. The HYDRUS-2D software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably-saturated
media. Version 2.0. IGWMC-TPS-53. International Ground Water Modeling Center, Colorado
School of Mines, Golden, CO.
Šimçnek, J., and M.Th. van Genuchten. 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour. Res. 32:2683–2696.

THE SOIL SOLUTION PHASE

1007

Šimçnek, J., and M.Th. van Genuchten. 1997. Estimating unsaturated soil hydraulic properties from multiple tension disc infiltrometer data. Soil Sci. 162:383–398.
Šimçnek, J., and M.Th. van Genuchten. 2000. The DISC computer software for analyzing tension disc
infiltrometer data by parameter estimation. Versions 1.0. Research Rep. 145. U.S. Salinity Laboratory, USDA-ARS, Riverside, CA.
Šimçnek, J., O. Wendroth, and M.Th. van Genuchten. 1998c. A parameter estimation analysis of theevaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.
Šimçnek, J., O. Wendroth, and M.Th. van Genuchten. 1999c. Estimating soil hydraulic properties fromlaboratory evaporation experiments by parameter estimation. p. 713–724. In M.Th. van Genuchten
et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous
media. University of California, Riverside, CA,Šimçnek, J., O. Wendroth, and M.Th. van Genuchten. 1999d. A parameter estimation analysis of thelaboratory tension disc experiment for determining soil hydraulic properties. Water Resour. Res.
35:2965–2979.
Šimçnek, J., O. Wendroth, N. Wypler, and M.Th. van Genuchten. 2001. Nonequilibrium water flow characterized from an upward infiltration experiment. Eur. J. Soil Sci. 52:13–24.
Smettem, K.R.J., and B.E. Clothier. 1989. Measuring unsaturated sorptivity and hydraulic conductivity using multiple disk permeameters. J. Soil Sci. 40:563–568.
Tamari, S., L. Bruckler, J. Halbertsma, and J. Chadoeuf. 1993. A simple method for determining soil
hydraulic properties in the laboratory. Soil Sci. Soc. Am. J. 57:642–651.
Toorman, A.F., P.J. Wierenga, and R.G. Hills. 1992. Parameter estimation of soil hydraulic properties
from one-step outflow data. Water Resour. Res. 28:3021–3028.
Topp, G.C., and J.L. Davis. 1985. Measurement of soil water content using time domain reflectometry(TDR). Soil Sci. Soc. Am. J. 46:19–24.
Vachaud, G., C. Dancette, S. Sonko, and J.L. Thony. 1978. Méthodes de caractérisation hydrodynamique in situ d’un sol non saturé. Application à deux types de sol du Sénégal en veu de ladétermination des termes du bilan hydrique. Ann. Agron. 29:1–36.
Valiantzas, J.D., and D.G. Kerkides. 1990. A simple iterative method vor the simultaneous determination of soil hydaulic properties from one-step outflow data. Water Resour. Res. 26:143–152.
van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1992. Evaluation of the inverse method for determining soil hydraulic functions from one-step outflow experiments. Soil Sci. Soc. Am. J.
56:1042–1050.
van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1994. Inverse method for determining soil hydraulicfunctions from multi-step outflow experiments. Soil Sci. Soc. Am. J. 58:647–652.
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.
van Genuchten, M.Th., and E.A. Sudicky. 1999. Recent advances in vadose zone flow and transport modeling. p. 155–193. In M.B. Parlange and J.W. Hopmans (ed.) Vadose zone hydrology: Cutting
across disciplines. Oxford University Press, Oxford, UK.
Vogel, T., M. Nahaei, and M. Císlerová. 1999. Description of soil hydraulic properties near saturationfrom the point of view of inverse modeling. p. 693–703. In M.Th. van Genuchten et al. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. University of California, Riverside, CA.
Vrugt, J., J.W. Hopmans, and J. Šimçnek. 2001a. Calibration of a two-dimensional root water uptake
model. Soil Sci. Soc. Am. J. 65:1027–1037.
Vrugt, J., A. Weerts, and W. Bouten. 2001b. Information content of data for identifying soil hydraulicproperties from outflow experiments. Soil Sci. Soc. Am. J. 65:19–27.
Wang, D., S.R. Yates, and F. F. Ernst. 1998. Determining soil hydraulic properties using tension infiltrometers, time domain reflectometry, and tensiometers. Soil Sci. Soc. Am. J. 62:318–325.
Wendroth, O., W. Ehlers, J.W. Hopmans, H. Kage, J. Halbertsma, and J.H.M. Wösten. 1993. Reevaluation of the evaporation method for determining hydraulic functions in unsaturated soils. Soil
Sci. Soc. Am. J. 57:1436–1443.
Whisler, F.D., and K.K. Watson. 1968. One-dimensional gravity drainage of uniform columns of porous
materials. J. Hydrol. 6:277–296.
White, I., and M.J. Sully. 1987. Macroscopic and microscopic capillary length and time scales from field
infiltration, Water Resour. Res. 23:1514–1522.
Wildenschild, D., J.W. Hopmans, and J. Šimçnek. 2001. Flow rate dependence of soil hydraulic characteristics. Soil Sci. Soc. Am. J. 65:35–48.
Wildenschild, D., K.H. Jensen, K.J. Hollenbeck, T.H. Illangasekare, D. Znidarcic, T. Sonnenborg, andM.B. Butts. 1997. A two-stage procedure for determining unsaturated hydraulic characteristics
using a syringe pump and outflow observations. Soil Sci. Soc. Am. J. 61:347–359.

1008

CHAPTER 3

Wind, G. P. 1968. Capillary conductivity data estimated by a simple method. p. 181–191. In P.E. Rijtema and H. Wassink (ed.) Water in the unsaturated zone. Proc. Wageningen Symp. June 1966.
Vol. 1. IASAH, Gentbrugge, Belgium.
Wooding, R.A. 1968. Steady infiltration from large shallow circular pond. Water Resour. Res.
4:1259–1273.
Zachmann, D.W., P.C. Duchateau, and A. Klute. 1981. The calibration of Richards flow equation for adraining column by parameter identification. Soil Sci. Soc. Am. J. 45:1012–1015.
Zachmann, D.W., P.C. Duchateau, and A. Klute. 1982. Simultaneous approximation of water capacity
and soil hydraulic conductivity by parameter identification. Soil Sci. 134:157–163.
Zijlstra, J., and J.H. Dane. 1996. Identification of hydraulic parameters in layered soils based on a quasiNewton method. J. Hydrol. 181:233–250.
Zurmühl, T. 1998. Capability of convection dispersion transport models to predict transient water andsolute movement in undisturbed soil columns. J. Contam. Hydrol. 30:99–126.
Zurmühl, T., and W. Durner. 1998. Determination of parameters for bimodal hydraulic functions by inverse modeling. Soil Sci. Soc. Am. J. 62:874–880.