Stochastic Fusion of Information for Characterizing and Monitoring the Vadose Zone
Tian-Chyi J. Yeh* and Jirka Šimůnek
ABSTRACT

are still often raised regarding parameter identifiability
and their uniqueness for particular methods.
While various laboratory and field methods for evaluating soil hydraulic properties are relatively well established, several major problems remain. Most laboratory
methods are applied to samples ranging from 100 toabout 500 cm3. The scale of field methods generally does
not extend beyond a plot of 1 m2 and depths of one toseveral meters. There is an urgent need to develop methods that characterize hydraulic properties of the vadosezone on a much larger scale. Recently developed geophysical methods such as electrical resistivity tomography (ERT) (Daily et al., 1992; Zhou et al., 2001) andground-penetrating radar (GPR) (Binley et al., 2001),
coupled with geostatistically based inverse methods (Liu,
2001; Yeh et al., 2002) appear to be promising tools formonitoring large-scale vadose zones. Combining these
new tools with a three-dimensional hydrologic inverse
model (e.g., Hughson and Yeh, 2000), a detailed threedimensional characterization of the vadose zone on alarge scale appears to become possible.
In this review we examine conditions required tomake inverse problems well posed for unsaturated flow
through homogeneous and heterogeneous soil columns.
Stochastic conceptualizations of inverse problems forthe vadose zone are then introduced. Simple examples
are used to illustrate the principles of geostatistically
based inverse approaches (GA), including cokriging, asequential linear estimator, and a successive linear estimator method. We discuss applications related to GA to
vadose zone inversion problems, hydraulic tomography,
and electrical resistivity tomography for vadose zone
characterization and monitoring. Finally, we introduce
a stochastic fusion of information concept to assimilate
information from soil physics, hydrology, geophysics,
and geology for characterizing and monitoring the vadose zone. Preliminary results of the information fusion
technology are presented. Our discussions, we hope,
will lead to better-designed laboratory and field experiments, as well as to vigorous research and developmentof integrative inversion approaches for characterizing
and monitoring the vadose zone.

Inverse problems for hydrological processes in the vadose zone
are often perceived as being ill posed and intractable. Consequently,
solutions to the inverse problems are frequently subject to skepticism.
In this paper, we examine the necessary and sufficient conditions forthe inverse problems to be well posed and discuss difficulties associated with solving the inverse problems. We subsequently explain theneed for a stochastic conceptualization of inverse problems of thevadose zone. Principles of geostatistically based inverse approaches,
which rely on stochastic concepts, are then illustrated, including cokriging, a sequential linear estimator, and a successive linear estimator.
We then discuss applications involved in the approaches to classicalvadose zone inversion problems (using observed pressure heads, moisture contents, concentrations, and arrival times), hydraulic tomography, and electrical resistivity tomography for vadose zone characterization and monitoring. Finally, we present a stochastic information
fusion technology that assimilates information from both unsaturated
hydraulic tomography and electrical resistivity tomography. Preliminary results suggest that this fusion technology is a promising toolfor effectively characterizing heterogeneity, monitoring processes inthe vadose zone, and quantifying uncertainties associated with vadosezone characterization and monitoring.

K

nowledge of unsaturated soil hydraulic properties is indispensable for successful predictions ofwater flow and solute transport in the vadose zone. A variety of laboratory and field methods are currently available for direct and indirect evaluations of unsaturated
hydraulic conductivity, K, as a function of pressure head,h, and/or water content, ␪. Popular laboratory methods
include the one- and multistep outflow methods (Koolet al., 1985; van Dam et al., 1994), the upward infiltration method (Hudson et al., 1996), and the evaporation
method (Gardner and Miklich, 1962, Šimůnek et al.,1998). Popular field methods include the instantaneous
profile method (Hillel et al., 1972), various unit-gradienttype approaches, sorptivity methods following ponded
infiltration, and the crust method based on steady water
flow. While existing field methods are relatively simple
in concept, these direct measurement methods have anumber of limitations that restrict their use in practice.
For example, most methods are very time-consuming
to execute because of the need to adhere to relatively
restrictive initial and boundary conditions. This is especially true for field gravity drainage experiments involving medium- and fine-textured soils. While most of theabove methods are widely used and accepted, questions

The Inverse Problem in Subsurface Hydrology:
An Ill- or Well-Posed Problem
Movement of water in variably saturated porous media, assuming isothermal conditions and neglecting the

T.J. Yeh, Department of Hydrology and Water Resources, The University of Arizona, Tucson, AZ 85721; J. Šimůnek, George E. Brown,
Jr. Salinity Laboratory, USDA-ARS, Riverside, CA 92507-4617. Received 12 Mar. 2002. *Corresponding author (yeh@hwr.arizona.edu).

Abbreviations: ERT, electrical resistivity tomography; GA, geostatistically based inverse approaches; GPR, ground-penetrating radar;
MAP, nonlinear maximum a posteriori [method]; SLE, successive linear estimator [method]; VG, van Genuchten model (see Eq. [2] and [3]).

Published in Vadose Zone Journal 1:207–221 (2002).

207

208

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

gas phase, may be described using a modified form of
the Richards equation
⳵h ⳵␪
␻Ss
⫹
⫽ 䉮·[K(h)䉮(h ⫹ z)]
⳵t
⳵t

[1]

where h is the pressure head, which is positive whenthe medium is fully saturated and negative when themedium is partially saturated, t is time, ␪ represents thevolumetric moisture content, and z is the positive upward vertical coordinate. The term Ss represents specificstorage and ␻ a transitioning parameter that is one whenh ⱖ 0 and zero when h is negative. To describe the ␪–hrelationship of unsaturated media, van Genuchten’s (1980)
model is often assumed:⫺m␪(h) ⫽ ␪r ⫹ (␪s ⫺ ␪r)共1 ⫹ ␣兩h兩n兲

[2]where 兩 兩 is the absolute value, ␪s is the saturated moisture content, ␪r is the moisture content at residual saturation and ␣, n, and m are shape-fitting parameters withm ⫽ 1 ⫺ 1/n. We further assume that the unsaturated
hydraulic conductivity function, K(h), follows Mualem’spore-size distribution model (van Genuchten, 1980);
that is,K(h) ⫽ Ks(1 ⫺ (␣兩h兩)n⫺1[1 ⫹ (␣兩h兩)n]⫺m)2[1 ⫹ (␣兩h兩)n]m/2[3]where Ks is the locally isotropic saturated hydraulic conductivity. Hereafter, Eq. [2] and [3] are referred to as
the VG model.
Generally, the boundary conditions associated with
Eq. [1] are given as: (i) K(h)䉮φ ⫽ n·q* at boundaries,⍀ 1 , where n is a normal unit vector and q* a specified
specific discharge, and (ii) φ(x, y, z) ⫽ φ* at boundaries,⍀ 2 , where φ* is the specified hydraulic head. The initialcondition is given as φ(x, y, z, t) ⫽ φ*(x, y, z, 0) at t ⫽ 0.
If we define Ks, Ss, ␣, n, ␪s, and ␪r as parameters or
primary variables, then φ, h, and ␪ are state variables,
secondary variables, or system responses. A forward problem refers to solving the flow equation for the pressure
head or moisture distribution in time and space withknown primary variables, and for given initial and boundary conditions. On the other hand, an inverse problem
refers to estimating values of the primary variables frominformation about excitations to the system and its response (secondary information) to those excitations.
A forward problem is well posed if the parameters
and initial and boundary conditions are completely specified in the solution domain so that the problem canhave a unique solution. It is ill posed and has an infinite
number of solutions otherwise. A well-posed problem,
however, does not necessarily warrant a solution. For
instance, a well-posed forward problem for variably saturated flow may still encounter convergence and stability issues related to solution techniques. In addition, awell-posed problem with incorrect information can leadto an erroneous solution. Similarly, an inverse problem
is ill posed if there is no unique solution to the inverse
problem. Next, we will examine conditions necessary for
an inverse problem to be well posed for both steadystate and transient unsaturated flow in porous media.

Steady-State Unsaturated Flow
The fact that the hydraulic conductivity is a functionof the pressure head or moisture content generally increases the difficulties of inverse problems for flow
through unsaturated porous media. Consider steady onedimensional vertical flow in an unsaturated soil column
and assume Darcy’s Law to be valid for the flow process:q ⫽ ⫺K(h)[dh/dz ⫹ 1]

[4]

In this case, the pressure head, h, along the column mustbe specified in addition to the specific discharge so thatthe unsaturated hydraulic conductivity at the given pressure head can be determined. Based on this principle, aunit gradient approach is frequently used to determine
the unsaturated hydraulic conductivity. Since hydraulicconductivity varies with pressure head, several unit gradient situations with different steady-state flow rates must
be created to determine accurately the shape of the hydraulic conductivity function. The number of steady flows
can be reduced if we prescribe a relationship between K
and h. For example, if we assume that K(h) ⫽ Ks exp(␤h),
where Ks was previously defined and ␤ is a pore-size
distribution parameter (Gardner, 1958), we haveq ⫽ ⫺Ks exp(␤h)[5]

