HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Yeh, T.-C. J. Articles by Simunek, J. Search for Related Content PubMed Articles by Yeh, T.-C. J. Articles by Simunek, J. GeoRef GeoRef Citation Agricola Articles by Yeh, T.-C. J. Articles by Simunek, J. Related Collections Spatial Distribution Vadose Zone Processes and Chemical Transport Soil Hydrology

Reviews and Analyses

Stochastic Fusion of Information for Characterizing and Monitoring the Vadose Zone

Tian-Chyi J. Yeh^*,a and Jirka imnek^b

^a Department of Hydrology and Water Resources, The University of Arizona, Tucson, AZ 85721
^b George E. Brown, Salinity Laboratory, USDA-ARS, Riverside, CA 92507-4617
* Corresponding author (yeh{at}hwr.arizona.edu)

Received 12 March 2002.

	ABSTRACT

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
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 for the inverse problems to be well posed and discuss difficulties associated with solving the inverse problems. We subsequently explain the need for a stochastic conceptualization of inverse problems of the vadose 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 classical vadose 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 tool for effectively characterizing heterogeneity, monitoring processes in the vadose zone, and quantifying uncertainties associated with vadose zone characterization and monitoring.

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])

	INTRODUCTION

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
KNOWLEDGE OF UNSATURATED SOIL hydraulic properties is indispensable for successful predictions of water 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 (Kool et al., 1985; van Dam et al., 1994), the upward infiltration method (Hudson et al., 1996), and the evaporation method (Gardner and Miklich, 1962, imnek et al., 1998). Popular field methods include the instantaneous profile method (Hillel et al., 1972), various unit-gradient type 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 a number 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 the above methods are widely used and accepted, questions 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 to about 500 cm³. The scale of field methods generally does not extend beyond a plot of 1 m² and depths of one to several meters. There is an urgent need to develop methods that characterize hydraulic properties of the vadose zone on a much larger scale. Recently developed geophysical methods such as electrical resistivity tomography (ERT) (Daily et al., 1992; Zhou et al., 2001) and ground-penetrating radar (GPR) (Binley et al., 2001), coupled with geostatistically based inverse methods (Liu, 2001; Yeh et al., 2002) appear to be promising tools for monitoring large-scale vadose zones. Combining these new tools with a three-dimensional hydrologic inverse model (e.g., Hughson and Yeh, 2000), a detailed three-dimensional characterization of the vadose zone on a large scale appears to become possible.

In this review we examine conditions required to make inverse problems well posed for unsaturated flow through homogeneous and heterogeneous soil columns. Stochastic conceptualizations of inverse problems for the vadose zone are then introduced. Simple examples are used to illustrate the principles of geostatistically based inverse approaches (GA), including cokriging, a sequential 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 development of integrative inversion approaches for characterizing and monitoring the vadose zone.

	The Inverse Problem in Subsurface Hydrology: An Ill- or Well-Posed Problem

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
Movement of water in variably saturated porous media, assuming isothermal conditions and neglecting the gas phase, may be described using a modified form of the Richards equation

[1]

where h is the pressure head, which is positive when the medium is fully saturated and negative when the medium is partially saturated, t is time, represents the volumetric moisture content, and z is the positive upward vertical coordinate. The term S_s represents specific storage and a transitioning parameter that is one when h 0 and zero when h is negative. To describe the –h relationship of unsaturated media, van Genuchten's (1980) model is often assumed:

