Geostatistical Inverse Research Articles

A New Inverse Modeling Approach for Hydraulic Conductivity Estimation Based on Gaussian Mixtures

AbstractThis study proposes a new inverse algorithm to estimate the hydraulic conductivity (K) distribution based on a Gaussian Mixture Model that significantly reduces the number of parameters to be estimated during the inversion process. Moreover, a new objective function that increases the sensitivity of parameters using the spatial derivatives of hydraulic heads is introduced, and the algorithm is further improved by including a Bayes estimator that takes advantage of different possible solutions. The developed approach is tested through multiple synthetic experiments consisting of 250 randomly generated K fields resulting in different levels of heterogeneity and the use of different number of pumping tests, with a total of 1,000 cases of two‐dimensional configuration. A large number of cases are considered to ensure that our findings and conclusions are not based on a single realization. Results revealed significant improvements to K estimates, computational time, and predictions of independently conducted tests not used in the calibration effort when compared to a geostatistical inverse approach. Overall, our results reveal that the Gaussian Mixture inversion approach is able to achieve similar or higher levels of accuracy using half of the pumping tests and 20% of the computational time compared to a geostatistical inversion approach.

3D geostatistical inversion of induced polarization data and its application to coal seam fires

We applied the principal component geostatistical approach (PCGA) to the inversion of time-domain induced polarization data in terms of resistivity and chargeability distributions. The PCGA presents two major advantages over standard methods: (1) It avoids the storage of the usually large covariance matrix, which contains the geostatistical information, by factorizing it in a product of low-rank matrices. (2) It does not assemble the Jacobian matrix per se. We determine the robustness of this approach with three examples. We first reconstruct the electrical conductivity and chargeability fields of two synthetic models generated using the geostatistical software Stanford Geostatistical Modeling Software. The PCGA approach performs better than the Tikhonov-based regularization approach when the true fields are very heterogeneous and the amount of data is limited. The third example is devoted to a field study over the former Lewis coal mine in Colorado (USA). We perform a 3D localization of the burning front of this coal seam fire by applying our geostatistical inverse methodology to a time-domain induced-polarization data set. In this case, the horizontal components of the semivariogram are determined from a self-potential map and the correlation length scale for the vertical component is determined from the known thickness of the coal bed. The tomogram presents a high normalized chargeability associated with the burning front. We evaluate the high normalized chargeability of the burning front in terms of the physical mechanism associated with the cation exchange capacity of the coal and the effect of temperature. This demonstrates the potential of the geostatistical inversion and its suitability for inverting geophysical data, especially when the data density is sparse. In the case of coal seam fires, we determine the suitability of the induced polarization method to localize the burning front and the effect of temperature on the normalized chargeability.

A Reduced‐Order Successive Linear Estimator for Geostatistical Inversion and its Application in Hydraulic Tomography

AbstractHydraulic tomography (HT) is a recently developed technology for characterizing high‐resolution, site‐specific heterogeneity using hydraulic data (nd) from a series of cross‐hole pumping tests. To properly account for the subsurface heterogeneity and to flexibly incorporate additional information, geostatistical inverse models, which permit a large number of spatially correlated unknowns (ny), are frequently used to interpret the collected data. However, the memory storage requirements for the covariance of the unknowns (ny × ny) in these models are prodigious for large‐scale 3‐D problems. Moreover, the sensitivity evaluation is often computationally intensive using traditional difference method (ny forward runs). Although employment of the adjoint method can reduce the cost to nd forward runs, the adjoint model requires intrusive coding effort. In order to resolve these issues, this paper presents a Reduced‐Order Successive Linear Estimator (ROSLE) for analyzing HT data. This new estimator approximates the covariance of the unknowns using Karhunen‐Loeve Expansion (KLE) truncated to nkl order, and it calculates the directional sensitivities (in the directions of nkl eigenvectors) to form the covariance and cross‐covariance used in the Successive Linear Estimator (SLE). In addition, the covariance of unknowns is updated every iteration by updating the eigenvalues and eigenfunctions. The computational advantages of the proposed algorithm are demonstrated through numerical experiments and a 3‐D transient HT analysis of data from a highly heterogeneous field site.