In this case, two independent equations are required touniquely determine the value of Ks and ␤, implying thatat least two steady unit-gradient experiments under different pressure head values must be undertaken. The
problem will otherwise be ill posed. Similarly, if one uses
the VG model (Eq. [3]), which has three parameters (Ks,␣, and n), three unit gradient conditions with different
flow rates must be implemented such that three independent equations exist to yield a unique solution. Although the inverse problem becomes well posed, the resultant nonlinear equations may be still difficult to solve.
Consider one-dimensional steady unsaturated infiltration into a column consisting of two layers of known
thickness. The unsaturated hydraulic conductivity ofeach layer is assumed to be described by the exponential
model, and the thickness of the layer is constant, ⌬z.
Therefore, we now have two equations:
Ksi exp[␤ih(z)]⫹ 1冥 ⫽ ⫺q
冤dh(z)dzz 僆 (zi, zi⫹1), i ⫽ 1,2[6]Each of these equations has two unknowns (Ks and␤) if q and h are specified. For unsaturated flow, hvaries nonlinearly within each material. Accordingly, adetailed h distribution for one steady flow rate is necessary to allow specification of h and dh/dz values atdifferent z values, such that two independent equations
can be formulated for each material. If only one h anddh/dz (e.g., unit gradient conditions if they exist) in eachmaterial is known, imposing two different flow rates isessential so that two independent equations for eachmaterial exist. To demonstrate this concept, Fig. 1A and1B show both the true and estimated values of Ks and␣ (after removing their means) in a synthetic soil profile
with 20 layers, respectively. The estimated values were209

www.vadosezonejournal.org

Fig. 1. Comparisons of true and estimated (A) Ks and (B) ␣ values (after removing their means) from one flux boundary condition. (C) and(D) compare the true and estimated values for two flux boundary conditions.

obtained using the successive linear estimator approach(discussed below) with one pressure head measurementin each layer for one steady unsaturated flux situation.
This inverse problem is ill posed since each layer has
two unknowns, Ks and ␣, but only one pressure headvalue, and thus only one independent equation. The
estimates therefore deviate from, yet resemble, the true
fields. If one fully saturated steady and one unsaturated
steady flux experiments are conducted, two differentpressure head and gradient values—thus two independent equations—become available for each layer. The
inverse problem then becomes well posed. Figure 1C
and 1D show the estimated Ks and ␣ values with thetwo fluxes. Notice that the estimates are identical to thetrue values—the inverse problem is well posed. The
parameter n in the VG model was assumed to be knowna priori in this example.
Transient Unsaturated Flow
For transient unsaturated flow problems, after integrating Richards equation (Eq. [1]) over the length of
the soil column (from z ⫽ 0, the bottom of the column,
to z ⫽ L, the top of the column), we have

⫺K(h)

冤

冥兩

⳵h
⳵
⫹1
⫽
⳵z
z⫽0
⳵t

冤

冥

冮 ␪(z, t) ⫺ K(h)冤⳵z ⫹ 1冥兩

z⫽L

z⫽0

⳵h

z⫽L

[7]

The term on the left-hand side of Eq. [7] is the outflow
evaluated at z ⫽ 0. The first term on the right-handside of Eq. [7] denotes changes in storage, while thesecond term is the specific discharge flowing into thecolumn at z ⫽ L. The sum of the two terms on the righthand side of Eq. [8] yields the specific discharge at theend of the column, q(0, t):⫺ K(h)

冤⳵h⳵z ⫹ 1冥兩

z⫽0⫽ q(0, t)[8]Equation [8] suggests that the pressure head, the pressure gradient, and the discharge at z ⫽ 0 at a specified
time must be known beforehand to uniquely define K(h)at the given pressure head. Again, the minimum numberof times that these measurements must be specified will
depend on the form of the invoked hydraulic conductivity model. For instance, if the VG model with three
unknown parameters (Ks, ␣, and n) is used, at least threedifferent discharge rates and pressure head profiles atthree times are needed to create three independent

210

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

equations. This scenario forms the basis of the wellknown instantaneous profile method (Rose et al., 1965;
Watson, 1966).
Knowledge of the pressure head and moisture content
distributions along the column at two specified times
allows evaluation of the pressure head gradient andthe change in moisture content along the length of thecolumn. However, this is not sufficient to make theinverse problem of Eq. [8] well posed. Since the unsaturated hydraulic conductivities K(h) at z ⫽ 0 and z ⫽ Lremain unknown, the number of unknowns is greaterthan the number of independent equations. Additional
information about the pressure head and moisture content distributions at different times does not help because more unknowns [i.e., K(h) at different h values]
are introduced. Consequently, to resolve the ill-posed
issue, in addition to the head and moisture content distributions at different times, one of the boundary fluxes
in Eq. [7] at different times must be specified. Noticethat if at different times measurements of both the pressure head and the moisture content at a point in a
homogeneous soil column are available, the water retention curve, ␪(h), can be constructed for the medium. If
we assume that the parameter values for the retention
curve are representative of the relative conductivity
curve based on the VG model, the relative conductivity–pressure head relationship can be derived without inverting the flow model.
For transient flow in an unsaturated heterogeneous
medium, the general rule is that the initial and boundary
conditions, and many sets of spatial steady-state head
distributions, and temporal head and moisture distributions during transient flows must be known. The inclusion of two or three fluxes is invaluable for making theinverse problem well posed.
These examples assume that the size of each soil blockis the same and known. If the size of each block varies
and the distribution of the blocks is unknown, the problem becomes more complex. While the previous discussion focused on one-dimensional flow problems, thesame principles are also applicable to multidimensional
flow problems. The hydraulic conductivity anisotropyof the equivalent homogeneous medium, or the anisotropy of each block, may have to be considered if they
are significant. This means that more parameters mayneed to be identified. Therefore, additional information
about the system’s spatial and temporal responses mustbe acquired.
The discussions above lead to the fact that if sufficient
information is available to yield enough independentequations, the inverse problem is always well posed andshould have a unique solution. A well-posed problem,
however, may still fail to yield a solution because of limitations of the adopted solution technique and erroneous
information. We should also emphasize that the discussions are restricted to inversions based on mathematicalmodels, which assume a one-to-one relationship between
the hydraulic property and the hydraulic response of asystem in both forward and inverse operators.
The inversion of hydraulic properties from actual laboratory and field experiments may encounter othercomplications. These complications include representativeness of mathematical models for actual physical processes (model error), inability to measure detailed hand ␪ distributions and fluxes (lack of information),precision of numerical models and computational devices (numerical errors), and noise in the measurements(measurement errors), in addition to consistency in scale
between measurements and the representative elementary volume. Thus, a theoretically well-posed problem,
complicated by these factors, may yield a unique but undesirable solution. In recent investigations, many soil
physicists (Toorman et al., 1992; Zurmühl, 1996; Šimůnek and van Genuchten, 1997; Šimůnek et al., 1998) havedeveloped practical approaches to circumvent these difficulties. Nevertheless, the well-posedness conditions discussed above are necessary to have a unique solution
to the inverse problem. These prerequisites must be considered a priori in the design of field and laboratory
experiments for the inverse problem.

Stochastic Conceptualizationof Inverse Problems
While inverse problems associated with laboratory experiments in general can be well posed if sufficient andnecessary information is collected, virtually all inverse
problems for field experiments are ill posed because theprerequisites are difficult to meet, especially when one
attempts to identify detailed heterogeneity. Many approaches have attempted to solve these types of inverse
problems in subsurface hydrology (e.g., Neuman, 1973;
Yeh, 1986; Sun, 1994; McLaughlin and Townley, 1996).
In spite of their differences, all methods, in theory,
should yield an exact solution if the inverse problem iswell posed. On the contrary, no unique solution will
exist if the problem is ill posed despite having a very
effective solution technique. The fact is that if many
global minima exist, most or all solution techniques will
yield large numbers of solutions to the inverse problem.
However, if some prior information about the conductivity distribution (e.g., the average block size, the maximum and minimum parameter values), in addition tosparsely observed responses of a system, is given, thesolution can be better constrained (e.g., Menke, 1989).
Even with this prior information, the estimated parameters still involve uncertainty because of spatial variability of the parameters.
Because of these considerations, a probabilistic description (or stochastic representation) of the hydraulicproperties becomes appropriate. That is to say, eachof the properties of a geological formation should beconsidered as a stochastic process with an infinite number of possible realizations, characterized by a joint probability distribution (Gelhar, 1993). In practice, the joint
probability distribution is seldom known but can beapproximated by the first and second moments of samples. The first moment, the mean, represents the most
likely value of the property. The second moment, thespatial covariance function, specifies the variance andcorrelation structure of the process, analogous to a description of maximum and minimum values of the prop-

211

www.vadosezonejournal.org