[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 with m = 1 - 1/n. We further assume that the unsaturated hydraulic conductivity function, K(h), follows Mualem's pore-size distribution model (van Genuchten, 1980); that is,

[3]

where K_s 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, ₁, where n is a normal unit vector and q* a specified specific discharge, and (ii) (x, y, z) = * at boundaries, ₂, where * is the specified hydraulic head. The initial condition is given as (x, y, z, t) = *(x, y, z, 0) at t = 0.

If we define K_s, S_s, , 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 with known primary variables, and for given initial and boundary conditions. On the other hand, an inverse problem refers to estimating values of the primary variables from information 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 can have 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, a well-posed problem with incorrect information can lead to 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 steady-state and transient unsaturated flow in porous media.

Steady-State Unsaturated Flow
The fact that the hydraulic conductivity is a function of the pressure head or moisture content generally increases the difficulties of inverse problems for flow through unsaturated porous media. Consider steady one-dimensional vertical flow in an unsaturated soil column and assume Darcy's Law to be valid for the flow process:

[4]

In this case, the pressure head, h, along the column must be specified in addition to the specific discharge so that the unsaturated hydraulic conductivity at the given pressure head can be determined. Based on this principle, a unit gradient approach is frequently used to determine the unsaturated hydraulic conductivity. Since hydraulic conductivity 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) = K_s exp(ßh), where K_s was previously defined and ß is a pore-size distribution parameter (Gardner, 1958), we have

[5]