Three-dimensional imaging of aquifer and aquitard heterogeneity via transient hydraulic tomography at a highly heterogeneous field site

Previous studies have shown that geostatistics-based transient hydraulic tomography (THT) is robust for subsurface heterogeneity characterization through the joint inverse modeling of multiple pumping tests. However, the hydraulic conductivity (K) and specific storage (Ss) estimates can be smooth or even erroneous for areas where pumping/observation densities are low. This renders the imaging of interlayer and intralayer heterogeneity of highly contrasting materials including their unit boundaries difficult. In this study, we further test the performance of THT by utilizing existing and newly collected pumping test data of longer durations that showed drawdown responses in both aquifer and aquitard units at a field site underlain by a highly heterogeneous glaciofluvial deposit. The robust performance of the THT is highlighted through the comparison of different degrees of model parameterization including: (1) the effective parameter approach; (2) the geological zonation approach relying on borehole logs; and (3) the geostatistical inversion approach considering different prior information (with/without geological data). Results reveal that the simultaneous analysis of eight pumping tests with the geostatistical inverse model yields the best results in terms of model calibration and validation. We also find that the joint interpretation of long-term drawdown data from aquifer and aquitard units is necessary in mapping their full heterogeneous patterns including intralayer variabilities. Moreover, as geological data are included as prior information in the geostatistics-based THT analysis, the estimated K values increasingly reflect the vertical distribution patterns of permeameter-estimated K in both aquifer and aquitard units. Finally, the comparison of various THT approaches reveals that differences in the estimated K and Ss tomograms result in significantly different transient drawdown predictions at observation ports.

Geostatistical seismic Amplitude‐versus‐angle inversion

ABSTRACTSeismic inversion plays an important role in reservoir modelling and characterisation due to its potential for assessing the spatial distribution of the sub‐surface petro‐elastic properties. Seismic amplitude‐versus‐angle inversion methodologies allow to retrieve P‐wave and S‐wave velocities and density individually allowing a better characterisation of existing litho‐fluid facies. We present an iterative geostatistical seismic amplitude‐versus‐angle inversion algorithm that inverts pre‐stack seismic data, sorted by angle gather, directly for: density; P‐wave; and S‐wave velocity models. The proposed iterative geostatistical inverse procedure is based on the use of stochastic sequential simulation and co‐simulation algorithms as the perturbation technique of the model parametre space; and the use of a genetic algorithm as a global optimiser to make the simulated elastic models converge from iteration to iteration. All the elastic models simulated during the iterative procedure honour the marginal prior distributions of P‐wave velocity, S‐wave velocity and density estimated from the available well‐log data, and the corresponding joint distributions between density versus P‐wave velocity and P‐wave versus S‐wave velocity. We successfully tested and implemented the proposed inversion procedure on a pre‐stack synthetic dataset, built from a real reservoir, and on a real pre‐stack seismic dataset acquired over a deep‐water gas reservoir. In both cases the results show a good convergence between real and synthetic seismic and reliable high‐resolution elastic sub‐surface Earth models.

On the importance of geological data for three-dimensional steady-state hydraulic tomography analysis at a highly heterogeneous aquifer-aquitard system

Hydraulic tomography (HT) has been shown to map subsurface heterogeneity accurately through the joint interpretation of multiple pumping tests. Previous research has shown that smooth hydraulic conductivity (K) estimates are obtained beyond where pumping/observation data are available using the geostatistical inversion approach, when the inversion begins with a homogeneous K and when data densities are not high. However, geological data are typically available through outcrops and borehole logs to provide geological variability. Therefore, we investigate the usefulness of geological data for HT analysis at a highly heterogeneous field site by: (1) comparing calibrated geological models of two different resolutions to two homogeneous and four highly parameterized geostatistical inverse models, in terms of both model calibration and validation performances as well as correspondence of estimated K values with permeameter-estimated K profiles along boreholes; and (2) using geological models as prior information for the geostatistical inversion approach. Results reveal that the simultaneous calibration of geological models to seven pumping test data yields K values that correctly reflect the general patterns of vertical distributions of permeameter-estimated K. We also find that the geostatistical inversion approach using a geological model as prior information performs better for both model calibration and validation than using a homogenous K as a prior, and more importantly, improves the correspondence of K estimates to permeameter test results along wells, as well as in preserving geological features where drawdown measurements are lacking. Overall, our results suggest the joint use of both geological and pumping test data for HT analysis when accurate geological data are available.