erties and the average block size in previous discussions,
respectively. If one adopts the stochastic representationof the hydraulic property, a corresponding response of
the formation to an excitation is then considered a stochastic process.
With limited secondary information, an inverse model
thus is best perceived as a means to produce property
and response fields that agree with properties and responses at sample locations. In addition, these fields
must satisfy the statistics (i.e., the mean and covariance)
describing their spatial variability, while the governing
equation must describe the underlying physical process.
In a conditional probability concept, this resultant field
is a conditional realization of the property or response
field, among many possible realizations in the ensemble.
While many possible realizations of such a conditional
field exist (i.e., nonunique solutions), the conditionalmean field is unique. This field also represents the most
likely solution to the inverse problem, even though thismay not necessarily be the true field of the soil profile
or geological formation. Its deviation from the true field
is quantified through the conditional variance (uncertainty). As more pieces of primary and secondary information are acquired, the conditional mean will graduallyresemble the true realization of the property field of
the given geological formation, and the uncertainty progressively diminishes. An ill-posed inverse problem is thus
best considered as a stochastic inverse problem, whereas
a well-posed problem is a deterministic inverse problem—a least-squares approach would be appropriate.

Geostatistically Based Inverse Approaches
If one accepts the stochastic conceptualization of theill-posed inversion problem, GA are appropriate inverse
methods. In this section, principles of cokriging, sequential conditioning, and a successive linear estimator approach are illustrated using some simple examples.
Cokriging
Cokriging in essence is a classical linear predictor thatconsiders spatial correlation structures of flow processes(such as pressure heads and velocities) and hydraulicproperties of geological media. In addition, cokriging
takes into account possible cross-correlation between
the flow processes and the hydraulic properties. Cokriging has been widely used to estimate transmissivities,
heads, velocities, and concentrations of pollutants inhighly heterogeneous aquifers (e.g., Kitanidis and Vomvoris, 1983; Hoeksema and Kitanidis, 1984, 1989; Rubin
and Dagan, 1987; Gutjahr and Wilson, 1989; Sun andYeh, 1992; Harvey and Gorelick, 1995; Yeh et al., 1995,
1996). It has also been used to estimate water content
distributions based on combined measurements of watercontent, soil water pressure head, soil surface temperature, and/or soil texture (e.g., Vauclin et al., 1983; Yates
and Warrick, 1987; and Mulla, 1988).
Consider estimating the hydraulic conductivity value,f 0, at location x0, using simple cokriging and a known
hydraulic conductivity, f 1, and pressure head, h2, at locations x1 and x2, respectively; that is,f̂ 0 ⫽ ␭01 f 1 ⫹ ␤02 h2[9]where the weights ␭01 and ␤02, represent the contribution
to the hydraulic conductivity estimate at x0 from the
known hydraulic conductivity, f 1, and the pressure head,h2, at locations x1 and x2, respectively. The weights canbe obtained by minimizing the mean-square error of
the estimate, E[(f 0 ⫺ f̂ 0)2]. The minimization leads to asystem of equations:

冤CC (x(x ,, xx )) CC (x(x ,, xx ))冥 冤␤␭ 冥 ⫽ 冤CC (x(x ,, xx ))冥
ff

1

1

fh

1

2

01

ff

0

1

fh

1

2

hh

2

2

02

fh

0

2

[10]

where Cgg (xi, xj) represents the spatial covariance ofthe stochastic processes, gi and gj (which can be either
conductivity or head) at locations xi and xj. Notice thatCfh(xi, xj) ⫽ Chf (xj, xi). The solution of Eq. [10] yields
the cokriging weights,
␭01 ⫽␳hh(x2, x2)␳ff (x0, x1) ⫺ ␳fh(x1, x2) ␳fh(x0, x2) ,␳hh(x2, x2) ⫺ ␳fh(x1, x2)2

␤02 ⫽

␳fh(x0, x2) ⫺ ␳fh(x1, x2) ␳fh (x0, x1)
␳hh(x2, x2) ⫺ ␳fh(x1, x2)2

[11]

if we define normalized covariance functions, that is,␳ff (xi, xj) ⫽ Cff (xi , xj)/Cff (xi , xi), ␳hh(xi , xj) ⫽ Chh(xi , xj)/
Cff (xi , xi), and ␳fh(xi , xj) ⫽ Cfh(xi , xj)/Cff (xi , xi).
Note from Eq. [11] that the weight does not dependon the variance of hydraulic conductivity. To evaluate
the weights, the covariance, Cff (xi , xj) must be specified
a priori. In theory, the covariances Chh(xi , xj) and Cfh(xi ,xj) can be derived from data, if sufficient data sets areavailable, but are generally calculated using Cff (xi , xj)
and a first-order analysis based on the flow model (e.g.,
Dettinger and Wilson, 1981). As a result, the flow process is implicitly considered in this cokriging technique.
Consider a simple linear forward model for flow, say,h ⫽ a21 f 1 ⫹ a22f 2 , where a21 and a22 are coefficients. We
are to estimate the hydraulic conductivity value, f 2 , atlocation x2 , using simple cokriging with known hydraulicconductivity, f 1 , and hydraulic head, h2, at location x1
and x2 . In this case, cokriging yields an exact solution.
The covariances needed are calculated from the forward
model and they are
Chh(x2, x2) ⫽ a221Cff (x1, x1) ⫹ 2a21a22Cff (x1, x2) ⫹ a222Cff (x2, x2)
Cfh(x1, x2) ⫽ a21Cff (x1, x1) ⫹ a22Cff (x1, x2)
Cfh(x2, x2) ⫽ a21Cff (x1, x2) ⫹ a22Cff (x2, x2)[12]Using these covariances, a solution to a matrix similar to Eq. [10] lead to the weights: ␭21 ⫽ ⫺a21/a22, and␤22 ⫽ 1/a22 . Notice that in this case the weights are
independent of the covariance function of hydraulicconductivity. This is because the problem is well posed(i.e., a deterministic inverse problem). On the other
hand, if the problem is ill posed (i.e., a stochastic inverse
problem), cokriging yields weights that depend on thecovariance functions (see Eq. [11]) and can produce aunique estimate that satisfies the minimum mean squareerror criterion. The influence of the covariance, however, decreases with an increase in the number of observations of either h or f.

212

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

Sequential Linear Estimator

Then the weight is given byInstead of estimating f 0 simultaneously using both f 1
and h2 as in the previous example, we now will use asequential linear estimator approach (Vargas-Guzman
and Yeh, 1999) that estimates f 0 , f 1 , and h2 using only
f 1 at first, and then improves the estimate by adding theobserved h2 information. Our first estimates are␪02 ⫽f̂ (1)f̂ (1)ĥ (1)
0 ⫽ ␭01 f 1 ,1 ⫽ ␭11 f 1 ,2 ⫽ ␤21 f 1[13]

where the superscript (1) indicates the first step of thesequential estimate andC (x , x )C (x , x )␭01 ⫽ ff 0 1 ⫽ ␳ff (x0, x1), ␭11 ⫽ ff 1 1 ⫽ ␳ff (x1, x1) ⫽ 1,
Cff (x1, x1)
Cff (x1, x1)
Cfh(x1, x2)
␤21 ⫽
[14]
⫽ ␳fh(x1, x2)
Cff (x1, x1)

Notice that the estimate of f 1 is exact because kriging
is an exact interpreter at the sample location (i.e., a conditioning approach).
Next, we will update our estimate of f 0 by using thedifference between the observed h2 and its estimate from
Eq. [13]:(1)f̂ (2)
0 ⫽ f̂ 0 ⫹ ␪02(h2 ⫺ ĥ2)[15]where the superscripted (1) and (2) indicate the orderof the sequential linear estimator. The weight, ␪02, in
Eq. [15] is obtained by minimizing the variance of the
new estimate:E[( f 1 ⫺ f̂ 0(2))2] ⫽ E兵[( f 0 ⫺ f̂ 0(1)) ⫺ ␪02(h2 ⫺ ĥ 2(1))]2其(1)
ε ff(2)(x0, x0) ⫽ ε ff(1)(x0, x0) ⫺ 2␪02ε fh(1)(x0, x2) ⫹ ␪ 202ε hh(x2, x2)[16]
where ε represents the conditional (or residual) covariance and cross-covariance. Differentiating Eq. [16] withrespect to ␪02 and setting the result equal to zero, wehave(1)␪02ε (1)hh (x2, x2) ⫽ ε fh (x0, x2)[17]in which the conditional covariance can be determined byE[(h2 ⫺ ĥ 2(1))2] ⫽ E[(h2 ⫺ ␤21f 1)2](1)
ε hh(x2, x2) ⫽ Chh(x2, x2) ⫺ 2␤21Cfh(x1, x2) ⫹ ␤ 221Cff (x1, x1)[18]
Using Eq. [14], the normalized conditional head covariance for Eq. [18] becomes2
␳(1)hh (x2, x2) ⫽ ␳hh(x2, x2) ⫺ ␳fh(x2, x1)[19]Equation [19] shows that the conditional variance of thehydraulic head thus becomes smaller than the unconditional head variance because of the effects of conditioning using f 1. Similarly, the conditional cross-covariancecan be derived:E[( f 0 ⫺ f̂ 0(1))(h2 ⫺ ĥ 2(1))] ⫽ E[( f 0 ⫺ ␭01f 1)(h2 ⫺ ␤21f 1)]
ε fh(1)(x0, x2) ⫽ Cfh(x0, x2) ⫺ ␤21Cff(x0, x1)⫺ ␭01Cfh(x1, x2) ⫹ ␭01␤21Cff (x1, x1)[20]
Written in terms of normalized covariances, Eq. [20]becomes␳(1)fh (x0, x2) ⫽ ␳fh(x0, x2) ⫺ ␳ff (x0, x1) ␳fh(x1, x2)