In this case, two independent equations are required to uniquely determine the value of K_s and ß, implying that at 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 (K_s, , 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 of each 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:

[6]

Each of these equations has two unknowns (K_s and ß) if q and h are specified. For unsaturated flow, h varies nonlinearly within each material. Accordingly, a detailed h distribution for one steady flow rate is necessary to allow specification of h and dh/dz values at different z values, such that two independent equations can be formulated for each material. If only one h and dh/dz (e.g., unit gradient conditions if they exist) in each material is known, imposing two different flow rates is essential so that two independent equations for each material exist. To demonstrate this concept, Fig. 1A and 1B show both the true and estimated values of K_s and (after removing their means) in a synthetic soil profile with 20 layers, respectively. The estimated values were obtained using the successive linear estimator approach (discussed below) with one pressure head measurement in each layer for one steady unsaturated flux situation. This inverse problem is ill posed since each layer has two unknowns, K_s and , but only one pressure head value, 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 different pressure 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 K_s and values with the two fluxes. Notice that the estimates are identical to the true values—the inverse problem is well posed. The parameter n in the VG model was assumed to be known a priori in this example.

View larger version (27K): [in this window] [in a new window]   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.

[7]

The term on the left-hand side of Eq. [7] is the outflow evaluated at z = 0. The first term on the right-hand side of Eq. [7] denotes changes in storage, while the second term is the specific discharge flowing into the column at z = L. The sum of the two terms on the right-hand side of Eq. [8] yields the specific discharge at the end of the column, 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 number of 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 (K_s, , and n) is used, at least three different discharge rates and pressure head profiles at three times are needed to create three independent equations. This scenario forms the basis of the well-known 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 and the change in moisture content along the length of the column. However, this is not sufficient to make the inverse problem of Eq. [8] well posed. Since the unsaturated hydraulic conductivities K(h) at z = 0 and z = L remain unknown, the number of unknowns is greater than 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. Notice that 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 the inverse problem well posed.

These examples assume that the size of each soil block is 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, the same principles are also applicable to multidimensional flow problems. The hydraulic conductivity anisotropy of the equivalent homogeneous medium, or the anisotropy of each block, may have to be considered if they are significant. This means that more parameters may need to be identified. Therefore, additional information about the system's spatial and temporal responses must be acquired.

The discussions above lead to the fact that if sufficient information is available to yield enough independent equations, the inverse problem is always well posed and should 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 mathematical models, which assume a one-to-one relationship between the hydraulic property and the hydraulic response of a system in both forward and inverse operators.

The inversion of hydraulic properties from actual laboratory and field experiments may encounter other complications. These complications include representativeness of mathematical models for actual physical processes (model error), inability to measure detailed h and 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; Zurmühl, 1996; imnek and van Genuchten, 1997; imnek et al., 1998) have developed 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 Conceptualization of Inverse Problems

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
While inverse problems associated with laboratory experiments in general can be well posed if sufficient and necessary information is collected, virtually all inverse problems for field experiments are ill posed because the prerequisites 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 is well 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 to sparsely observed responses of a system, is given, the solution 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 hydraulic properties becomes appropriate. That is to say, each of the properties of a geological formation should be considered 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 be approximated by the first and second moments of samples. The first moment, the mean, represents the most likely value of the property. The second moment, the spatial covariance function, specifies the variance and correlation structure of the process, analogous to a description of maximum and minimum values of the properties and the average block size in previous discussions, respectively. If one adopts the stochastic representation of 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 conditional mean field is unique. This field also represents the most likely solution to the inverse problem, even though this may 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 gradually resemble 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

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
If one accepts the stochastic conceptualization of the ill-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 that considers spatial correlation structures of flow processes (such as pressure heads and velocities) and hydraulic properties 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 in highly heterogeneous aquifers (e.g., Kitanidis and Vomvoris, 1983; Hoeksema and Kitanidis, 1984, 1989; Rubin and Dagan, 1987; Gutjahr and Wilson, 1989; Sun and Yeh, 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 water content, 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₀, at location x₀, using simple cokriging and a known hydraulic conductivity, f₁, and pressure head, h₂, at locations x₁ and x₂, respectively; that is,

[9]

where the weights ₀₁ and ß₀₂, represent the contribution to the hydraulic conductivity estimate at x₀ from the known hydraulic conductivity, f₁, and the pressure head, h₂, at locations x₁ and x₂, respectively. The weights can be obtained by minimizing the mean-square error of the estimate, E[(f₀ - ₀)²]. The minimization leads to a system of equations:

[10]

where C_gg (x_i, x_j) represents the spatial covariance of the stochastic processes, g_i and g_j (which can be either conductivity or head) at locations x_i and x_j. Notice that C_fh(x_i, x_j) = C_hf(x_j, x_i). The solution of Eq. [10] yields the cokriging weights,

[11]

if we define normalized covariance functions, that is, _ff(x_i, x_j) = C_ff(x_i, x_j)/C_ff(x_i, x_i), _hh(x_i, x_j) = C_hh(x_i, x_j)/C_ff(x_i, x_i), and _fh(x_i, x_j) = C_fh(x_i, x_j)/C_ff(x_i, x_i).

Note from Eq. [11] that the weight does not depend on the variance of hydraulic conductivity. To evaluate the weights, the covariance, C_ff(x_i, x_j) must be specified a priori. In theory, the covariances C_hh(x_i, x_j) and C_fh(x_i, x_j) can be derived from data, if sufficient data sets are available, but are generally calculated using C_ff(x_i, x_j) 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 = a₂₁f₁ + a₂₂f₂, where a₂₁ and a₂₂ are coefficients. We are to estimate the hydraulic conductivity value, f₂, at location x₂, using simple cokriging with known hydraulic conductivity, f₁, and hydraulic head, h₂, at location x₁ and x₂. In this case, cokriging yields an exact solution. The covariances needed are calculated from the forward model and they are

[12]

Using these covariances, a solution to a matrix similar to Eq. [10] lead to the weights: ₂₁ = -a₂₁/a₂₂, and ß₂₂ = 1/a₂₂. Notice that in this case the weights are independent of the covariance function of hydraulic conductivity. 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 the covariance functions (see Eq. [11]) and can produce a unique estimate that satisfies the minimum mean square error criterion. The influence of the covariance, however, decreases with an increase in the number of observations of either h or f.

Sequential Linear Estimator
Instead of estimating f₀ simultaneously using both f₁ and h₂ as in the previous example, we now will use a sequential linear estimator approach (Vargas-Guzman and Yeh, 1999) that estimates f₀, f₁, and h₂ using only f₁ at first, and then improves the estimate by adding the observed h₂ information. Our first estimates are

[13]

where the superscript (1) indicates the first step of the sequential estimate and

[14]

Notice that the estimate of f₁ 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₀ by using the difference between the observed h₂ and its estimate from Eq. [13]:

[15]

where the superscripted (1) and (2) indicate the order of the sequential linear estimator. The weight, ₀₂, in Eq. [15] is obtained by minimizing the variance of the new estimate:

[16]

where represents the conditional (or residual) covariance and cross-covariance. Differentiating Eq. [16] with respect to ₀₂ and setting the result equal to zero, we have

[17]

in which the conditional covariance can be determined by

[18]

Using Eq. [14], the normalized conditional head covariance for Eq. [18] becomes

[19]

Equation [19] shows that the conditional variance of the hydraulic head thus becomes smaller than the unconditional head variance because of the effects of conditioning using f₁. Similarly, the conditional cross-covariance can be derived:

[20]

Written in terms of normalized covariances, Eq. [20] becomes

[21]

Then the weight is given by

[22]

The estimate (Eq. [15]) thus becomes

[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 that uses all of the information simultaneously. The sequential approach nevertheless reduces computational efforts. Consider there are N_f measurements of hydraulic conductivity and N_h measurements of the hydraulic head. The computational cost of inverting the matrix in Eq. [10] is (N_f + N_h)³, whereas the cost of the sequential estimator is N³_f + N³_h. 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 be emphasized 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 ₁₁ and ß₁₂ for Location 1, independently from solving the equation for obtaining ₂₁ and ß₂₂ for Location 2. Therefore, the algorithm is most suitable for parallel computing, and is highly efficient, since the inverse of the left-hand side of Eq. [11] has to be computed only once. This feature of the GA approach offers a similar advantage to inverse modeling of variably saturated flows where several parameters (e.g., K_s, , n, _s, and _r of the VG model) require estimation—the GA can estimate these parameters independently.

In the case of stochastic inverse problems where the number of observations, M, is always much smaller than the number of parameters to be estimated, N, which is the scenario in real-world problems, the GA method offers an even greater advantage over other classical inverse methods. For inverse problems of variably saturated flow, the GA also permits sequential inclusion of pressure head and moisture content measurements. Moreover, the GA method allows sequential inclusion of observations from multiple experiments such as in hydraulic 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 and secondary variables of the vadose zone are generally highly nonlinear. Therefore, the GA approach cannot fully exploit available secondary information. To overcome this limitation, Zhang and Yeh (1997) adapted a successive linear estimator (SLE) technique (Yeh et al., 1996) to the vadose zone inverse problem. The principle of 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 is available. 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 field from Step 1 is used to simulate the hydraulic head by means of a forward flow model. In Step 4, the conditional covariances of the secondary information and their 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 to compute new weights. Step 6 is where the new weights, along with the difference between simulated and observed system responses at observation locations, are used 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 probability of the primary variable in a consistent manner using a first-order approximation approach. For a well-posed inverse problem (saturated or unsaturated flow), the SLE converges to the true solution rapidly (as shown in Fig. 1). Also, notice that the SLE can account for measurement errors in the estimation by addition of error variances to the diagonal of the left-hand side of matrix of (10) or (17).

	Applications

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
Unsaturated Inverse Problems
Traditional inverse modeling approaches typically consider 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 (1996) to identify unsaturated hydraulic parameters in heterogeneous 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 head enhance the prediction of the flow paths of solutes in the vadose zone. Yeh and Zhang (1996) reported that the unsaturated hydraulic parameters of a heterogeneous vadose zone could be identified if sufficient pressure head and degree of saturation information was available. 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 a hypothetical, heterogeneous vadose zone subject to steady-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, and to the linear predictor nature of cokriging.

Among pressure head, solute concentration, and solute arrival time measurements, Li and Yeh (1999) found that pressure head measurements of steady-state flow fields were the most useful secondary information for estimating the K_s field using cokriging. They attributed this finding to several factors. First, the nonlinear relationship between head and K_s may have been relatively mild for the cases they studied. Additionally, the assumption of ergodicity was approximately satisfied for steady-state flow. In other words, the ensemble mean head distribution evaluated with mean values of the parameters closely approximates the spatial mean head distribution. Such an existence of ergodicity reduces the variance of the head, and consequently improves the linearity between head and conductivity. Conversely, the ergodicity assumption cannot be easily satisfied for the solute transport case. Only when a solute plume has traveled over enough correlation lengths, has sampled enough heterogeneity, and has become a Fickian process, will the ergodicity assumption be satisfied. Because of 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 is a function of not only the hydraulic conductivity but also the hydraulic gradient and the moisture content. Cross-correlation between the concentration and the hydraulic conductivity is thus not as strong as the relationship between pressure head and hydraulic conductivity. This suggests that for estimating the hydraulic conductivity, concentration measurements are generally less effective. The same argument applies to arrival times of a tracer concentration. However, Li and Yeh (1999) found that estimates of the conductivity based on head measurements could still be improved by incorporating additional concentration information.

Of course, the examples assumed that the scales of measurements 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, a practical solution is to collect "enough" measurements in space (e.g., along a borehole), smooth them with a moving 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 for the measurements in the SLE such that contributions from the measurements to estimates are reduced. As a consequence, estimated fields are expected to be smooth, but they may retain essential heterogeneity patterns. Lastly, one must recognize that the propagation process of pressure head is highly diffusive, moisture content is less, 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 vadose zone inverse problems, Zhang and Yeh (1997) showed that the SLE approach could yield more detailed images of unsaturated hydraulic parameters than cokriging. This implies that the SLE, which successively incorporates the nonlinear relation between the primary and secondary 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 to allow pressure head and water content data obtained at different times to be sequentially included in the inversion. Using this inverse model, they investigated the 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 vadose zone 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 operations 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; Yeh and 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, and monitoring the steady-state head response at the other wells 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 the system 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 hydraulic tomography 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 approach and concluded that the tomography is a viable technology for delineate aquifer heterogeneity and the sequential SLE algorithm is a promising tool for inverting the vast amount of data generated with tomography tests.

While their investigations dealt with saturated flow problems, the tomographic concept can be equally well applied to unsaturated infiltration tests, that is, unsaturated hydraulic tomography. The unsaturated hydraulic tomography method involves the sequential injection of water at different depths and locations, followed by monitoring the response of the vadose zone at selected locations. Using an appropriate inverse modeling technique, the vast amount of secondary information can thus 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 the permeability and porosity of unsaturated soil and fractured rocks. In principle, air injection tests are very similar to their hydraulic counterparts conducted in fully water-saturated media. Air is either injected into or withdrawn from sections of boreholes isolated by means of 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 and porosity through analytical techniques or numerical inverse modeling. Pneumatic tomography (i.e., cross-hole pneumatic injection tests conducted sequentially and corresponding 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 and Walker, 1994; LeCain, 1996, 1998; Wang et al., 1998; Huang et al., 1999), in Box Canyon, Idaho (Benito et al., 1998, 1999), and at the Apache Leap Research Site near 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 et al., 1998; Illman and Neuman, 2001), which cause the effectiveness 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.

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 be practical 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 electrical potential at other locations. By moving the DC source location, one can generate many electrical potential fields and source pairs, from which a three-dimensional image of changes in the resistivity can be derived using an inverse model. Assuming changes in resistivity mimic changes in moisture content, water flow in the vadose 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 geological material 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 pattern of change in the resistivity and, thus, the pattern of change in the moisture content. Furthermore, the resolution of the ERT image can be easily improved by deploying a large number of sensors within boreholes and at the soil surface. The relative inexpensiveness of the sensor 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 of Archie's Law exhibit profound spatial variability, although they are spatially correlated. This variability compounds the difficulty in translating resistivity to moisture content. In other words, at some location, a small change in resistivity may indicate a large change in water content, while at other locations a substantial change in the resistivity may correspond to only a small change in the water content (Yeh et al., 2002). Consequently, the ERT's ability to yield an accurate image of change in moisture content, or more critically the moisture content itself, remains to be proven.

	Stochastic Fusion of Information

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
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 volume of 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 collect a vast amount of information pertinent to hydrological inversion; however, significant spatial variation in the relations between various hydrologic and geophysical properties diminishes the value of the information for its intended applications.
   

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 and cons, an integration of their individual strengths may facilitate a better way to monitor and characterize the vadose zone and quantify uncertainty.

During the integration, not only the hydrological and geophysical 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 information and 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 all available information. Consequently, conditional means, the best and unique estimates with minimal uncertainty, can be obtained and their conditional variances can quantify 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 to unsaturated hydraulic tomography or ERT surveys to independently 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 hydraulic properties 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 measurements of the electric field induced by the ERT survey, point measurements of resistivity and moisture content, and parameters of resistivity–moisture relations. In addition, the prior information about covariance structures of resistivity, 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 and their associated covariance functions. For the ERT, conditional mean resistivity and moisture content fields and mean parameter fields of the resistivity–moisture content relation, in addition to their conditional covariance structures, are products of the Level 1 fusion.

Level 2 fusion aims to honor the reciprocal nature of hydraulic and ERT inversion. It takes an iterative approach to arrive at the best estimate of the primary variable fields and uncertainty. The flow chart shown in Fig. 2 depicts a general concept of the fusion process. Specifically, during an unsaturated hydraulic tomography experiment, several water infiltration tests at different locations (packed-off intervals) could be conducted sequentially. 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 improved by using Level l fusion to integrate measured hydrologic and geophysical information into the ERT inversion process along with some prior information. The prior information, the mean and covariance structure of the resistivity and moisture content fields during the unsaturated hydraulic tomography test, can be derived from a forward simulation of the hydraulic tomography test with mean hydraulic parameters. With this information, the ERT inversion can yield a reasonable result and is able to estimate not only changes in resistivity and the water content, but also the moisture content itself and their conditional moments.

View larger version (28K): [in this window] [in a new window]   Fig. 2. Flow chart of Level 2 iterative stochastic fusion.

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 and cost-effective technology for characterizing the vadose zone. Meanwhile, an ERT may become both a reliable and cost-effective monitoring tool for moisture contents.

Example
To demonstrate the promise of the stochastic fusion technology, some preliminary results are presented below. During the past few years, we have developed a method that fuses information obtained from ERT surveys, point hydrological measurements, and unsaturated hydraulic inversion modeling to better characterize and monitor 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 properties of 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, K_s, , and n, were considered as random fields with geometric means of 0.0063 cm s^-1, 0.028 1/cm, and 2.0, respectively. The variances of lnK_s, ln , and ln n were 0.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 scale of 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 center of the vadose zone (from x = 80 to 120 cm, y = 0 to 20 cm, and z = 200 cm) was treated as a constant head boundary 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 was assumed 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 shown in Fig. 3B.

View larger version (31K): [in this window] [in a new window]   Fig. 3. True (A) ln and (B) fields of the synthetic vadose zone.

View larger version (38K): [in this window] [in a new window]   Fig. 4. (A) Estimated fields using six measurements and a kriging technique; (B) estimated ln field from an unsaturated inverse model using 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 the using the ln field from Fig. 4E.

Similar to Case 1, 34 moisture contents were sampled from 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 estimated ln 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 in this 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 an even better estimate of the ln field and, in turn, a better simulation of the moisture content distribution than Figure 4F.

	Final Remarks

TOP ABSTRACT INTRODUCTION The Inverse Problem in... Stochastic Conceptualization of... Geostatistically Based Inverse... Applications Stochastic Fusion of Information Final Remarks REFERENCES

 
Under field conditions, stochastic conceptualization of the inverse problem is inevitable. Geostatistically based inverse procedures appear promising because they utilize prior statistical information of the geological structure, point observations, and flow processes to yield conditional hydraulic parameters and information regarding their uncertainty. To further reduce uncertainty of the conditional hydraulic parameters, and to make the 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 to accomplish this task.

	ACKNOWLEDGMENTS

 
The first author is funded in part by a DOE EMSP96 grant through 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 Science and 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

TOP ABSTRACT INTRODUCTION <