Integration of well data into geostatistical seismic amplitude variation with angle inversion for facies estimation

We have developed a new iterative geostatistical seismic amplitude variation with angle (AVA) inversion algorithm that inverts prestack seismic data, sorted by angle gathers, directly for high-resolution density, P-wave velocity, S-wave velocity, and facies models. This novel iterative geostatistical inverse procedure is based on two key main principles: the use of stochastic sequential simulation and cosimulation as the perturbation technique of the model parameter space and a global optimizer based on a crossover genetic algorithm to converge the simulated earth models toward an objective function, in this case, the mismatch between the recorded and synthetic prestack seismic data. As a geostatistical approach, all the elastic models simulated during the iterative procedure honors the well-log data at its own locations, the marginal prior distributions of P-wave velocity and S-wave velocity, and density estimated from the available well-log data, and the corresponding joint distributions between density versus P-wave velocity and P-wave versus S-wave velocity. We successfully tested and implemented this new algorithm on a synthetic prestack data set that mimicked the main properties of a real reservoir, and on a real seismic data set acquired over a deepwater turbidite oil reservoir. In both cases, the results showed a good convergence between the recorded and synthetic seismic. The synthetic example showed high-resolution inverted petroelastic models that reproduced the true petroelastic models. The inverted petroelastic models retrieved from the real case study found high resolution and do agree with previous seismic reservoir characterization studies.

Should hydraulic tomography data be interpreted using geostatistical inverse modeling? A laboratory sandbox investigation

AbstractThe robust performance of hydraulic tomography (HT) based on geostatistics has been demonstrated through numerous synthetic, laboratory, and field studies. While geostatistical inverse methods offer many advantages, one key disadvantage is its highly parameterized nature, which renders it computationally intensive for large‐scale problems. Another issue is that geostatistics‐based HT may produce overly smooth images of subsurface heterogeneity when there are few monitoring interval data. Therefore, some may question the utility of the geostatistical inversion approach in certain situations and seek alternative approaches. To investigate these issues, we simultaneously calibrated different groundwater models with varying subsurface conceptualizations and parameter resolutions using a laboratory sandbox aquifer. The compared models included: (1) isotropic and anisotropic effective parameter models; (2) a heterogeneous model that faithfully represents the geological features; and (3) a heterogeneous model based on geostatistical inverse modeling. The performance of these models was assessed by quantitatively examining the results from model calibration and validation. Calibration data consisted of steady state drawdown data from eight pumping tests and validation data consisted of data from 16 separate pumping tests not used in the calibration effort. Results revealed that the geostatistical inversion approach performed the best among the approaches compared, although the geological model that faithfully represented stratigraphy came a close second. In addition, when the number of pumping tests available for inverse modeling was small, the geological modeling approach yielded more robust validation results. This suggests that better knowledge of stratigraphy obtained via geophysics or other means may contribute to improved results for HT.

Fast iterative implementation of large-scale nonlinear geostatistical inverse modeling