[21]

ε (1)
␳ (x , x ) ⫺ ␳ff (x0, x1) ␳fh(x1, x2)
fh (x0, x2)⫽ fh 0 2(1)
ε hh (x2, x2)
␳hh(x2, x2) ⫺ ␳2fh(x2, x1)[22]

The estimate (Eq. [15]) thus becomes
f̂

(2)
0

⫽ f̂ (1)0 ⫹ ␪02(h2 ⫺ ĥ2) ⫽ (␭01 ⫺ ␪02␤21) f 1 ⫹ ␪02h2⫽

␳hh(x2, x2) ␳fh(x0, x2) ⫺ ␳ff (x0, x1) ␳fh(x1, x2)
f1
␳hh(x2, x2) ⫺ ␳2fh(x1, x2)⫹

␳fh(x0, x2) ⫺ ␳ff (x0, x1) ␳fh(x1, x2)h2
␳hh(x2, x2) ⫺ ␳2fh(x1, x2)[23]This result is essentially Eq. [9] with the weights given
in Eq. [12]. Therefore, linear estimation using different
information sequentially is equivalent to cokriging thatuses all of the information simultaneously. The sequential approach nevertheless reduces computational efforts. Consider there are Nf measurements of hydraulicconductivity and Nh measurements of the hydraulichead. The computational cost of inverting the matrix in
Eq. [10] is (Nf ⫹ Nh)3, whereas the cost of the sequential
estimator is N f3 ⫹ N h3 . Suppose the measurements are further subdivided into subgroups and then used sequentially in the estimation. The savings in terms of computational efforts can be even more significant. It must beemphasized that during sequential estimation it is essential to update the unconditional covariance after conditioning with available information.
Notice that the weights of the GA method for different estimation locations also can be obtained independently. That is, one can solve Eq. [10] to obtain ␭11
and ␤12 for Location 1, independently from solving theequation for obtaining ␭21 and ␤22 for Location 2. Therefore, the algorithm is most suitable for parallel computing, and is highly efficient, since the inverse of the lefthand side of Eq. [11] has to be computed only once. Thisfeature of the GA approach offers a similar advantage
to inverse modeling of variably saturated flows where
several parameters (e.g., Ks, ␣, n, ␪s, and ␪r of the VG
model) require estimation—the GA can estimate theseparameters independently.
In the case of stochastic inverse problems where thenumber of observations, M, is always much smaller than
the number of parameters to be estimated, N, which isthe scenario in real-world problems, the GA method
offers an even greater advantage over other classicalinverse methods. For inverse problems of variably saturated flow, the GA also permits sequential inclusionof pressure head and moisture content measurements.
Moreover, the GA method allows sequential inclusionof observations from multiple experiments such as inhydraulic tomography and sequential fusion of information of different processes, as discussed below.
Successive Linear Estimator
While the GA approach is powerful, it still is only
a linear predictor. Relationships between primary andsecondary variables of the vadose zone are generally
highly nonlinear. Therefore, the GA approach cannot
fully exploit available secondary information. To over-www.vadosezonejournal.org

come this limitation, Zhang and Yeh (1997) adapted a
successive linear estimator (SLE) technique (Yeh et al.,1996) to the vadose zone inverse problem. The principleof the SLE is similar to the sequential linear estimator.
However, instead of using different pieces of information sequentially, SLE uses the same information successively to account for the nonlinear relation, and thus to
improve the estimate.
Detailed SLE algorithms for vadose zone inverse
problems can be found in Hughson and Yeh (2000).
Here we present a brief outline of the steps involved.
The SLE consists of seven steps. Step 1 starts with cokriging, which integrates primary and secondary information to estimate the value of the primary variable
at locations where no information of the variable isavailable. In Step 2, the covariance of the primary variable is updated to reflect the effects of available information. In Step 3, the newly estimated variable fieldfrom Step 1 is used to simulate the hydraulic head bymeans of a forward flow model. In Step 4, the conditional covariances of the secondary information andtheir cross-covariances with the primary values are modified using a first-order analysis as in the sequential
linear estimator approach. For Step 5, these newly evaluated covariances and cross-covariances are used tocompute new weights. Step 6 is where the new weights,
along with the difference between simulated and observed system responses at observation locations, areused to improve the estimate of the primary variable,
which is similar to the sequential linear estimator approach. In Step 7, the weights are used to update conditional covariances for the next iteration. This newly
updated primary variable field and new conditional covariances are used again in Steps 3 and 4, followed by
Steps 5 through 7. In essence, Steps 3 through 7 are
repeated until no improvement in the estimate of the
primary variable is found (i.e., the solution converges).
Similar to a nonlinear maximum a posteriori (MAP)
method (see McLaughlin and Townley, 1996), the SLE
is also based on a Bayesian framework but differs from
MAP in many aspects (Kitanidis, 1986, 1997). More
importantly, during nonlinear iterations, our SLE updates the second moment of the posteriori probabilityof the primary variable in a consistent manner using afirst-order approximation approach. For a well-posed
inverse problem (saturated or unsaturated flow), theSLE converges to the true solution rapidly (as shownin Fig. 1). Also, notice that the SLE can account formeasurement errors in the estimation by addition oferror variances to the diagonal of the left-hand side ofmatrix of (10) or (17).

Applications
Unsaturated Inverse Problems
Traditional inverse modeling approaches typicallyconsider secondary information collected during one excitation to the system only (e.g., an infiltration event)to identify all primary variables. Based on this approach,a cokriging technique was developed by Harter and Yeh

213(1996) to identify unsaturated hydraulic parameters inheterogeneous vadose zones during steady-state uniform flow and by Yeh and Zhang (1996) for nonuniform
flow. Later, Hughson and Yeh (1998) extended the approach to transient unsaturated flow. Harter and Yeh(1996) showed that measurements of hydraulic headenhance the prediction of the flow paths of solutes inthe vadose zone. Yeh and Zhang (1996) reported thatthe unsaturated hydraulic parameters of a heterogeneous vadose zone could be identified if sufficient pressure head and degree of saturation information wasavailable. They found that, under relatively wet conditions, measurements of the pressure head improved estimates of the saturated hydraulic conductivity. Under
relatively dry conditions, measurements of saturation
enhanced the estimate of the pore-size distribution parameter, ␤, of the Gardner model.
Li and Yeh (1999) used cokriging to estimate the hydraulic conductivity from pressure head, solute concentration, and solute arrival time measurements in ahypothetical, heterogeneous vadose zone subject tosteady-state infiltration at different degrees of saturation. Their analysis showed that the performance of
cokriging deteriorates as the medium becomes less saturated when either head or concentration measurements
are used. They attributed this result to an increase in
nonlinearity between head or concentration and the
conductivity as the medium becomes less saturated, andto the linear predictor nature of cokriging.
Among pressure head, solute concentration, and solute arrival time measurements, Li and Yeh (1999) foundthat pressure head measurements of steady-state flow
fields were the most useful secondary information for
estimating the Ks field using cokriging. They attributed
this finding to several factors. First, the nonlinear relationship between head and Ks may have been relatively
mild for the cases they studied. Additionally, the assumption of ergodicity was approximately satisfied forsteady-state flow. In other words, the ensemble meanhead distribution evaluated with mean values of theparameters closely approximates the spatial mean head
distribution. Such an existence of ergodicity reduces thevariance of the head, and consequently improves thelinearity between head and conductivity. Conversely,
the ergodicity assumption cannot be easily satisfied forthe solute transport case. Only when a solute plume has
traveled over enough correlation lengths, has sampledenough heterogeneity, and has become a Fickian process, will the ergodicity assumption be satisfied. Becauseof a lack of ergodicity, the variance of a concentration
perturbation can be very large and the cokriging estimation using concentration measurements can be unsatisfactory.
In theory, concentration distributions of a tracer are
the result of both convection and dispersion processes.
Suppose we neglect dispersion and assume that convection is the dominant process. The convection process isa function of not only the hydraulic conductivity but
also the hydraulic gradient and the moisture content.
Cross-correlation between the concentration and thehydraulic conductivity is thus not as strong as the rela-

214

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

