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 Similar articles in ISI Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (12) Citing Articles via Google Scholar Google Scholar Articles by Liu, S. Articles by Yeh, T.-C. J. Search for Related Content PubMed Articles by Liu, S. Articles by Yeh, T.-C. J. GeoRef GeoRef Citation Agricola Articles by Liu, S. Articles by Yeh, T.-C. J. Related Collections Water Content Tomography

ORIGINAL RESEARCH

An Integrative Approach for Monitoring Water Movement in the Vadose Zone

^a Burgess & Niple, 5025 E. Washington St., Phoenix, AZ 85034
^b Dep. of Hydrology and Water Resour., John W. Harshbarger Bldg., The Univ. of Arizona, Tucson, AZ 85721
^* Corresponding author (yeh{at}hwr.arizona.edu).

Received 18 July 2003.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

 
Electrical resistivity tomography (ERT), during the past few years, has emerged as a potentially cost-effective, noninvasive tool for imaging changes of moisture content in the vadose zone. The accuracy of ERT surveys, however, has been the subject of debate because of its nonunique inverse solution and spatial variability in the constitutive relation between resistivity and moisture content. In this paper, an integrative inverse approach for ERT, based on a stochastic information fusion concept of Yeh and imnek, was developed to derive the best unbiased estimate of the moisture content distribution. Unlike classical ERT inversion approaches, this new approach assimilates prior information about the geological and moisture content structures in a given geological medium, as well as sparse point measurements of the moisture content, electrical resistivity, and electric potential. Using these types of data and considering the spatial variability of the resistivity–moisture content relation, the new approach directly estimates three-dimensional moisture content distributions instead of simply changes in moisture content in the vadose zone. Numerical experiments were conducted to investigate the effect of uncertainties in the prior information on the estimate. The effects of spatial variability in the constitutive relation were then examined on the interpretation of the change in moisture content, based on the change in electrical resistivity from the ERT survey. Finally, the ability of the integrative approach was tested by directly estimating moisture distributions in three-dimensional, heterogeneous vadose zones. Results show that the integrative approach can produce accurate estimates of the moisture content distributions and that incorporating some measurements of the moisture content is essential to improve the estimate.

Abbreviations: ERT, electrical resistivity tomography

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

 
ELECTRICAL RESISTIVITY SURVEYS are increasingly used to collect extensive electric current and potential data in multiple dimensions to image subsurface electrical resistivity distribution (Daily et al., 1992; Ellis and Oldenburg, 1994; Li and Oldenburg, 1994; Zhang et al., 1995). Recently, ERT surveys have found their way into subsurface hydrology applications. This is attributed to the fact that knowledge of the spatial distribution of the electrical properties of subsurface media can provide valuable information for characterizing waste sites and monitoring flow and contaminant movement in the vadose zone. For example, during an infiltration event, the moisture content of a geological medium is generally the only factor that undergoes dramatic changes, and the changes in resistivity can be related to changes in moisture content. Tracking the changes in resistivity through time, therefore, has been found to be useful for detecting temporal changes in the moisture content of the vadose zone (Daily et al., 1992; Zhou et al., 2001).