[1] In nonlinear geostatistical inverse problems, it often takes a significant amount of computational cost to form linear geostatistical inversion systems by linearizing the forward model. More specifically, the storage cost associated with the sensitivity matrix H (m × n, where m and n are the numbers of measurements and unknowns, respectively) is high, especially when both m and n are large in for instance, 3-D tomography problems. In this research, instead of explicitly forming and directly solving the linear geostatistical inversion system, we use MINRES, a Krylov subspace method, to solve it iteratively. During each iteration in MINRES, we only compute the products Hx and HTx for any appropriately sized vectors x, for which we solve the forward problem twice. As a result, we reduce the memory requirement from O(mn) to O(m)+O(n). This iterative methodology is combined with the Bayesian inverse method in Kitanidis (1996) to solve large-scale inversion problems. The computational advantages of our methodology are demonstrated using a large-scale 3-D numerical hydraulic tomography problem with transient pressure measurements (250,000 unknowns and ∼100,000 measurements). In this case, ∼200 GB of memory would otherwise be required to fully compute and store the sensitivity matrix H at each Newton step during optimization. The CPU cost can also be significantly reduced in terms of the total number of forward simulations. In the end, we discuss potential extension of the methodology to other geostatistical methods such as the Successive Linear Estimator.

Role of model selection criteria in geostatistical inverse estimation of statistical data‐ and model‐parameters

We analyze theoretically the ability of model quality (sometimes termed information or discrimination) criteria such as the negative log likelihood NLL, Bayesian criteria BIC and KIC and information theoretic criteria AIC, AICc, and HIC to estimate (1) the parameter vector of the variogram of hydraulic log conductivity (Y = ln K), and (2) statistical parameters and proportional to head and log conductivity measurement error variances, respectively, in the context of geostatistical groundwater flow inversion. Our analysis extends the work of Hernandez et al. (2003, 2006) and Riva et al. (2009), who developed nonlinear stochastic inverse algorithms that allow conditioning estimates of steady state and transient hydraulic heads, fluxes and their associated uncertainty on information about conductivity and head data collected in a randomly heterogeneous confined aquifer. Their algorithms are based on recursive numerical approximations of exact nonlocal conditional equations describing the mean and (co)variance of groundwater flow. Log conductivity is parameterized geostatistically based on measured values at discrete locations and unknown values at discrete “pilot points.” Optionally, the maximum likelihood function on which the inverse estimation of Y at pilot points is based may include a regularization term reflecting prior information about Y. The relative weight assigned to this term and its components and , as well as are evaluated separately from other model parameters to avoid bias and instability. This evaluation is done on the basis of criteria such as NLL, KIC, BIC, HIC, AIC, and AICc. We demonstrate theoretically that, whereas all these six criteria make it possible to estimate , KIC alone allows one to estimate validly and (and thus ). We illustrate this discriminatory power of KIC numerically by using a differential evolution genetic search algorithm to minimize it in the context of a two‐dimensional steady state groundwater flow problem. We find that whereas , and the integral scale of Y can be estimated on the basis of a zero‐order mean flow equation, the sill of the Y‐variogram is estimated more accurately by a second‐order approximation of flow. This notwithstanding, KIC prefers the simpler zero‐order moment over the more complex second‐order version.

Exact sensitivity matrix and influence of the number of pilot points in the geostatistical inversion of moment equations of groundwater flow

We present novel equations for the exact sensitivity matrix of the (ensemble) mean hydraulic head under steady state groundwater flow conditions. These equations are embedded in a geostatistical inverse procedure to condition approximations of stochastic moment equations of flow on measured hydraulic conductivities and heads. Our formulation allows considerable improvement of the methodology proposed by Hernandez et al. (2003, 2006) and renders the inversion of moment equations feasible for a large number of unknown hydraulic parameters. The spatial distribution of the natural logarithm, Y, of conductivity is parameterized within the pilot points framework. Whereas prior values of Y at pilot points are obtained by a variant of kriging, posterior estimates at pilot points are obtained through a maximum likelihood fit of computed to measured heads. The maximum likelihood function also includes a regularization term. By means of a synthetic example and upon adopting formal model information criteria we explore the influence of (1) the number of pilot points and (2) the order of approximation of the governing mean flow equation on our ability to properly estimate the log conductivity and head fields and identify the relative weight of the regularization term and the parameters of the underlying Y variogram. We find that none of the adopted information criteria can identify the optimum number of pilot points and the plausibility weight and variogram parameters values can be determined by the Kashyap's Bayesian measure.