tionship between pressure head and hydraulic conductivity. This suggests that for estimating the hydraulicconductivity, concentration measurements are generally
less effective. The same argument applies to arrival timesof a tracer concentration. However, Li and Yeh (1999)
found that estimates of the conductivity based on headmeasurements could still be improved by incorporating
additional concentration information.
Of course, the examples assumed that the scales ofmeasurements of hydraulic properties, pressure head,
concentration, and arrival time are the same as the element size of the model. In practice, this assumption will
unlikely be met and thus, results in biased estimates.
While issues of different scales remain to be solved, apractical solution is to collect “enough” measurementsin space (e.g., along a borehole), smooth them with amoving average method based on the size of the element, and then use the smoothed values at measurement
locations as the conditioning values. Alternatively, one
can add nugget effects to the covariance functions forthe measurements in the SLE such that contributions
from the measurements to estimates are reduced. As aconsequence, estimated fields are expected to be smooth,but they may retain essential heterogeneity patterns.
Lastly, one must recognize that the propagation processof pressure head is highly diffusive, moisture content isless, and concentration and travel time are the least. Therefore, effects of the scale are expected to be minimal if
the pressure head information is used.
Because of the highly nonlinear nature of most vadosezone inverse problems, Zhang and Yeh (1997) showedthat the SLE approach could yield more detailed imagesof unsaturated hydraulic parameters than cokriging.
This implies that the SLE, which successively incorporates the nonlinear relation between the primary andsecondary variables, can maximize the utility of available secondary information.
The SLE principle was adopted by Hughson and Yeh(2000) to develop an inverse model for three-dimensional, transient flow in heterogeneous vadose zones.
They also used the sequential estimator approach toallow pressure head and water content data obtainedat different times to be sequentially included in the
inversion. Using this inverse model, they investigatedthe efficacy of estimating the VG parameters using pressure and moisture content measurements at relatively
early, intermediate, and late time periods. They concluded that for the cases investigated, late time data
provide the best estimates.
Hydraulic Tomography
Results of the previous applications of the GA methods to vadose zone inverse problems demonstrate a
simple fact: in spite of the robustness of GA based
inverse models, accurate identification of a large number of primary variables in a three-dimensional vadosezone still requires a large amount of secondary information. To satisfy this requirement using traditional sampling means, installation of numerous soil sampling sites
and/or boreholes is necessary, except that such opera-

tions are often impractical and cost-prohibitive. Collecting copious secondary information with a minimum
number of boreholes and a limited budget then becomes
an important and practical issue for improving the inverse solution.
Hydraulic or pneumatic tomography appears to be
one possible technique to accomplish this goal. Hydraulic or pneumatic tomography can generate many pressure measurements with a fixed number of boreholes(Gottlieb and Dietrich, 1995; Butler et al., 1999; Yehand Liu, 2000; Liu et al., 2002; Vesselinov et al., 2001a,b).
For example, using packers, two fully screened wells in
an aquifer can be partitioned into many intervals. By
sequentially pumping water at different intervals, andmonitoring the steady-state head response at the otherwells and packed-off intervals, many pairs of head–
discharge data sets can be obtained with only two wells.
This vast amount of information may yield additional
independent equations for the inverse process—analogous to a multiple pumping test to determine aquifer
anisotropy (Neuman et al., 1984).
Interpretation of the vast amount of data from tomographic tests demands an efficient inverse algorithm.
A combined SLE and sequential method offers great
promise to accumulate the high-density secondary information from the tomographic test and maintains thesystem of equations at a manageable size so that they
can be solved with the least numerical difficulties. Using
this sequential SLE approach, Yeh and Liu (2000)
showed that tomographic aquifer tests have significant
advantages over traditional pumping tests. Such tests
can provide a detailed image of three-dimensional aquifer heterogeneity with the same number of wells used in
traditional pumping tests. They also examined network
design issues and reported that the effects of the spatial
covariance function diminish as the number of hydraulictomography data sets increases. They nonetheless emphasized accuracy of hydraulic head measurements.
Subsequently, Liu et al. (2002) conducted sand box experiments to test the hydraulic tomography approachand concluded that the tomography is a viable technology for delineate aquifer heterogeneity and the sequential SLE algorithm is a promising tool for inverting thevast amount of data generated with tomography tests.
While their investigations dealt with saturated flow
problems, the tomographic concept can be equally wellapplied to unsaturated infiltration tests, that is, unsaturated hydraulic tomography. The unsaturated hydraulictomography method involves the sequential injectionof water at different depths and locations, followed bymonitoring the response of the vadose zone at selected
locations. Using an appropriate inverse modeling technique, the vast amount of secondary information canthus be processed to image spatial distributions of the
primary variables of the vadose zone.
Pumping or injection air tests in the unsaturated zone
are one of many possible approaches to determine thepermeability and porosity of unsaturated soil and fractured rocks. In principle, air injection tests are very
similar to their hydraulic counterparts conducted in fullywater-saturated media. Air is either injected into or

215

www.vadosezonejournal.org

withdrawn from sections of boreholes isolated by meansof inflatable packers, with the pressure responses monitored in observation wells and packed-off intervals. The
pressure response in observation wells can be related
to pneumatic flow parameters such as permeability andporosity through analytical techniques or numerical inverse modeling. Pneumatic tomography (i.e., cross-holepneumatic injection tests conducted sequentially andcorresponding inverse interpretation) has recently been
proposed as a method for characterizing subsurface heterogeneity (Vesselinov et al., 2001a,b). In general, pneumatic injection and gaseous tracer experiments in fractured rocks are not widely performed. Much of current
experience has been gained during pneumatic injection
tests in tuffs at Yucca Mountain, Nevada (LeCain andWalker, 1994; LeCain, 1996, 1998; Wang et al., 1998;
Huang et al., 1999), in Box Canyon, Idaho (Benito etal., 1998, 1999), and at the Apache Leap Research Sitenear Superior, AZ (Trautz 1984; Yeh et al., 1988; Rasmussen et al., 1990, 1993; Guzman et al., 1996; Illman
et al., 1998).
Despite their ability to create many data sets to constrain the inverse solution, hydraulic and pneumatic tomography also have several limitations. The effectiveness of tomography has been found to decrease rapidly,
indicating that excessive sequential excitation often produces only redundant information (Yeh and Liu, 2000)
because secondary information is always collected at
the same locations. In addition, data are often strongly
affected by barometric pressure fluctuations (Illman etal., 1998; Illman and Neuman, 2001), which cause theeffectiveness of the data to decline rapidly. As a result,acquisition of high-density and accurate secondary information throughout the vadose zone remains the only
viable means to enhance the ability of an inverse model
to produce high-resolution subsurface images.

zone can be monitored. Furthermore, assuming universal moisture content and resistivity relations (i.e., Archie’s Law), a three-dimensional image of the moisture
content distribution over a large volume of geologicalmaterial can thus be obtained.
Electrical resistivity tomography also relies on inversion of the potential equation. Yeh et al. (2002) extended the sequential SLE algorithm to three-dimensional ERT inversions such that both electric potential
and point measurements of the electrical conductivity
can be included in the inversion. They showed that electric potential sampling arrays perpendicular to bedding
are more effective than those parallel to bedding. They
also showed that the ERT can detect the general patternof change in the resistivity and, thus, the pattern ofchange in the moisture content. Furthermore, the resolution of the ERT image can be easily improved by deploying a large number of sensors within boreholes andat the soil surface. The relative inexpensiveness of thesensor makes such a deployment feasible, thus making
ERT a highly desirable monitoring tool for vadose zone
investigations. However, recent field studies (Baker,
2001; Yeh et al., 2002) indicated that parameters ofArchie’s Law exhibit profound spatial variability, although they are spatially correlated. This variability
compounds the difficulty in translating resistivity tomoisture content. In other words, at some location, asmall change in resistivity may indicate a large change
in water content, while at other locations a substantialchange in the resistivity may correspond to only a smallchange in the water content (Yeh et al., 2002). Consequently, the ERT’s ability to yield an accurate imageof change in moisture content, or more critically the
moisture content itself, remains to be proven.

Stochastic Fusion of Information
Electrical Resistivity Tomography
To attain high-density secondary information needed
for hydrological inverse modeling, geophysical surveys
appear to be viable and cost-effective technologies.
Whereas geophysical surveys may not provide the primary information for hydrologic modeling, they can be
extremely cost-effective, indirect tools for detecting geological structures (e.g., Rubin et al., 1992; Hyndman and
Gorelick, 1996; Rea and Knight, 1998; Hubbard et al.,1999, 2001) and for monitoring hydrological processes.
Recently, ERT surveys have been demonstrated to bepractical tools for collecting high-density moisture content data in the vadose zone without excessive invasive
sampling (e.g., Daily et al., 1992, Zhou et al., 2001,
Brainard et al., 2001).
Similar to hydraulic tomography, ERT emits DC currents at a point in space and monitors the electricalpotential at other locations. By moving the DC source
location, one can generate many electrical potential fields
and source pairs, from which a three-dimensional imageof changes in the resistivity can be derived using an inverse model. Assuming changes in resistivity mimicchanges in moisture content, water flow in the vadose