While ERT surveys are useful, Yeh et al. (2002) showed that uncertainties in tracking moisture movement in the vadose zone on the basis of ERT surveys alone can be quite significant. They reported that both inverse modeling of the ERT surveys and the spatial variability of the parameters in the constitutive resistivity–moisture content relation contribute to the overall uncertainties. Depending on the number of parameters to be inverted and the quantity of data available (Yeh and imnek, 2002), the ERT survey inverse problem is often ill posed, and may have no unique solution. For ill-posed problems, most ERT inversion approaches employ optimization algorithms with some type of regularization (e.g., Tikhonov, 1963). While the regularization algorithm yields smooth estimates, there is no guarantee that it will produce the best unbiased estimate of the resistivity field, which reflects flow processes and underlying geological structures. Neither does the regularization provide a meaningful way to quantify the uncertainty associated with the spatial variability. Additionally, when converting changes in resistivity to changes in moisture content, it is often assumed that the parameters in the constitutive relation (e.g., Archie's Law) are constant across the entire domain. Recent studies by Baker (2001) and Yeh et al. (2002) reported pronounced spatial variability of these parameters in the field. Neglecting this spatial variability and using a single resistivity–moisture calibration curve can add to the level of uncertainty in the final interpreted change in moisture content (in an unquantifiable way).

Besides the uncertainties inherent in the interpretation of ERT surveys, most current three-dimensional inverse models demand significant computational resources to process the large data sets typically collected during a survey. Furthermore, a change in moisture content only provides qualitative information about water movement in the vadose zone. The actual moisture content values at each point, which are vital to hydrological investigations or hydrological inverse modeling (e.g., Yeh and imnek, 2002), remain unknown. For instance, since the unsaturated hydraulic conductivity is a nonlinear function of moisture content, changes in the moisture content alone do not provide enough information to uniquely characterize the unsaturated hydraulic properties of the vadose zone. As a consequence, an inverse approach is needed that can account for spatial variability of the resistivity–moisture content relation, efficiently process the large number of data sets, and produce detailed moisture content distributions with the least uncertainty.

While the physical process of electric current flow is different from that of groundwater flow, the governing equation of the electric current and the potential fields created during ERT surveys is analogous to that of steady flow in a confined aquifer. The mathematical treatment of the inversion of an ERT survey is therefore similar to that used in hydraulic tomography (Yeh and Liu, 2000; Bohling et al., 2002). Using this similarity, the sequential, successive linear estimator approach for hydraulic tomography (Yeh and Liu, 2000), which has been validated with sandbox experiments (Liu et al., 2002), was extended to ERT surveys by Yeh et al. (2002). However, the work by Yeh et al. (2002) did not directly estimate moisture content distributions, nor consider the spatial variability in the constitutive resistivity–moisture content relation.

We developed an integrative algorithm based on a stochastic information fusion concept (Yeh and imnek, 2002) to provide the best unbiased estimate of moisture content distributions in the vadose zone by fusing both hydrological and geophysical information. The hydrological information includes point measurements of the moisture content and prior information about moisture content distributions (i.e., mean, variance, and correlation structure). The geophysical information consists of point measurements of electric potentials, parameters of the resistivity–moisture content relation, and prior information about spatial variability of these parameters. Two-dimensional numerical experiments were first conducted to evaluate the effect of uncertainties associated with the prior information. Subsequently, numerical experiments were used to demonstrate the robustness of the integrative inverse approach for delineating transient moisture content distributions during a nonuniform infiltration event in a three-dimensional heterogeneous vadose zone.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

 
Governing Equation for the Electric Potential
In a geological formation, the electric current flow induced by an electrical resistivity survey in general can be described by

[1]

where is electric potential (V), I represents the electric current source density per volume (A m^–3), and is the electrical conductivity (S m^–1). Electrical conductivity, , is the reciprocal of the electrical resistivity, (m), which is assumed to be locally isotropic. The boundary conditions associated with Eq. [1] are

[2]

and

[3]

where * is the electric potential specified at boundary ₁, i denotes the electrical current density per unit area (A m^–2), and n is the unit vector normal to the boundary ₂.

Constitutive Resistivity and Moisture Content Relation
In this study, a power law relation was used to relate resistivity to moisture content (e.g., Yeh et al., 2002):

[4]

where is bulk electrical resistivity, ₀ is a fitting parameter that is related to the electrical resistivity of pore water, m is a dimensionless fitting parameter, and denotes volumetric moisture content. We assume that ₀ does not change during an infiltration event in a field. Using Eq. [4], the linear relation between log resistivity before and after an infiltration can be expressed as

[5]

According to this equation, the change of log resistivity (ln) is linearly proportional to the change in log moisture content (ln). If m is spatially invariant, the pattern of ln(x) directly corresponds to the pattern of ln(x) in the entire field. On the other hand, if m is spatially variable and independent of ln(x), then the pattern of change seen in ln(x) does not directly correspond to the pattern of ln (x).

Furthermore, Eq. [5] is derived with the assumption that during an infiltration event ₀ remains constant, which may not always be valid. The resistivity of a porous medium can be highly variable, depending on the degree of saturation and the type of ions present in the pore water (Sharma, 1997). The spatial variability of ₀ and m may directly correspond to the pore water chemistry. Since silica, which comprises most mineral grains (except metallic ores and clays), is an insulator, the observed electrical conduction in porous media is mainly through interstitial pore water. When clay minerals are present, a relatively large number of ions may flow into or out of solution, through ion exchange, thus significantly changing the electrical conductivity of the fluid. During an infiltration event, other chemical reactions or processes may occur due to differences in water chemistry in the infiltrated water, thus altering the composition of ions present in pores and changing the nature of the pore electrolytes. This adds another level of spatial variability to the resistivity distribution.

Baker (2001) measured electrical resistivity as a function of moisture content for core samples collected from a bore hole at the Sandia-Tech Vadose Zone infiltration field site in Socorro, NM. A total of 25 samples were collected from eight 1.52-m (5-foot) continuous cores. The electrical resistivity values of the samples at several moisture contents were determined using an impedance analyzer. Equation [4] was subsequently fitted to the measured resistivity and moisture data to determine the values for ₀ and m. On the basis of analysis of the data set, Yeh et al. (2002) reported that both ln₀ and lnm were approximately normally distributed. The geometric mean of ₀ was 7.036 m and the variance, standard deviation, and coefficient of variation for ln₀ were 0.633, 0.796, and 40.8%, respectively. For m, the geometric mean was 1.336 (dimensionless), while variance, standard deviation, and coefficient of variation for lnm were 0.034, 0.185, and 63.7%, respectively. In addition, they found that both parameters are not entirely disordered in space but correlated over short distances. For ln₀, an exponential variogram model was used to describe its spatial variability with sill, range, and nugget values of 0.8, 3.5 m, and 0.08, respectively. Similarly, an exponential variogram model was used for lnm. The sill, range, and nugget values for lnm were 0.043, 3.5 m, and 0.01, respectively. No significant correlation between ln₀ and lnm was reported.

The observed spatial variability of ln₀ and lnm implies that equivalent changes in moisture content at different locations in the medium may lead to different changes in the measured electrical resistivity. As a result, the pattern of change in resistivity, detected by ERT surveys in a field, may not necessarily reflect the true pattern of change in the moisture content. To overcome this difficulty in interpreting ERT field surveys, we developed an integrative inverse algorithm that is based on the concept of stochastic information fusion developed by Yeh and imnek (2002).

Inverse (or Estimation) Algorithm
Linear Estimator
Since the parameters, ₀, , and m, vary spatially, and their specification at every point in space is practically impossible, they will be treated as random fields (e.g., Yeh, 1992 and 1998), which are characterized by their statistical moments (i.e., means, variances, and correlation scales). Notice that the fields can be either stationary or nonstationary stochastic processes. Accordingly, and the electric potential are also regarded as stochastic processes. In the following analysis, it is assumed that ln ₀ = F + f, ln = A + a, lnm = N + n and = H + h, where F, A, N, and H are the mean values, and f, a, n, and h are the perturbations. The primary goal of the inverse (or estimation) algorithm is to estimate and at any point in the three-dimensional geologic medium, although it can be used to estimate other parameters (i.e., ₀, m, and ) as well. The estimation algorithm integrates point parameter measurements at several locations (including ₀, m, and ) and ERT measurements that consist of and the transmitted current density, I. Prior information about the means and covariances of ₀, m, and is also included; this information could be estimated from core samples, geological well logs, or outcrops (Yeh et al., 2002). The effects of uncertainty in the prior information will be studied below in the Effects of Uncertainties in the Prior Information section.

Using a first-order analysis, a state variable, such as , can be expanded in a Taylor series about the mean values of parameters. Neglecting the second- and higher-order terms of the Taylor series leads to a linear relation between the state variable and the parameters, _i:

[6]

For ERT surveys, _i = _i – _i represents the zero mean perturbation of a natural log transformed parameter, such as f = ln₀ – F, a = ln – A, and n = lnm – N; the zero mean perturbation of the state variable (i.e., electric potential) is h = – = – H; /_i represent the sensitivity derivatives of electric potential with respect to the parameters, and are computed using an adjoint state method. Details of the derivation of these sensitivities can be found in Sun and Yeh (1992), Li and Yeh (1998) and Hughson and Yeh (2000). The sensitivity of the electric potential at location i to a perturbation in a parameter at location k can be generally expressed as

[7]

Specifically, the sensitivities for the parameters in this study are:

[8]

[9]

[10]

where _k is the domain of the element containing node k if a finite element approach is used and represents the adjoint state variable, which can be solved for using the adjoint state equation:

[11]

where is the Dirac delta function and x_k is the measurement location of the electric potential. Notice that the mean electric potential, , is needed to evaluate the sensitivities (see Eq. [8], [9], and [10]). The mean electric potential is derived by solving a mean equation that is of the same form as the Eq. [1] (Yeh et al., 2002), with a mean electrical resistivity and moisture content relation that is the same as Eq. [4], but with the parameters set to their mean values.

Once the mean electric potential field is derived, the above sensitivity equations are used to calculate covariance of h and the cross covariance between h and . Rewriting Eq. [6] in a matrix form in terms of perturbations yields

[12]

where {} indicates the vector of the discretized variable, J_h is a jacobian matrix representing the derivatives of the potential with respect to the parameters (i.e., /_i), which can be obtained using Eq. [8] through [10], and has dimensions of n_h x nelem, where the number of voltage measurement locations is given as n_h, and nelem is the total number of elements in the domain. Multiplying Eq. [12] by the transpose of {} (i.e., {f}, {a}, {n}) and {h}, and then taking the expectations on both sides, yields

[13]

where superscript T denotes transpose; R_h represents the cross-covariance functions between h and f, h, and a, or h and n, with dimensions of n_h x nelem; R denotes the covariance functions of f, a, and n with dimensions of nelem x nelem, and they are given a priori. A nugget can be added to the covariance to represent measurement errors, or variations within the sample scale if they are known. Because little information is available about the parameters f, a, and n in the field, these parameters are assumed to be independent of each other in our study. Such an assumption merely represents the worst scenario in which knowledge of one variable does not provide any information about the others. R_hh in Eq. [13] represents the covariance function of h, which has dimensions of n_h x n_h. Similarly, a nugget representing errors or variability due to scale disparity in potential measurements can be added to this covariance function.

Using these cross-covariance and covariance functions, a first-order estimate of the perturbations of the log-transformed parameters is obtained by using the conditional expectation given observed primary information * (i.e., ₀, m, and ) and secondary information h* (i.e., ) collected during an ERT survey (Priestley, 1989; Yeh and Zhang, 1996):

[14]

where is a nelem long vector of the estimated perturbation of parameters, ₀, m, and ; * and h* are perturbations of measured parameters and electrical potential, respectively; and _h are weights (or cokriging weights in geostatistics) for the measurements, * and h*, respectively. They are evaluated as follows:

[15]

where C represents the covariance of between measurement locations; C_hh represents the covariance of h between measurement locations; and C_h denotes cross covariance between h and at measurement locations. The right-hand side of Eq. [15] consists of the covariance of between measurement locations and the estimate position, and the cross covariance of h at measurement locations and at the estimate location. In other words, the weights used in the estimation are not only directly related to spatial correlation structures of parameters and state variables at measurement locations, but also the cross correlations of the parameters and state variables, and correlation and cross correlation between measurements and the location where the parameter is to be estimated. These matrices all are subsets of the covariance and cross covariance matrices in Eq. [13].

Successive Linear Estimator
Equation [14] approximates the nonlinear relation between the parameter to be estimated and the measured electric potential by means of a linear first-order approximation. Thus, the equation cannot fully exploit electric potential measurements. To circumvent this problem, a successive linear estimator similar to that used by Yeh et al. (1996), Zhang and Yeh (1997), Hanna and Yeh (1998), Vargas-Guzmán and Yeh (1999), Hughson and Yeh (2000), and Vargas-Guzmán and Yeh (2002) is employed. That is,

[16]

where ^(r+1) and ^(r) represent the parameter estimated at iteration r + 1 and r, ^r is the electric potential at the measurement locations calculated from the forward simulation using parameters estimated at iteration r, and _h is the weight at iteration r, which is determined from the following:

[17]

The solution to Eq. [17] requires knowledge of _hh and _h; they are estimated using the following approximations at each iteration:

[18]

where J_h is the sensitivity matrix of n_h x nelem at iteration r, and superscript T stands for the transpose. At iteration r = 0, is given by

[19]

where is a subset of R. For r 1, the residual covariances are evaluated according to

[20]

Notice that these residual covariances represent first-order approximates of the conditional covariances.

Once all three parameters, f, a, and n, are estimated by fully utilizing the potential data and direct measurement of the parameters (if any), the electrical conductivity is then estimated using Eq. [4]. Afterwards, the mean electric potential equation is solved again with the newly estimated electrical conductivity for a new electric potential field. Then, the maximum change of ² (the variance of the estimated parameters of f, a, and n) and the change in the largest potential misfit among all measurement locations between two successive iterations are evaluated. If both changes are smaller than prescribed tolerances, the iteration stops. If not, new values for _h and _hh are evaluated using Eq. [18]. Equation 17 is then solved to obtain a new set of weights that are used in Eq. [16] with [* – ^(r)] to obtain a new estimate of the parameters. A theoretical proof of the convergence of the successive linear estimator is given in Vargas-Guzmán and Yeh (2002).

Sequential Estimation Approach
The previous section describes the ERT inversion algorithm for only one set of primary and secondary information obtained in one DC transmission. This algorithm can simultaneously include all potential measurements collected during all DC transmissions in an ERT survey. However, the system of equations (Eq. [15] and [17]) can become extremely large and ill conditioned, in which case stable solutions to the equations are difficult to obtain (Hughson and Yeh, 2000). To avoid numerical difficulties in solving the large system of equations, the voltage data sets are included sequentially. The sequential algorithm used is similar to the one developed for use in hydraulic tomography inversion (Yeh and Liu, 2000). In essence, the proposed sequential approach uses the estimated electrical conductivity field, the ₀, m, and fields, and their covariances and cross covariances, as prior information for the next estimation using new sets of current–voltage data from DC transmissions at different locations. Vargas-Guzmán and Yeh (1999) and Yeh and imnek (2002) gave an illustrative example of the sequential approach and explained the necessity for updating the covariances and cross covariances. The sequential inclusion of data sets from different DC transmissions continues until all data sets have been utilized. All data sets are fully processed by propagating the conditional first and second moments from one data set to another. Such a sequential approach allows accumulation of high-density secondary information obtained from an ERT survey, while maintaining the covariance matrix at a manageable size that can be solved with minimal numerical difficulties.

Inversions of ERT surveys for environmental applications are generally ill posed since the number of parameters to be estimated is often much greater than the number of measurements of the state variable (see Yeh and imnek, 2002 for a discussion about necessary and sufficient conditions for a well-posed inverse problem). An ill-posed problem has an infinite number of global minima and solutions. Classical inverse algorithms (e.g., regularized least-squares approach, Tikhonov, 1963) can only derive an estimated parameter field that produces an electric potential field honoring measurements at sampling locations and a smooth estimate at other locations. This smooth field, however, does not necessarily honor the characteristics of the spatial variability of the true parameters (e.g., mean, variance, and correlation structure). On the other hand, our sequential–successive linear estimator estimates conditional effective parameters (Yeh et al., 1996; Hanna and Yeh, 1998; Yeh, 1998; Yeh and imnek, 2002). It aims to yield a parameter field that produces not only parameter values and state variables observed at measurement locations, but also conditional effective parameter values at locations where no measurements are available. The effective parameters are our estimates based on the spatial statistics of the parameter fields and their cross correlations with state variables. These cross correlations are evaluated using the governing equation for the electric potential and fluid flow processes. Discussions of advantages of the approach can be found in Yeh and imnek (2002).

	NUMERICAL EXPERIMENTS

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

 
Effects of Uncertainties in the Prior Information
The proposed inverse approach requires input for the prior information such as the mean, variance, and correlation structures of the parameters to be estimated (i.e., ₀, m, and ). In practice, these inputs have to be all estimated and involve uncertainty. The following numerical experiments are devoted to investigating the effect of uncertainties in the mean and covariance functions of on the estimates of the distribution. In these experiments, the other two parameters, ₀ and m, are spatial variables but their spatial distributions are here assumed to be known precisely. This assumption is necessary to reduce the number of cases to be analyzed. The results of this investigation should be applicable to cases where the prior information of these two parameters also involve uncertainty.

The numerical experiments were conducted for a hypothetical, two-dimensional, vertical vadose zone with dimensions of 200 by 200 cm, discretized into 10 horizontal and 20 vertical elements, 200 cm² each. The unsaturated hydraulic properties of each element were assumed to be described by the Mualem–van Genuchten model (van Genuchten, 1980):

[21]

where K() is the unsaturated hydraulic conductivity as a function of the pressure head, ; K_s is the saturated hydraulic conductivity; and ß are shape factors; = (1 – 1/ß); _s is the saturated moisture content; and _r is the residual moisture content. To represent heterogeneity, the parameters of Eq. [21] were assumed to be stochastic processes. Since spatial variations in _s and _r are generally negligible, both were treated as deterministic constants with values of 0.366 and 0.029, respectively. The parameters K_s, , and ß for each element in the simulation domain were generated using a method by Gutjahr (1989) with specified means, variances, and correlation structures as listed in Table 1. These generated random fields are then denoted as the true distributions of the hydraulic parameters for the hypothetical site.

View larger version (33K): [in this window] [in a new window]   Fig. 1. (a) The mean moisture content distribution simulated using homogeneous hydraulic parameters; (b) the true moisture content distribution simulated using heterogeneous hydraulic parameters.

View larger version (35K): [in this window] [in a new window]   Fig. 2. (a) Generated random field for ln0; (b) generated random field for lnm.

View larger version (48K): [in this window] [in a new window]   Fig. 3. (a) True moisture content distribution; (b) estimated moisture content distribution assuming moisture content has a geometric mean and exponential covariance function; (c) estimated moisture content distribution assuming moisture content has a true mean and exponential covariance function; (d) estimated moisture content distribution assuming a true mean and computed covariance function. Circles indicate current source locations, triangles indicate voltage measurement locations, and squares indicate the moisture content measurement location.

Notice that the variance of the parameters appears on both sides of the system of equations (i.e., Eq. [15] and [17]) and will be canceled out. Therefore, the value of the variance has no effect on the solution of the system of equations and in turn, the estimate of parameters. That is, the estimation procedure relies on only the correlation and cross correlation. The variance, however, does affect the predicted uncertainty associated with the estimated parameters.

The first scenario approximates the mean distribution using the geometric mean of the true distribution. Regardless of the fact that the covariance structure of the field is different from those of the hydraulic properties due to flow processes (Yeh et al., 1985a, 1985b, 1985c), in this scenario they were assumed to be the same (i.e., an exponential model with horizontal and vertical correlation scale of 240 and 20 cm, respectively). Using these approximations, together with other necessary input, we estimated perturbations of the distribution, which were then added to the geometric mean to derive the final estimated distribution (Fig. 3b). A comparison between the true and estimated distribution is shown using a scatter plot (Fig. 4a) .

View larger version (16K): [in this window] [in a new window]   Fig. 4. (a) Scatter plot corresponding to Fig. 3b; (b) scatter plot corresponding to Fig. 3c; (c) scatter plot corresponding to Fig. 3d.

The third scenario employed the ensemble mean distribution as input, but used a correct covariance function of instead of that of the hydraulic properties. The covariance function was computed using the first-order approximation:

[22]

where R is the covariance function of the moisture content; R denotes the covariance functions of the log-transformed perturbations of the hydraulic properties of K_s, , and ß, which were assumed known; and J is the Jacobian matrix (details of computing J can be found in Hughson and Yeh, 2000). The final estimated distribution and its scatter plot are shown in Fig. 3d and 4c, respectively.

Comparing Fig. 3b, 3c, and 3d with Fig. 3a reveals that the estimated distributions for the three scenarios resemble the general pattern of the true distribution, despite the different assumptions used in the prior information. Comparisons of scatter plots (Fig. 4) also lead to the same conclusion. To quantify this conclusion, the L1 and L2 norms of these scenarios were evaluated as follows:

[23]

[24]

where _i and _i represent the true and estimated moisture content, respectively. We found that values of the L1 and L2 norms for the three cases differed only slightly: the norms from the use of the ensemble mean are slightly smaller than those from the use of the geometric mean. We therefore conclude that choices of mean and covariance function of moisture content have only a minor effect on the estimate if a sufficient number of potential measurements are available and the ₀ and m distributions are known precisely. This result agrees with the findings by Liu et al., (2002) for hydraulic tomography surveys. That is, abundant data sets collected during tomographic surveys (either ERT or hydraulic tomography surveys) can overcome the uncertainty associated with prior information such as mean, variance, and correlation structure. This result is generally true in these synthetic cases with small domains. Under field conditions, where a site is three-dimensional and encompasses a much larger domain than in the example, and where the ₀ and m fields are unknown, the effects of uncertainty in the input parameters are expected to be greater (Yeh and imnek, 2002). Nevertheless, the reliability of inverse modeling of tomography surveys increases with the spatial density of the secondary information. A dense deployment of secondary information sensors therefore will yield small uncertainty in the result. On the other hand, a sparse deployment would require accurate prior information to facilitate statistically unbiased estimates.

Direct Estimation of Three-Dimensional Moisture Content Fields
In the following numerical examples, we tested our inverse algorithm for transient infiltration and flow through three-dimensional vadose zones. Water movement was simulated in a three-dimensional hypothetical vadose zone at 1000 and 50000 min from the commencement of an infiltration event. The simulated moisture content distributions at these two times, ₁₀₀₀ and ₅₀₀₀₀, were used as the true moisture content fields; and their difference was denoted as the true moisture content change in the following analysis. Following generation of random fields of ₀ and m, the true resistivity fields at 1000 and 50000 min were calculated using Eq. [4] with the generated _0, m, ₁₀₀₀, and ₅₀₀₀₀ fields.

Next, ERT surveys were simulated using these two resistivity fields. After collecting voltage–current data sets, two different inverse approaches were used to interpret water flow due to the infiltration event. The first approach used the Yeh et al. (2002) inverse model based on electric potential measurements only to derive resistivity fields at times of 1000 and 50000 min. Next, the change in resistivity from 1000 to 50000 min is computed. Finally, the estimated change in resistivity is used to interpret the change in moisture content. The second approach estimated moisture distribution at 50000 min directly using our new integrative method.

The hypothetical site was assumed to be a cube, 200 cm on each side, consisting of 2000 elements of 20 by 20 by 10 cm. Random values of the heterogeneous hydraulic parameters (K_s, , and ß) were generated and assigned to each element using the spectral method with the parameters shown in Table 2. Figures 5a, 5b, and 5c show the generated lnK_s, ln, and lnß fields, respectively.

View larger version (64K): [in this window] [in a new window]   Fig. 5. (a) Generated lnKs field; (b) generated ln field; (c) generated lnß field; (d) true moisture content field at t = 1000 min.

View larger version (35K): [in this window] [in a new window]   Fig. 7. (a) True ln from 1000 to 5000 min; (b) estimated ln from 1000 to 50000 min.

View larger version (38K): [in this window] [in a new window]   Fig. 6. (a) Generated true ln0 field; (b) generated true lnm field.

View larger version (53K): [in this window] [in a new window]   Fig. 8. (a) True moisture content at t = 50000 min; (b) estimated at t = 50000 min without measurements; (c) estimated at t = 50000 min with 20 measurements. (d) Schematic diagram for three-dimensional electrical resistivity tomography experiments. Top center square indicates the infiltration area, right triangles indicate voltage measurement locations, circles show current source locations, and squares show the locations of 20 measurements. (e) The scatter plot corresponding to (b); (f) the scatter plot corresponding to (c).

Direct Estimation of Moisture Content Distribution
Figures 8a and b show the true moisture content distribution at 50000 min and the corresponding estimated moisture content distribution using 111 x 5 voltage measurements, and 80 pairs of ₀ and m measurements. The estimated field using the same number of voltage, ₀ and m measurements, but also including 20 direct measurements, is illustrated in Fig. 8c. Comparing the three figures we find that the proposed inverse algorithm reproduces the general pattern of the simulated true moisture content distribution, even though the constitutive relation between resistivity and moisture varies spatially. The inclusion of the 20 measurements, in particular, greatly improves the estimates of the distribution. Figures 8e and 8f are scatter plots corresponding to Fig. 8b and 8c, respectively; a 45° line indicates perfect estimation. The goodness of fit was also evaluated using L1 and L2 norms. The reduction in the L1 and L2 norms from Fig. 8e to 8f suggests that the addition of moisture content measurements dramatically improves the estimate. The conditional variance of the estimate from our stochastic fusion approach can be used to assess the uncertainty associated with the estimate—a smaller conditional variance indicates less uncertainty in the estimate. The conditional variances corresponding to the estimates by using zero and 20 measurements are shown on Fig. 9a and 9b , respectively. Small conditional variances are located close to the four bore holes where secondary information is measured. At locations where moisture content measurements were collected, the conditional variance is zero, indicating that these observations are honored in the inverse model and that no uncertainty exists.

View larger version (37K): [in this window] [in a new window]   Fig. 9. (a) Conditional variance when no measurements are used at t = 50000 min; (b) conditional variance when 20 measurements are used at t = 50000 min. Circles indicate measurement locations.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

Numerical examples illustrate that interpretations of water movement in the subsurface, based only on the estimated resistivity changes, can be misleading because of spatial variability in the – constitutive relation. Numerical examples also demonstrate the ability of the integrative ERT survey inverse model for estimating moisture contents directly. The model yields good estimates at the locations where primary and secondary information is measured. Primary information (i.e., the moisture measurements) contributes significantly to the accuracy of the estimated moisture content distribution. A large number of potential measurements are useful, but they did not dramatically improve the estimate. This further supports the finding by Yeh et al. (2002) that potential measurements alone are inadequate to characterize water flow in vadose zone. Finally, we conclude that fusing geophysical measurements with hydrological information using a stochastic approach is necessary to yield hydrologically realistic results under field conditions and to quantify uncertainty associated with the results.

	ACKNOWLEDGMENTS

 
This research was funded in part by a DOE EMSP96 grant through Sandia National Laboratories (Contract AV-0655#1), a DOE EMSP99 grant through University of Wisconsin (A019493), and NSF and SERDP grant EAR-0229717. We also thank Kris Kuhlman for his technical editing. Many thanks are extended to two reviewers for providing critical comments to improve the paper. In particular, we are greatly indebted to Dr. Andreas Kemna, who provided a meticulous, insightful, constructive, and open-minded review of our manuscript.

	REFERENCES

TOP ABSTRACT INTRODUCTION THEORY NUMERICAL EXPERIMENTS CONCLUSIONS REFERENCES

 

• Baker, K. 2001, Analysis of hydrological and electrical properties at the Sandia-Tech Vadose Zone Facility. M.S. thesis. New Mexico Tech., Socorro, NM.
• Bohling, G.C., X. Zhan, J. Butler, Jr., and L. Zheng. 2002. Steady shape analysis of tomographic pumping tests for characterization of aquifer heterogeneities. Water Resour. Res. 38(12). 1324. doi:10.1029/2001WR001176.
• Daily, W., A. Ramirez, D. LaBrecque, and J. Nitao. 1992. Electrical resistivity tomography of vadose water movement. Water Resour. Res. 28:1429–1442.
• Ellis, R.G., and S.W. Oldenburg. 1994. The pole-pole 3D DC resistivity inverse problem: A conjugate gradient approach. Geophys. J. Int. 119:187–194.
• Gutjahr, A. 1989. Fast Fourier transforms for random field generation. Project Rep. Los Alamos Grant 4-R58-2690R. Dep. Mathematics, New Mexico Tech, Socorro, NM.
• Hanna, S., and T.-C.J. Yeh. 1998. Estimation of co-conditional moments of transmissivity, hydraulic head, and velocity fields. Adv. Water Resour 22:87–93.
• Hughson, D.L., and T.-C.J. Yeh. 2000. An inverse model for three-dimensional flow in variably saturated porous media. Water Resour. Res. 36:829–839.
• Li, B., and T.-C.J. Yeh. 1998. Sensitivity and moment analysis of head in variably saturated regimes. Adv. Water Resour. 21:477–485.
• Li, Y., and D.W. Oldenburg. 1994. Inversion of 3D dc-resistivity data using an approximate inverse mapping. Geophys. J. Int. 116:527–537.
• Liu, S., T.-C. J. Yeh, and R. Gardiner. 2002. Effectiveness of hydraulic tomography: Numerical and sandbox experiments. Water Resour. Res. 38(4). doi:10.1029/2001WR000338.
• Priestley, M.B. 1989. Spectral analysis and time series. Academic Press, San Diego, CA.
• Sharma, P.V. 1997. Environmental and engineering geophysics. Cambridge University Press, Cambridge, UK.
• Srivastava, R., and T.-C.J. Yeh. 1992. A three-dimensional numerical model for water flow and transport of chemically reactive solute through porous media under variably saturated conditions. Adv. Water Resour. 15:275–287.
• Sun, N.-Z., and W.W.-G. Yeh. 1992. A stochastic inverse solution for transient groundwater flow: Parameter identification and reliability analysis. Water Resour. Res. 28:3269–3280.
• Tikhonov, A.N. 1963. Resolution of ill-posed problems and the regularization method. (In Russian.) Dokl. Akad. Nauk SSSR 151:501–504.
• 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.[Abstract/Free Full Text]
• Vargas-Guzmán, J.A., and T.-C.J. Yeh. 1999. Sequential kriging and cokriging: Two powerful geostatistical approaches. Stochastic Environ. Res. Risk Assess. 13:416–435.
• Vargas-Guzmán, J.A., and T.-C.J. Yeh. 2002. The successive linear estimator: A revisit. Adv. Water Resour. 25:773–781.
• Yeh, T.-C.J. 1992. Stochastic modeling of groundwater flow and solute transport in aquifers. J. Hydrol. Processes 6:369–395.
• Yeh, T.-C.J. 1998. Scale issues of heterogeneity in vadose-zone hydrology. In G. Sposito (ed.) Scale dependence and scale invariance in hydrology. Cambridge University Press, Cambridge, UK.
• Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils: A. statistically isotropic media. Water Resour. Res. 21:447–456.
• Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985b. Stochastic analysis of unsaturated flow in heterogeneous soils: B. statistically anisotropic media. Water Resour. Res. 21:457–464.
• Yeh, T.-C.J., L.W. Gelhar, and A.L. Gutjahr. 1985c. Stochastic analysis of unsaturated flow in heterogeneous soils: C. observations and applications. Water Resour. Res. 21:465–471.
• Yeh, T.-C.J., 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.Y. Liu. 2000. Hydraulic tomography: Development of 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 geostatistically based inverse model for electrical resistivity surveys and its applications to vadose zone hydrology. Water Resour. Res. 38(12):1278. doi:10.1029/2001WR001024
• Yeh, T.-C.J., and J. imnek. 2002. Stochastic fusion of information for characterizing and monitoring the vadose zone. Available at www.vadosezonejournal.org. Vadose Zone J. 1:207–221.[Abstract/Free Full Text]
• Yeh, T.-C.J., and J. Zhang. 1996. A geostatistical inverse method for variably saturated flow in the vadose zone. Water Resour. Res. 32:2757–2766.
• Zhang, J., R.L. Mackie, and T. Madden. 1995. 3-D resistivity forward modeling and inversion using conjugate gradients. Geophysics 60:1313–1325.[ISI]