Stochastic characterization of the Montalto Uffugo research site (Italy) by geostatistical inversion of moment equations of groundwater flow

Summary We assess the applicability and performance of a methodology of inverting stochastic mean groundwater flow equations to characterize the spatial variability of (natural) log-transmissivity (Y) of a heterogeneous aquifer. The methodology, originally proposed by Hernandez et al., 2003 , Hernandez et al., 2006 , relies on a nonlinear geostatistical inverse algorithm for recursive approximations of steady-state (ensemble) mean groundwater flow that allows estimating jointly the spatial variability of Y, the underlying variogram parameters, and the variance–covariance of the estimates. Estimates of prediction errors of hydraulic heads and fluxes are then calculated a posteriori, upon solving equations satisfied by the corresponding covariances. Here, we extend the methodology to quasi-steady state flow conditions and present its first field application by using information collected during a pumping test performed at the Montalto Uffugo research site (Italy). Log-transmissivity is parameterized geostatistically on the basis of an available measured value and a set of unknown values at discrete pilot points. Best estimates of Y at the measurement location and at the pilot points are obtained by a maximum likelihood fit of computed and measured heads. These posterior estimates are then projected onto the computational grid by kriging. Information on head drawdowns is provided through self-potential signals recorded by 47 surface electrodes during the test. The maximum likelihood-based objective function includes a regularization term reflecting prior information about Y. The relative weight assigned to this term is evaluated separately from other model parameters to avoid bias and instability. We explore the effectiveness of both a zero- and a second-order closure of the mean flow equation at providing a proper geostatistical characterization of the log-transmissivity field. The parameters of the variogram of Y are estimated a posteriori using formal model selection criteria. The adoption of a second-order mean flow model renders the sharpest definition of the regularization term and of the Y variogram parameters.

Use of steady-state concentration measurements in geostatistical inversion

In geostatistical inverse modeling, hydrogeological parameters, such as hydraulic conductivity, are estimated as spatial fields. Upon discretization this results in several thousand (log-)hydraulic conductivity values to be estimated. Common inversion schemes rely on gradient-based parameter estimation methods which require the sensitivity of all measurements with respect to all parameters. Point-like measurements of steady-state concentration in aquifers are generally not well suited for gradient-based methods, because typical plumes exhibit only a very narrow fringe at which the concentration decreases from a maximal value to zero. Only here the sensitivity of concentration with respect to hydraulic conductivity significantly differs from zero. Thus, if point-like measurements of steady-state concentration do not lie in this narrow fringe, their sensitivity with respect to hydraulic conductivity is zero. Observations of concentrations averaged over a larger control volume, by contrast, show a more regular sensitivity pattern. We thus suggest artificially increasing the sampling volume of steady-state concentration measurements for the evaluation of sensitivities in early stages of an iterative parameter estimation scheme. We present criteria for the extent of artificially increasing the sampling volume and for decreasing it when the simulation results converge to the measurements. By this procedure, we achieve high stability in geostatistical inversion of steady-state concentration measurements. The uncertainty of the estimated parameter fields is evaluated by generating conditional realizations.

Geostatistical inversing for large-contrast transmissivity fields