Our previous inverse modeling efforts have brought
about several aspects of characterization and monitoring of vadose zones:
1. Geologic materials exhibit random spatial variations in hydro- or geophysical properties on a multiplicity of scales. Small-scale direct sampling of
the entire vadose zone is in most or all cases virtually impossible.
2. Reliable inversions of variably saturated flow
problems require a large amount of secondary information at high densities over the entire volumeof the investigation, and many constraints prohibit
the collection of such information.
3. Tomography approaches can produce additional
data sets to constrain the inversion, but their efficacy can be limited.
4. Geophysical surveys have the potential to collecta vast amount of information pertinent to hydrological inversion; however, significant spatial variation in the relations between various hydrologic
and geophysical properties diminishes the value ofthe information for its intended applications.

216

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

Approach
Because of these problems and uncertainties, assimilation of different types of information using a stochastic
approach appears to be the only plausible solution. In
other words, while geological characterization, point measurements of hydrological and geophysical properties,
hydraulic tomography, and ERT have their pro andcons, an integration of their individual strengths may
facilitate a better way to monitor and characterize thevadose zone and quantify uncertainty.
During the integration, not only the hydrological andgeophysical inversion methods need to include available
primary and secondary information to better condition
the ensemble of their primary variables, but the reciprocal nature of hydrologic and geophysical informationand inversions must be recognized as well. For example,hydrologic information can provide useful constraints
for an ERT inversion, while, on the other hand, an ERT
inversion can furnish a vast amount of water content
data for hydrological inverse modeling. Such a reciprocal relation between hydrologic and geophysical information and inversions thus demands a joint inversion
that requires an iterative approach to fully utilize allavailable information. Consequently, conditional means,the best and unique estimates with minimal uncertainty,
can be obtained and their conditional variances canquantify their uncertainty. This is the basis of the stochastic fusion of information concept described below.
The stochastic fusion approach comprises two levels.
Level 1 fusion aims to include information related tounsaturated hydraulic tomography or ERT surveys toindependently enhance interpretations of spatial distributions of their own primary variables. For this purpose,the sequential SLE algorithm is most suitable.
The information required for inversions of unsaturated hydraulic tomography consists of spatial covariance structures of the primary variables, and hydraulicproperties from core samples, as well as in situ pressure
head and moisture content measurements during the
tomography. For the electrical tomography, information
to be fused may include electric potential measurementsof the electric field induced by the ERT survey, point
measurements of resistivity and moisture content, andparameters of resistivity–moisture relations. In addition,
the prior information about covariance structures ofresistivity, moisture content, and parameters of resistivity–moisture relation will be included.
Outputs from this level of fusion of hydraulic tomography are conditional means of hydraulic parameters,pressure head, and moisture content distributions andtheir associated covariance functions. For the ERT, conditional mean resistivity and moisture content fields andmean parameter fields of the resistivity–moisture content relation, in addition to their conditional covariancestructures, are products of the Level 1 fusion.
Level 2 fusion aims to honor the reciprocal natureof hydraulic and ERT inversion. It takes an iterative
approach to arrive at the best estimate of the primaryvariable fields and uncertainty. The flow chart shownin Fig. 2 depicts a general concept of the fusion process.

Fig. 2. Flow chart of Level 2 iterative stochastic fusion.

Specifically, during an unsaturated hydraulic tomography experiment, several water infiltration tests at different locations (packed-off intervals) could be conductedsequentially. During each test, some point measurements of the pressure head and moisture content are
taken, while an ERT is deployed to monitor water movement. The resistivity image from the ERT is improvedby using Level l fusion to integrate measured hydrologic
and geophysical information into the ERT inversion
process along with some prior information. The priorinformation, the mean and covariance structure of theresistivity and moisture content fields during the unsaturated hydraulic tomography test, can be derived froma forward simulation of the hydraulic tomography testwith mean hydraulic parameters. With this information,
the ERT inversion can yield a reasonable result and isable to estimate not only changes in resistivity and thewater content, but also the moisture content itself andtheir conditional moments.
Estimates of the moisture content distributions, together with their conditional moments, from the inversion of the ERT are subsequently used as input to thehydraulic tomography inversion model, which estimates
distributions of the unsaturated hydraulic parameters, aswell as their conditional moments. With the additional
information from the ERT regarding moisture content
at locations where no samples were available, the hydrological inversion of unsaturated media properties thus
improves, and a more detailed image of the hydraulicwww.vadosezonejournal.org

217

Fig. 3. True (A) ln ␣ and (B) ␪ fields of the synthetic vadose zone.

heterogeneity emerges. Simulations using the newly derived heterogeneity image then produce better pressure
head and moisture content estimates at locations where
no moisture content measurements are available. The
improved moisture content simulations and their covariances are subsequently returned to an ERT inversion toenhance estimates of resistivity and, in turn, moisture
content distributions. The iterative process continues until no further improvement in hydraulic heterogeneity,moisture content, and resistivity estimates can be gained.
This iterative procedure is repeated for each unsaturated hydraulic tomographic test. With such an iterative
fusion of information, data collected from both hydraulic and geophysical tests can then be fully used. Unsaturated hydraulic tomography may thus be a viable andcost-effective technology for characterizing the vadosezone. Meanwhile, an ERT may become both a reliable
and cost-effective monitoring tool for moisture contents.
ExampleTo demonstrate the promise of the stochastic fusion
technology, some preliminary results are presented below. During the past few years, we have developed amethod that fuses information obtained from ERT surveys, point hydrological measurements, and unsaturated
hydraulic inversion modeling to better characterize andmonitor the vadose zone. The effectiveness of this information fusion technology is illustrated below using a
synthetic vadose zone. This synthetic vadose zone consisted of 200 elements, with each element having a dimension of 20 cm in both horizontal directions and 10
cm in the vertical. The unsaturated hydraulic propertiesof each element were assumed to be described by the
VG model. The variability in ␪s, and ␪r , were assumed

negligible; both were treated as deterministic constants
with values of 0.366 and 0.029, respectively. The parameters, Ks , ␣, and n, were considered as random fields withgeometric means of 0.0063 cm s⫺1, 0.028 1/cm, and 2.0,
respectively. The variances of lnKs , ln ␣, and ln n were0.1, 0.1, and 0.01, respectively. It was also assumed that
all three parameters possessed the same exponential
covariance function with a horizontal correlation scaleof 240 cm and a vertical correlation scale of 20 cm.
Following the generation of random hydraulic parameter fields, a hydrostatic negative pressure head distribution, with zero pressure head at the bottom, was assigned
to the vadose zone as the initial condition. Next, a steady
infiltration event was simulated using a finite element
flow model (Srivastava and Yeh, 1992). The top centerof the vadose zone (from x ⫽ 80 to 120 cm, y ⫽ 0 to
20 cm, and z ⫽ 200 cm) was treated as a constant headboundary with a pressure head of ⫺80 cm. The remainder of the surface and the two sides of the domain
were considered as no-flux boundaries; the bottom wasassumed to be a water table. Figure 3A depicts the ln␣ field for the vadose zone, while the simulated “true”
␪ field corresponding to the steady infiltration is shownin Fig. 3B.
Two cases were examined in which the ln ␣ fields
were estimated using ␪ information derived from two
different approaches. In both cases, the other unsaturated hydraulic parameter fields were assumed known
for the vadose zone. In Case 1, point measurements of␪ were assumed available at six locations (black circles
in Fig. 3B). Using a theoretical moisture content covariance, the six measurements, and a kriging technique, anestimated ␪ field (Fig. 4A) was derived. The theoretical
moisture content covariance was calculated using a first-

218

VADOSE ZONE J., VOL. 1, NOVEMBER 2002

Fig. 4. (A) Estimated ␪ fields using six ␪ measurements and a kriging technique; (B) estimated ln ␣ field from an unsaturated inverse modelusing 34 ␪ values from Fig. 4A; (C) simulated ␪ fields using the ln ␣ field from Fig. 4B. (D) Estimated ␪ fields using six ␪ measurements
and an ERT survey; (E) estimated ln ␣ field from an unsaturated inverse model using 34 ␪ values from Fig. 4D; (F) simulated ␪ fields theusing the ln ␣ field from Fig. 4E.

219

www.vadosezonejournal.org

order analysis of the flow model with mean values andcovariances of the parameters. From this kriged ␪ field,
34 ␪ values were subsequently sampled at locations indicated by red circles. Afterward, our hydraulic inversion
model (Hughson and Yeh, 2000) was employed using
these ␪ measurements in conjunction with one point
measurement of the ln ␣ parameter (red circle in Fig.
4A) to estimate the ln ␣ field over the entire vadosezone. With this estimated ln ␣ field, a forward simulation
yielded a new ␪ distribution corresponding to steady
infiltration (Fig. 4C).
In Case 2, an ERT survey was conducted over thevadose zone with a synthetic resistivity field that wasconverted from the “true” ␪ field assuming a spatiallyconstant resistivity–moisture relation. A geostatistically
based ERT inversion model for ␪ (Liu, 2001) was thenused in which the six ␪ measurements together with thecovariance function of ␪ were used to constrain the
inversion. Figure 4D shows the estimated ␪ field fromthis inversion. In comparison with Figure 3B, Figure
4D reveals a striking similarity with the true ␪ field,
indicating the advantage of fusing point measurementsof ␪ and an ERT survey.
Similar to Case 1, 34 moisture contents were sampledfrom the estimated ␪ field. Combining these samples
with the six directly measured ␪ values and one ␣ sample, our hydraulic inversion model produced an estimated ln ␣ field (Fig. 4E). As expected, the estimatedln ␣ field showed better agreement with the true ln ␣field (Fig. 3A) than the one obtained from Case 1. This
is attributed to the fact that more accurate ␪ information,acquired from fusion of information, was included inthis inversion than in Case 1. Figure 4F shows a plot of
the simulated ␪ distribution from the estimated ln ␣field, which is only slightly better than Figure 4D. This
was attributed to greater effect of the other prescribed
unsaturated hydraulic parameters. Nonetheless, if more␪ information from the ERT survey (Fig. 4D) had been
used, the hydraulic inversion would have yielded aneven better estimate of the ln ␣ field and, in turn, abetter simulation of the moisture content distribution
than Figure 4F.