The estimation of field parameters, such as transmissivity, is an important part of groundwater modeling. This work deals with the quasilinear geostatistical inverse approach to the estimation of the transmissivity fields from hydraulic head measurements. The standard quasilinear approach is an iterative method consisting of successive linearizations. We examine a synthetic case to evaluate the basic methodology and some modifications and extensions. The first objective is to evaluate the performance of the quasilinear approach when applied to strongly heterogeneous (or “high-contrast”) transmissivity fields and, when needed, to propose improvements that allow the solution of such problems. For large-contrast cases, the standard quasilinear method often fails to converge. However, by introducing a derivative-free line search as a polishing step after each Gauss–Newton iteration, we have found that convergence can be practically assured. Another issue is that the quasilinear procedure, which uses linearization about the best estimate to evaluate estimation variances, may lead to inaccurate estimation of the variance of the estimated variable. Our numerical results suggest that this may not be a particularly serious problem, though it is hard to say whether this conclusion will apply to other cases. Nevertheless, since the quasilinear approach is an approximation, we propose a potentially more accurate but computer-intensive Markov Chain Monte Carlo (MCMC) procedure based on conditional realizations generated through the quasilinear approach and accepted or rejected according to the Metropolis–Hastings algorithm. Six transmissivity fields with increasing contrast were generated and one thousand conditional realizations were computed for each studied case. The MCMC procedure proposed in this work gives an overall more accurate picture than the quasilinear approach but at a considerably higher computational cost.

Two‐dimensional characterization of hydraulic heterogeneity by multiple pumping tests

The conventional analysis of pumping tests by type‐curve methods is based on the assumption of a homogeneous aquifer. Applying these techniques to pumping test data from real heterogeneous aquifers leads to estimates of the hydraulic parameters that depend on the choice of the pumping and observation well positions. In this paper, we test whether these values may be viewed as pseudo‐local values of transmissivity and storativity, which can be interpolated by kriging. We compare such estimates to those obtained by geostatistical inverse modeling, where heterogeneity is assumed in all stages of estimation. We use drawdown data from multiple pumping tests conducted at the test site in Krauthausen, Germany. The geometric mean values of transmissivity and storativity determined by type‐curve analysis are very close to those obtained by geostatistical inversion, but the conventional approach failed to resolve the spatial variability of transmissivity. In contrast, the estimate from geostatistical inversion reveals more structure. This indicates that the estimates of the type‐curve approaches can not be treated as pseudo‐local values. Concerning storativity, both analysis methods show strong fluctuations. Because the variability of all terms making up the storativity is small, we believe that the estimated variability of storativity is biased. We examine the influence of measurement error on estimating structural parameters of covariance functions in the inversion. We obtain larger correlation lengths and smaller prior variances if we trust the measured data less.

A Bayesian geostatistical transfer function approach to tracer test analysis

Reactive transport modeling is often used in support of bioremediation and chemical treatment planning and design. There remains a pressing need for practical and efficient models that do not require (or assume attainable) the high level of characterization needed by complex numerical models. We focus on a linear systems or transfer function approach to the problem of reactive tracer transport in a heterogeneous saprolite aquifer. Transfer functions are obtained through the Bayesian geostatistical inverse method applied to tracer injection histories and breakthrough curves. We employ nonparametric transfer functions, which require minimal assumptions about shape and structure. The resulting flexibility empowers the data to determine the nature of the transfer function with minimal prior assumptions. Nonnegativity is enforced through a reflected Brownian motion stochastic model. The inverse method enables us to quantify uncertainty and to generate conditional realizations of the transfer function. Complex information about a hydrogeologic system is distilled into a relatively simple but rigorously obtained function that describes the transport behavior of the system between two wells. The resulting transfer functions are valuable in reactive transport models based on traveltime and streamline methods. The information contained in the data, particularly in the case of strong heterogeneity, is not overextended but is fully used. This is the first application of Bayesian geostatistical inversion to transfer functions in hydrogeology but the methodology can be extended to any linear system.

Inverse stochastic moment analysis of steady state flow in randomly heterogeneous media

Nonlocal stochastic moment equations have been used successfully to analyze steady state and transient flow in randomly heterogeneous media conditional on measured values of medium properties. We present a nonlinear geostatistical inverse algorithm for steady state flow that makes it possible to further condition such analyses on measured values of state variables. Our approach accounts for all scales of spatial variability resolvable by the computational grid. It is based on recursive finite element approximations of exact nonlocal first and second conditional moment equations. Hydraulic conductivity is parameterized geostatistically on the basis of measured values at discrete locations and unknown values at discrete “pilot points.” Prior estimates of pilot point values are obtained (optionally) by generalized kriging. Posterior estimates at pilot points and (optionally) at measurement points are obtained by calibrating mean flow against measured values of head. The parameters are projected onto a computational grid via kriging. Maximum likelihood calibration allows estimating not only hydraulic but also (optionally) unknown variogram parameters with or without prior information about the former. The approach yields covariance matrices for parameter estimation as well as head and flux prediction errors, the latter being obtained a posteriori. We implement our inverse approach on highly and mildly heterogeneous media under superimposed mean uniform and convergent flows in a bounded two‐dimensional domain. Our examples illustrate how conductivity and head data act separately and jointly to reduce parameter estimation errors and to model predictive uncertainty. We also evaluate the functional form of the log conductivity variogram and its parameters using likelihood‐based model discrimination criteria.

Conditioning mean steady state flow on hydraulic head and conductivity through geostatistical inversion

Nonlocal moment equations allow one to render deterministically optimum predictions of flow in randomly heterogeneous media and to assess predictive uncertainty conditional on measured values of medium properties. We present a geostatistical inverse algorithm for steady-state flow that makes it possible to further condition such predictions and assessments on measured values of hydraulic head (and/or flux). Our algorithm is based on recursive finite-element approximations of exact first and second conditional moment equations. Hydraulic conductivity is parameterized via universal kriging based on unknown values at pilot points and (optionally) measured values at other discrete locations. Optimum unbiased inverse estimates of natural log hydraulic conductivity, head and flux are obtained by minimizing a residual criterion using the Levenberg-Marquardt algorithm. We illustrate the method for superimposed mean uniform and convergent flows in a bounded two-dimensional domain. Our examples illustrate how conductivity and head data act separately or jointly to reduce parameter estimation errors and model predictive uncertainty.

Geostatistical Inversion of Cross‐Hole Pumping Tests for Identifying Preferential Flow Channels Within a Shear Zone

AbstractCross‐hole pumping tests within a subvertical shear zone in granite were interpreted using a geostatistical inverse approach. The objective was to identify flow channels that influence tracer experiments. The spatial correlation structure of transmissivity (T) was unknown. Therefore, a wide variety of geostatistical hypotheses were made for estimating a total of 40 T fields, that honor equally well the pumping tests, the static head measurements, and the point T data. All estimated T fields reveal channels with similar topology and with a preferred horizontal orientation, possibly due to vertical displacements within the subvertical shear zone. This shows that inverse geostatistical modeling of abundant hydraulic measurements can be useful to reveal a coarse heterogeneity structure, even when data are insufficient to identify the geostatistical structure. However, numerical simulations of a field tracer test were highly sensitive to the geostatistical hypotheses. In summary, T fields leading to successful hydraulic data fitting do not necessarily lead to successful transport predictions, even when they appear to reproduce the coarse heterogeneity structure.

A Geostatistical Inverse Method for Variably Saturated Flow in the Vadose Zone

A geostatistical inverse technique utilizing both primary and secondary information is developed to estimate conditional means of unsaturated hydraulic conductivity parameters (saturated hydraulic conductivity and pore size distribution parameters) in the vadose zone. Measurements of saturated hydraulic conductivity and pore size distribution parameters are considered as the primary information, while measurements of steady state flow processes (soil‐water pressure head and degree of saturation) are regarded as the secondary information. This inverse approach relies on the classical linear predictor (cokriging) theory and takes the advantage of the spatial cross correlation between the soil‐water pressure head and each of the following: degree of saturation, saturated hydraulic conductivity, and a pore size distribution parameter. Using an approximate perturbation solution for steady, variably saturated flow under general boundary conditions, the cross covariances between the primary and secondary information are derived. The approximate solution is formulated on the basis of a first‐order Taylor series expansion of a discretized finite element equation. The sensitivity matrix in the solution is evaluated by an adjoint state sensitivity approach for flow in heterogeneous media under variably saturated conditions. Through several numerical examples the inverse model demonstrates its ability to improve the estimates of the spatial distribution of saturated hydraulic conductivity and pore size distribution parameters using the secondary information.