Final Remarks
Under field conditions, stochastic conceptualizationof the inverse problem is inevitable. Geostatistically
based inverse procedures appear promising because
they utilize prior statistical information of the geologicalstructure, point observations, and flow processes to yield
conditional hydraulic parameters and information regarding their uncertainty. To further reduce uncertaintyof the conditional hydraulic parameters, and to makethe parameters more representative of reality, assimilating a broad range of hydrologic, geologic, and geophysical information in the inversion becomes essential. We
believe that the stochastic fusion of information approach presented in this paper offers much promise toaccomplish this task.

ACKNOWLEDGMENTS
The first author is funded in part by a DOE EMSP96 grantthrough Sandia National Laboratories (contract AV-0655#1)
and a DOE EMSP99 grant through the University of Wisconsin, A019493, and in part by an EPA grant R-827114-01-0.
This work does not necessarily reflect the views of DOE and
EPA, and no official endorsement should be inferred. The
study was also partially supported by the SAHRA Scienceand Technology Center as part of NSF grant EAR-9876800.
We are grateful to Alexandre Desbarats, Roger Beckie, Martha P.L. Whitaker and an anonymous reviewer, for their reviews of the manuscript and useful suggestions and comments.

REFERENCES
Baker, K. 2001. Investigation of direct and indirect hydraulic property
laboratory characterization methods for heterogeneous alluvial deposits: Application to the Sandia-Tech vadose zone infiltration testsite. Master’s thesis, Geosciences. New Mexico Tech, Socorro, NM.
Benito, P.H., P. Cook, B. Faybishenko, B. Freifeld, and C. Doughty.
1998. Box Canyon pneumatic connectivity study: Preliminary data
analysis. Report LBNL-42359. Lawrence Berkeley National Laboratory, Berkeley, CA.
Benito, P.H., P. Cook, B. Faybishenko, B. Freifeld, and C. Doughty.
1999. Crosswell air-injection packer tests for the assessment ofpneumatic connectivity in fractured, unsaturated basalt. p. 843–850.
In Amadei et al. (ed.) Rock mechanics for industry.
Binley, A., P. Winship, R. Middleton, M. Pakar, and J. West. 2001.
High-resolution characterization of vadose zone dynamics using
cross-borehole radar. Water Resour. Res. 37:2639–2652.
Brainard, J.R., R.J. Glass, D.L. Alumbaugh, L. Paprocki, D. LaBrecque, X. Yang, T.-C.J. Yeh, K.E. Baker, and C.A. Rautman.
2001. The Sandia-Tech Vadose Zone Facility: Experimental designand data report of a constant flux infiltration experiment. Internal
Report. Sandia National Laboratories, Albuquerque, NM.
Butler, J.J., Jr., C.D. McElwee, and G.C. Bohling. 1999. Pumping
tests in networks of multilevel sampling wells: Methodology andimplications for hydraulic tomography. Water Resour. Res. 35:
3553–3560.
Daily, W., A. Ramirez, D.J. LaBrecque, and J. Nitao. 1992. Electricalresistivity tomography of vadose water movement. Water Resour.
Res. 28:1429–1442.
Dettinger, M.D., and J.L. Wilson. 1981. First-order analysis of uncertainty in numerical models of groundwater flow. Part 1. Mathematical development. Water Resour. Res. 17:149–161.
Gardner, W.R. 1958. Some steady state solutions of the unsaturated
moisture flow equation with application to evaporation from awater table. Soil Sci. 85:228–232.
Gardner, W.R., and F.J. Miklich. 1962. Unsaturated conductivity anddiffusivity measurements by a constant flux method. Soil Sci.
93:271–274.
Gelhar, L.W. 1993. Stochastic subsurface hydrology. Prentice Hall,Englewood Cliffs, NJ.
Gottlieb, J., and P. Dietrich. 1995. Identification of the permeability
distribution in soil by hydraulic tomography. Inverse Problems
11:353–360.
Gutjahr, A.L., and J.L. Wilson. 1989. Cokriging for stochastic models.
Transp. Porous Media 4:585–598.
Guzman, A.G., A.M. Geddis, M.J. Henrich, C.F. Lohrstorfer, andS.P. Neuman. 1996. Summary of air permeability data from singlehole injection tests in unsaturated fractured tuffs at the Apache
Leap Research Site: Results of steady-state test interpretation.
NUREG/CR-6360. U.S. Nuclear Regulatory Commission, Washington, DC.
Harter, Th., and T.-C.J. Yeh. 1996. Conditional stochastic analysis ofsolute transport in heterogeneous, variably saturated soils. Water
Resour. Res. 32:1597–1610.
Harvey, C., and S.M. Gorelick. 1995. Mapping hydraulic conductivity:
Sequential conditioning with measurements of solute arrival time,
hydraulic head, and local conductivity. Water Resour. Res. 31:
1615–1626.
Hillel, D., V.D. Krentos, and Y. Stylianou. 1972. Procedure and test

220

VADOSE ZONE J., VOL. 1, NOVEMBER 2002of an internal drainage method for measuring soil hydraulic characteristics in situ. Soil Sci. 114:395–400.
Hoeksema, R.J., and P.K. Kitanidis. 1984. An application of the geostatistical approach to the inverse problem in two-dimensional
groundwater modeling. Water Resour. Res. 20:1003–1020.
Hoeksema, R.J., and P.K. Kitanidis. 1989. Prediction of transmissivities, heads, and seepage velocities using mathematical modelingand geostatistics. Adv. Water Resour. Res. 12:90–101.
Huang, K., Y.W. Tsang, and G.S. Bodvarsson. 1999. Simultaneous
inversion of air-injection tests in fractured unsaturated tuff at Yucca
Mountain. Water Resour. Res. 35:2375–2386.
Hubbard, S. J. Chen, J. Peterson, E.L. Majer, K.H. Williams, D.J.
Swift, B. Mailloux, and Y. Rubin. 2001. Hydrological characterization of the South Oyster bacterial transport site using geophysicaldata. Water Resour. Res. 37:2431–2456.
Hubbard, S. Y. Rubin, and E. Majer. 1999. Spatial correlation structure estimation using geophysical data. Water Resour. Res. 35:
1809–1825.
Hudson, D.B., P.J. Wierenga, and R.G. Hills. 1996. Unsaturated hydraulic properties from upward flow into soil cores. Soil Sci. Soc.
Am. J. 60:388–396.
Hughson, D.L., and T.-C. J. Yeh. 1998. A geostatistically-based inverse
model for three-dimensional variably saturated flow. Stoch. Hydrol. Hydraul. 12:285–298.
Hughson, D.L., and T.-C.J. Yeh. 2000. An inverse model for threedimensional flow in variably saturated porous media. Water Resour. Res. 36:829–839.
Hyndman, D.W., and S.M. Gorelick. 1996. Estimating lithologic andtransport properties in three-dimensions using seismic and tracer
data. Water Resour. Res. 32:2659–2670.
Illman, W.A., and S.P. Neuman. 2001. Type-curve interpretation ofa cross-hole pneumatic injection test in unsaturated fractured tuff.
Water Resour. Res. 37:583–603.
Illman, W.A., D.L. Thompson, V.V. Vesselinov, and S.P. Neuman.
1998. Single-hole and cross-hole pneumatic tests in unsaturated
fractured tuffs at the Apache Leap Research Site: Phenomenology,
spatial variability, connectivity and scale. Rep. NUREG/CR-5559.
U.S. Nuclear Regulatory Commission, Washington, DC.
Kitanidis, P.K. 1986. Parameter uncertainty in estimation of spatial
functions: Bayesian analysis. Water Resour. Res. 22:499–507.
Kitanidis, P.K. 1997. Comment on “A reassessment of the groundwater inverse problem” by D. McLaughlin and L.R. Townley. Water Resour. Res. 33:2199–2202.
Kitanidis, P.K., and E.G. Vomvoris. 1983. A geostatistical approach tothe inverse problem in groundwater modeling and one-dimensional
simulations. Water Resour. Res. 19:677–690.
Kool, J.B., J.C. Parker, and M.Th. van Genuchten. 1985. Determining
soil hydraulic properties for one-step outflow experiments by parameter estimation. I. Theory and numerical studies. Soil Sci. Soc.
Am. J. 49:1348–1354.
LeCain, G.D. 1996. Air-injection testing in vertical boreholes inwelded and nonwelded tuff, Yucca Mountain, Nevada. U.S. Geol.
Survey Water-Resources Investigations Report 96–4262. U.S. Geol.
Survey, Reston, VA.
LeCain, G.D., Results from air-injection and tracer testing in theUpper Tiva Canyon, Bow Ridge Fault, and Upper Paintbrush Contact Alcoves of the Exploratory Studies Facility, August 1994
through July 1996, Yucca Mountain, Nevada. U.S. Geol. Survey
Water-Resources Investigations Report 98–4058. U.S. Geol. Survey, Reston, VA.
LeCain, G.D., and J.N. Walker. 1994. Results of air-permeability
testing in a vertical borehole at Yucca Mountain, Nevada. p. 2782–
2788. Proceedings of the Fifth Annual International Conference
on High-Level Radioactive Waste Management. 22–26 May 1994.
Las Vegas, NV. ASCE, Reston, VA.
Li, B., and T.-C.J. Yeh. 1999. Cokriging estimation of the conductivity
field under variably saturated flow conditions. Water Resour. Res.
35:3663–3674.
Liu, S. 2001. A sequential inverse approach for hydraulic tomography and electrical resistivity tomography: An effective method
for site characterization. Ph.D. diss., The University of Arizona,
Tucson, AZ.
Liu, S., T.-C.J. Yeh, and R. Gardiner. 2002. Effectiveness of hydrau-

lic tomography: sandbox experiments. Water Resour. Res. 38(4):10.1029/2001WR000338.
McLaughlin, D., and L.R. Townley. 1996. A reassessment of thegroundwater inverse problem. Water Resour. Res. 32:1131–1161.
Menke, W. Geophysical data analysis: Discrete inverse theory. Vol.
45. International Geophysics Series. Acad. Press, San Diego, CA.
Mulla, D.J. 1988. Estimating spatial patterns in water content, matricsuction, and hydraulic conductivity. Soil Sci. Soc. Am. J. 52:
1547–1553.
Neuman, S.P. 1973. Calibration of distributed parameter groundwaterflow models viewed as a multiple-objective decision process under
uncertainty. Water Resour. Res. 9:1006–1021.
Neuman, S.P., G.R. Walter, H.W. Bentley, J.J. Ward, and D.D. Gonzalez. 1984. Determination of horizontal aquifer anisotropy with threewells. Ground Water 22:66–72.
Rasmussen, T.C., D.D. Evans, P.J. Sheets, and J.H. Blanford. 1990.
Unsaturated fractured rock characterization methods and data sets
at the Apache Leap Tuff Site. NUREG/CR-5596. U.S. Nuclear
Regulatory Commission, Washington DC.
Rasmussen, T.C., D.D. Evans, P.J. Sheets, and J.H. Blanford. 1993.
Permeability of Apache Leap tuff: Borehole and core measurements using water and air. Water Resour. Res. 29:1997–2006.
Rea, J., and R. Knight. 1998. Geostatistical analysis of ground-penetrating radar data: A means of describing spatial variation in thesubsurface. Water Resour. Res. 34:329–339.
Rose, C.W., W.R, Stern, and J.E. Drummond. 1965. Determinationof hydraulic conductivity as a function of depth and water content
for soil in situ. Aust. J. Soil Res. 3:1–9.
Rubin, Y., and G. Dagan. 1987. Stochastic identification of transmissivity and effective recharge in steady groundwater flow. 1. Theory.
Water Resour. Res. 23:1185–1192.
Rubin, Y., G. Mavko, and J. Harris. 1992. Mapping permeability inheterogeneous aquifers using hydrological and seismic data. Water
Resour. Res. 28:1809–1816.
Šimůnek, J., and M.Th. van Genuchten. 1997. Estimating unsaturated
soil hydraulic properties from multiple tension disc infiltrometerdata. Soil Sci. 162:383–398.
Šimůnek, J., O. Wendroth, and M.Th. van Genuchten. 1998. A parameter estimation analysis of the evaporation method for determining
soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.
Srivastava, R., and T.-C.J. Yeh. 1992. A three dimensional numericalmodel for water flow and transport of chemically reactive solute
through porous media under variably saturated conditions. Adv.
Water Resour. Res. 15:275–287.
Sun, N.Z. 1994. Inverse problems in groundwater modeling. Kluwer
Acad., Norwell, MA.
Sun, N.-Z., and W.W.-G. Yeh. 1992. A stochastic inverse solution fortransient groundwater flow: Parameter identification and reliability
analysis. Water Resour. Res. 28:3269–3280.
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.
Trautz, R.C. 1984. Rock fracture aperture and gas conductivity measurements in-situ. M.S. thesis. Dep. of Hydrol. and Water Resour.
Univ. of Ariz., Tucson, AZ.
van Dam, J.C., J.N.M. Stricker, and P. Droogers. 1994. Inverse methodfor determining soil hydraulic functions 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.
Vargas-Guzman, A.J. and T.-C.J. Yeh. 1999. Sequential kriging andcokriging: Two powerful geostatistical approaches. Stoch. Environ.
Res. Risk Assess. 13:416–435.
Vauclin, M., S.R. Vieira, G. Vachaud, and D.R. Nielsen. 1983. The
use of cokriging with limited field soil observations. Soil Sci. Soc.
Am. J. 47:175–184.
Vesselinov, V.V., S.P. Neuman, and W.A. Illman. 2001a. Threedimensional numerical inversion of pneumatic cross-hole tests inunsaturated fractured tuff. 1. Methodology and borehole effects.
Water Resour. Res. 37:3001–3017.
Vesselinov, V.V., S.P. Neuman, and W.A. Illman. 2001b. Threedimensional numerical inversion of pneumatic cross-hole tests inunsaturated fractured tuff. 1. Equivalent parameters, high-resolu-

www.vadosezonejournal.org

tion stochastic imaging and scale effects. Water Resour. Res. 37:
3019–3041.
Wang, J.S.Y., P.J. Cook, R.C. Trautz, R. Salve, A.L. James, S. Finsterle, T.K. Tokunaga, R. Solbau, J. Clyde, A.L. Flint, and L.E. Flint.
1998. Field testing and observation of flow paths. In Niches: Phase
1 Status Report of the Drift Seepage Test and Niche Moisture
Study, Level 4 milestone SPC314M4 for WBS 1.2.3.3.1.2.4. Yucca
Mountain Site Characterization Project.
Watson, K.K. 1966. An instantaneous profile method for determining
the hydraulic conductivity unsaturated porous materials. Water
Resour. Res. 2:709–715.
Yates, S.R., and A.W. Warrick. 1987. Estimating soil water content
using cokriging. Soil Sci. Soc. Am. J. 51:23–30.
Yeh, T.-C.J., A.L. Gutjahr, and M. Jin. 1995. An iterative cokriginglike technique for groundwater flow modeling. Ground Water 33:
33–41.
Yeh, T.-C., M. Jin, and S. Hanna. 1996. An iterative stochastic inverse
method: Conditional effective transmissivity and hydraulic head
fields. Water Resour. Res. 32:85–92.
Yeh, T.-C.J., and S. Liu. 2000. Hydraulic tomography: Developmentof a new aquifer test method. Water Resour. Res. 36:2095–2105.
Yeh, T.-C.J., S. Liu, R.J. Glass, K. Baker, J.R. Brainard, D. Alumbaugh, and D. LaBrecque. 2002. A geostastically-based inverse

221

model for electrical resistivity surveys and its applications to vadosezone hydrology. Water Resour. Res. In press.
Yeh, T.C., T.C. Rasmussen, and D.D. Evans. 1988. Simulation ofliquid and vapor movement in unsaturated fractured rock at the
Apache Leap Tuff Site—Models and strategies. NUREG/CR-5097.
U.S. Nuclear Regulatory Commission, Washington DC.
Yeh, T.-C.J., and J. Zhang. 1996. A geostatistical inverse method forvariably saturated flow in the vadose zone. Water Resour. Res. 32:
2757–2766.
Yeh, W.W.-G. 1986. Review of parameter identification proceduresin groundwater hydrology: The inverse problem. Water Resour.
Res. 22:95–108.
Zhang, J., and T.-C.J. Yeh. 1997. An iterative geostatistical inverse
method for steady flow in the vadose zone. Water Resour. Res.
33:63–71.
Zhou, Q.Y., J. Shimada, and A. Sato. 2001. Three-dimensional spatial
and temporal monitoring of soil water content using electrical resistivity tomography. Water Resour. Res. 37:273–285.
Zurmühl, T. 1996. Evaluation of different boundary conditions forindependent determination of hydraulic parameters using outflow
methods. p. 81–84. In W. Durner et al. (ed.) Advanced methods
to determine hydraulic properties of soils. Dep. of Hydrology University of Bayreuth, Germany.