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.This work is supported in part by NSF/ITR Grant EAR-0110289. The first author was additionally supported by scholarships from CONACYT and Instituto de Investigaciones Electricas of Mexico. Additional support was provided by the European Commission under Contract EVK1-CT-1999-00041 (W-SAHaRA-Stochastic Analysis of Well Head Protection and Risk Assessment).     

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.     

On the geostatistical approach to the inverse problem

The geostatistical approach to the inverse problem is discussed with emphasis on the importance of structural analysis. Although the geostatistical approach is occasionally misconstrued as mere cokriging, in fact it consists of two steps: estimation of statistical parameters (“structural analysis”) followed by estimation of the distributed parameter conditional on the observations (“cokriging” or “weighted least squares”). It is argued that in inverse problems, which are algebraically undetermined, the challenge is not so much to reproduce the data as to select an algorithm with the prospect of giving good estimates where there are no observations. The essence of the geostatistical approach is that instead of adjusting a grid-dependent and potentially large number of block conductivities (or other distributed parameters), a small number of structural parameters are fitted to the data. Once this fitting is accomplished, the estimation of block conductivities ensues in a predetermined fashion without fitting of additional parameters. Also, the methodology is compared with a straightforward maximum a posteriori probability estimation method. It is shown that the fundamental differences between the two approaches are: (a) they use different principles to separate the estimation of covariance parameters from the estimation of the spatial variable; (b) the method for covariance parameter estimation in the geostatistical approach produces statistically unbiased estimates of the parameters that are not strongly dependent on the discretization, while the other method is biased and its bias becomes worse by refining the discretization into zones with different conductivity.     

A geostatistical approach to an inverse problem: Identification of geometry and estimate of equivalent conductivities for highly heterogeneous porous media with the differential system method

A methodology for identifying the geometry of different materials in highly heterogeneous porous media in discrete inverse problems (DIP) is described. It applies a geostatistical approach within the differential system method (DSM). DSM calculates conductivity values along an integration path beginning at a point with known conductivity. In aquifers with zero source terms, DSM completely describes the conductivity field through a spatially distributed parameter depending on hydraulic head gradients and integration path. A factor analysis of the structural components of this parameter (i.e. coregionalisation analysis) was carried out to identify the geometry of different materials, corresponding to distinct statistically homogeneous areas. The equivalent conductivity values for homogeneous areas were estimated.This approach was applied for a synthetic aquifer. The identification of geometry was accurate and the estimates of equivalent parameters were good, compared with reference values. The accuracy of the results depended on errors in hydraulic gradients, compared with conductivity gradients.     

A geostatistically based inverse model for three-dimensional variably saturated flow

We present a geostatistically based inverse model for characterizing heterogeneity in parameters of unsaturated hydraulic conductivity for three-dimensional flow. Pressure and moisture content are related to perturbations in hydraulic parameters through cross-covariances, which are calculated to first-order. Sensitivities needed for covariance calculations are derived using the adjoint state sensitivity method. Approximations of the conditional mean parameter fields are then obtained from the cokriging estimator. Correlation between parameters and pressure – moisture content perturbations is seen to be strongly dependent on mean pressure or moisture content. High correlation between parameters and pressure data was obtained under saturated or near saturated flow conditions, providing accurate estimation of saturated hydraulic conductivity, while moisture content measurements provided accurate estimation of the pore size distribution parameter under unsaturated flow conditions.     

A geostatistically based inverse model for three-dimensional variably saturated flow

We present a geostatistically based inverse model for characterizing heterogeneity in parameters of unsaturated hydraulic conductivity for three-dimensional flow. Pressure and moisture content are related to perturbations in hydraulic parameters through cross-covariances, which are calculated to first-order. Sensitivities needed for covariance calculations are derived using the adjoint state sensitivity method. Approximations of the conditional mean parameter fields are then obtained from the cokriging estimator. Correlation between parameters and pressure – moisture content perturbations is seen to be strongly dependent on mean pressure or moisture content. High correlation between parameters and pressure data was obtained under saturated or near saturated flow conditions, providing accurate estimation of saturated hydraulic conductivity, while moisture content measurements provided accurate estimation of the pore size distribution parameter under unsaturated flow conditions.     

Estimation of anisotropic hydraulic conductivity using geophysical data in a coastal aquifer of Karnataka,India

Estimation of hydraulic parameters is essential to understand the interaction between groundwater flow and seawater intrusion. Though several studies have addressed hydraulic parameter estimation, based on pumping tests as well as geophysical methods, not many studies have addressed the problem with clayey formations being present. In this study, a methodology is proposed to estimate anisotropic hydraulic conductivity and porosity values for the coastal aquifer with unconsolidated formations. For this purpose, the one-dimensional resistivity of the aquifer and the groundwater conductivity data are used to estimate porosity at discrete points. The hydraulic conductivity values are estimated by its mutual dependence with porosity and petrophysical parameters. From these estimated values, the bilinear relationship between hydraulic conductivity and aquifer resistivity is established based on the clay content of the sampled formation. The methodology is applied on a coastal aquifer along with the coastal Karnataka, India, which has significant clayey formations embedded in unconsolidated rock. The estimation of hydraulic conductivity values from the established correlations has a correlation coefficient of 0.83 with pumping test data, indicating good reliability of the methodology. The established correlations also enable the estimation of horizontal hydraulic conductivity on two-dimensional resistivity sections, which was not addressed by earlier studies. The inventive approach of using the established bilinear correlations at one-dimensional to two-dimensional resistivity sections is verified by the comparison method. The horizontal hydraulic conductivity agrees with previous findings from inverse modelling. Additionally, this study provides critical insights into the estimation of vertical hydraulic conductivity and an equation is formulated which relates vertical hydraulic conductivity with horizontal. Based on the approach presented, the anisotropic hydraulic conductivity of any type aquifer with embedded clayey formations can be estimated. The anisotropic hydraulic conductivity has the potential to be used as an important input to the groundwater models.     

Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements1

Spectral induced polarization as well as complex electrical measurements are used to estimate, on a non-invasive basis, hydraulic permeability in aquifers. Basic laboratory measurements on a variety of shaly sands, silts and clays showed that the main feature of their conductivity spectra in the frequency range from 10^-3 to 10³ Hertz is a nearly constant phase angle. Thus, a constant-phase-angle model of electrical conductivity is applied to interpret quantitatively surface and borehole spectral induced polarization measurements. The model allows for the calculation of two independent electrical parameters from only one frequency scan and a simple separation of electrical volume and interface effects. The proposed interpretation algorithm yields the true formation factor, the cation exchange capacity and the surface-area-to-porosity ratio, which corresponds to the inverse hydraulic radius. Using a Kozeny–Carman-like equation, the estimation of hydraulic permeability is possible.     

Multiple parameterization for hydraulic conductivity identification

Hydraulic conductivity identification remains a challenging inverse problem in ground water modeling because of the inherent nonuniqueness and lack of flexibility in parameterization methods. This study introduces maximum weighted log-likelihood estimation (MWLLE) along with multiple generalized parameterization (GP) methods to identify hydraulic conductivity and to address nonuniqueness and inflexibility problems in parameterization. A scaling factor for information criteria is suggested to obtain reasonable weights of parameterization methods for the MWLLE and model averaging method. The scaling factor is a statistical parameter relating to a desired significance level in Occam's window and the variance of the chi-squares distribution of the fitting error. Through model averaging with multiple GP methods, the conditional estimate of hydraulic conductivity and its total conditional covariances are calculated. A numerical example illustrates the issue arising from Occam's window in estimating model weights and shows the usefulness of the scaling factor to obtain reasonable model weights. Moreover, the numerical example demonstrates the advantage of using multiple GP methods over the zonation and interpolation methods because GP provides better models in the model averaging method. The methodology is applied to the Alamitos Gap area, California, to identify the hydraulic conductivity field. The results show that the use of the scaling factor is necessary in order to incorporate good parameterization methods and to avoid a dominant parameterization method.     

Estimation of soil hydraulic properties and their uncertainty through the Beerkan infiltration experiment

The Beerkan method based on in situ single‐ring water infiltration experiments along with the relevant specific Beerkan estimation of soil transfer parameters (BEST) algorithm is attractive for simple soil hydraulic characterization. However, the BEST algorithm may lead to erroneous or null values for the saturated hydraulic conductivity and sorptivity especially when there are only few infiltration data points under the transient flow state, either for sandy soil or soils in wet conditions. This study developed an alternative algorithm for analysis of the Beerkan infiltration experiment referred to as BEST‐generalized likelihood uncertainty estimation (GLUE). The proposed method estimates the scale parameters of van Genuchten water retention and Brooks–Corey hydraulic conductivity functions through the GLUE methodology. The GLUE method is a Bayesian Monte Carlo parameter estimation technique that makes use of a likelihood function to measure the goodness‐of‐fit between modelled and observed data. The results showed that using a combination of three different likelihood measurements based on observed transient flow, steady‐state flow and experimental steady‐state infiltration rate made the BEST‐GLUE procedure capable of performing an efficient inverse analysis of Beerkan infiltration experiments. Therefore, it is more applicable for a wider range of soils with contrasting texture, structure, and initial and saturated water content. Copyright © 2015 John Wiley & Sons, Ltd.     

Radial basis function neural network for hydrologic inversion: an appraisal with classical and spatio-temporal geostatistical techniques in the context of site characterization

This paper investigates three techniques for spatial mapping and the consequential hydrologic inversion, using hydraulic conductivity (or transmissivity) and hydraulic head as the geophysical parameters of concern. The data for the study were obtained from the Waste Isolation and Pilot Plant (WIPP) site and surrounding area in the remote Chihuahuan Desert of southeastern New Mexico. The central technique was the Radial Basis Function algorithm for an Artificial Neural Network (RBF-ANN). An appraisal of its performance in light of classical and temporal geostatistical techniques is presented. Our classical geostatistical technique of concern was Ordinary Kriging (OK), while the method of Bayesian Maximum Entropy (BME) constituted an advanced, spatio-temporal mapping technique. A fusion technique for soft or inter-dependent data was developed in this study for use with the neural network. It was observed that the RBF-ANN is capable of hydrologic inversion for transmissivity estimation with features remaining essentially similar to that obtained from kriging. The BME technique, on the other hand, was found to reveal an ability to map localized lows and highs that were otherwise not as apparent in OK or RBF-ANN techniques.     

An interactive multi-objective optimization framework for groundwater inverse modeling

The groundwater inverse problem of estimating heterogeneous groundwater model parameters (hydraulic conductivity in this case) given measurements of aquifer response (such as hydraulic heads) is known to be an ill-posed problem, with multiple parameter values giving similar fits to the aquifer response measurements. This problem is further exacerbated due to the lack of extensive data, typical of most real-world problems. In such cases, it is desirable to incorporate expert knowledge in the estimation process to generate more reasonable estimates. This work presents a novel interactive framework, called the ‘Interactive Multi-Objective Genetic Algorithm’ (IMOGA), to solve the groundwater inverse problem considering different sources of quantitative data as well as qualitative expert knowledge about the site. The IMOGA is unique in that it looks at groundwater model calibration as a multi-objective problem consisting of quantitative objectives – calibration error and regularization – and a ‘qualitative’ objective based on the preference of the geological expert for different spatial characteristics of the conductivity field. All these objectives are then included within a multi-objective genetic algorithm to find multiple solutions that represent the best combination of all quantitative and qualitative objectives. A hypothetical aquifer case-study (based on the test case presented by Freyberg [Freyberg DL. An exercise in ground-water model calibration and prediction. Ground Water 1988;26(3)], for which the ‘true’ parameter values are known, is used as a test case to demonstrate the applicability of this method. It is shown that using automated calibration techniques without using expert interaction leads to parameter values that are not consistent with site-knowledge. Adding expert interaction is shown to not only improve the plausibility of the estimated conductivity fields but also the predictive accuracy of the calibrated model.     

Hydraulic Gradient Comparison Method to Estimate Aquifer Hydraulic Parameters Under Steady-State Conditions

The hydraulic gradient comparison method is an inverse method for estimation of aquifer hydraulic conductivity (or trans-missivity) and boundary conductance for a ground water flow model under steady-state conditions. This method, following formal optimization techniques, defines its objective function to minimize differences between interpreted (observed) and simulated hydraulic gradients, which results in minimization of differences between observed and simulated hydraulic heads. The key features of this method are that (1) the derived optimality conditions have an explicit form with a clear hydrology concept that is con-sistent with Darcy's law, and (2) the derived optimality conditions are spatially independent as they are a function of only local hydraulic conductivity and local hydraulic gradient. This second feature allows a multidimensional optimization problem to be solved by many one-dimensional optimization procedures simultaneously, which results in a substantial reduction in computation time. The results of the numerical performance testing on a heterogeneous hypothetical case confirm that minimizing gradient residuals in the entire model domain leads to minimizing head residuals. Application of the method in real-world projects requires rigorous conceptual model development, use of a global calibration target, and an iterative calibration proess. The conceptual model development includes interpretation of a potentiometric surface and estimation of other hydrologic parameters. This method has been applied to a wide range of real-world modeling projects, including the Rocky Mountain Arsenal and Rocky Flats sites in Colorado, which demonstrates that the method is efficient and practical.     

Effect of heterogeneity on radionuclide retardation in the alluvial aquifer near Yucca Mountain, Nevada

The U.S. Department of Energy is currently studying Yucca Mountain, Nevada, as a potential site for a geological high-level waste repository. In the current conceptual models of radionuclide transport at Yucca Mountain, part of the transport path to pumping locations would be through an alluvial aquifer. Interactions with minerals in the alluvium are expected to retard the downstream migration of radionuclides, thereby delaying arrival times and reducing ground water concentrations. We evaluate the effectiveness of the alluvial aquifer as a transport barrier using the stochastic Lagrangian framework. A transport model is developed to account for physical and chemical heterogeneities and rate-limited mass transfer between mobile and immobile zones. The latter process is caused by small-scale heterogeneity and is thought to control the macroscopic-scale retardation in some field experiments. A geostatistical model for the spatially varying sorption parameters is developed from a site-specific database created from hydrochemical measurements and a calibrated modeling approach (Turner and Pabalan 1999). Transport of neptunium is considered as an example. The results are sensitive to the rate of transfer between mobile and immobile zones, and to spatial variability in the hydraulic conductivity. Chemical heterogeneity has only a small effect, as does correlation between hydraulic conductivity and the neptunium distribution coefficient. These results illustrate how general sensitivities can be explored with modest effort within the Lagrangian framework. Such studies complement and guide the application of more detailed numerical simulations.     

Multi-functional heat pulse probe measurements of coupled vadose zone flow and transport

Simultaneous measurement of coupled water, heat, and solute transport in unsaturated porous media is made possible with the multi-functional heat pulse probe (MFHPP). The probe combines a heat pulse technique for estimating soil heat properties, water flux, and water content with a Wenner array measurement of bulk soil electrical conductivity (EC_bulk). To evaluate the MFHPP, we conducted controlled steady-state flow experiments in a sand column for a wide range of water saturations, flow velocities, and solute concentrations. Flow and transport processes were monitored continuously using the MFHPP. Experimental data were analyzed by inverse modeling of simultaneous water, heat, and solute transport using an adapted HYDRUS-2D model. Various optimization scenarios yielded simultaneous estimation of thermal, solute, and hydraulic parameters and variables, including thermal conductivity, volumetric water content, water flux, and thermal and solute dispersivities. We conclude that the MFHPP holds great promise as an excellent instrument for the continuous monitoring and characterization of the vadose zone.     

Multiscale adaptive estimation of the conductivity field using pump test data

An inverse problem is posed in terms of log-conductivities which are decomposed into macroscale deterministic and microscale stochastic components. The macroscale and microscale conductivities conceptualize hierarchical, scale-dependent aquifer parameters. A deterministic parameter estimation scheme divides a flow domain into a limited number of macroscale constant conductivity zones. A stochastic microscale parameter estimation scheme is used to obtain fluctuations about the macroscale averages in terms of geostatistical models. Both the macroscale and the microscale conductivities are estimated via maximum likelihood, adjoint-state methodologies. Monte Carlo-type approaches are used to examine the distribution of macroscale and microscale conductivity estimates.     

Scaling of flow and transport behavior in heterogeneous groundwater systems

Three-dimensional numerical simulations using a detailed synthetic hydraulic conductivity field developed from geological considerations provide insight into the scaling of subsurface flow and transport processes. Flow and advective transport in the highly resolved heterogeneous field were modeled using massively parallel computers, providing a realistic baseline for evaluation of the impacts of parameter scaling. Upscaling of hydraulic conductivity was performed at a variety of scales using a flexible power law averaging technique. A series of tests were performed to determine the effects of varying the scaling exponent on a number of metrics of flow and transport behavior. Flow and transport simulation on high-performance computers and three-dimensional scientific visualization combine to form a powerful tool for gaining insight into the behavior of complex heterogeneous systems.Many quantitative groundwater models utilize upscaled hydraulic conductivity parameters, either implicitly or explicitly. These parameters are designed to reproduce the bulk flow characteristics at the grid or field scale while not requiring detailed quantification of local-scale conductivity variations. An example from applied groundwater modeling is the common practice of calibrating grid-scale model hydraulic conductivity or transmissivity parameters so as to approximate observed hydraulic head and boundary flux values. Such parameterizations, perhaps with a bulk dispersivity imposed, are then sometimes used to predict transport of reactive or non-reactive solutes. However, this work demonstrates that those parameters that lead to the best upscaling for hydraulic conductivity and head do not necessarily correspond to the best upscaling for prediction of a variety of transport behaviors. This result reflects the fact that transport is strongly impacted by the existence and connectedness of extreme-valued hydraulic conductivities, in contrast to bulk flow which depends more strongly on mean values. It provides motivation for continued research into upscaling methods for transport that directly address advection in heterogeneous porous media.An electronic version of this article is available online at the journal's homepage at http://www.elsevier.nl/locate/advwatres or http://www.elsevier.com/locate/advwatres (see “Special section on vizualization”. The online version contains additional supporting information, graphics, and a 3D animation of simulated particle movement.©1998 Elsevier Science Limited. All rights reserved     

Numerical method of moments for solute transport in a nonstationary flow field conditioned on hydraulic conductivity and head measurements

The unconditional stochastic studies on groundwater flow and solute transport in a nonstationary conductivity field show that the standard deviations of the hydraulic head and solute flux are very large in comparison with their mean values (Zhang et al. in Water Resour Res 36:2107–2120, 2000; Wu et al. in J Hydrol 275:208–228, 2003; Hu et al. in Adv Water Resour 26:513–531, 2003). In this study, we develop a numerical method of moments conditioning on measurements of hydraulic conductivity and head to reduce the variances of the head and the solute flux. A Lagrangian perturbation method is applied to develop the framework for solute transport in a nonstationary flow field. Since analytically derived moments equations are too complicated to solve analytically, a numerical finite difference method is implemented to obtain the solutions. Instead of using an unconditional conductivity field as an input to calculate groundwater velocity, we combine a geostatistical method and a method of moment for flow to conditionally simulate the distributions of head and velocity based on the measurements of hydraulic conductivity and head at some points. The developed theory is applied in several case studies to investigate the influences of the measurements of hydraulic conductivity and/or the hydraulic head on the variances of the predictive head and the solute flux in nonstationary flow fields. The study results show that the conditional calculation will significantly reduce the head variance. Since the hydraulic head measurement points are treated as the interior boundary (Dirichlet boundary) conditions, conditioning on both the hydraulic conductivity and the head measurements is much better than conditioning only on conductivity measurements for reduction of head variance. However, for solute flux, variance reduction by the conditional study is not so significant.     

Stochastic inverse method for estimation of geostatistical representation of hydrogeologic stratigraphy using borehole logs and pressure observations

An approach is presented for identifying statistical characteristics of stratigraphies from borehole and hydraulic data. The approach employs a Markov-chain based geostatistical framework in a stochastic inversion. Borehole data provide information on the stratigraphy while pressure and flux data provide information on the hydraulic performance of the medium. The use of Markov-chain geostatistics as opposed to covariance-based geostatistics can provide a more easily interpreted model geologically and geometrically. The approach hinges on the use of mean facies lengths (negative inverse auto-transition rates) and mean transition lengths (inverse cross-transition rates) as adjustable parameters in the stochastic inversion. Along with an unconstrained Markov-chain model, simplifying constraints to the Markov-chain model, including (1) proportionally-random and (2) symmetric spatial correlations, are evaluated in the stochastic inversion. Sensitivity analyses indicate that the simplifying constraints can facilitate the inversion at the cost of spatial correlation model generality. Inverse analyses demonstrate the feasibility of this approach, indicating that despite some low parameter sensitivities, all adjustable parameters do converge for a sufficient number of ensemble realizations towards their “true” values. This paper extends the approach presented in Harp et al. (doi:, 2008) to (1) statistically characterize the hydraulic response of a geostatistical model, thereby incorporating an uncertainty analysis directly in the inverse method, (2) demonstrate that a gradient-based optimization strategy is sufficient, thereby providing relative computational efficiency compared to global optimization strategies, (3) demonstrate that the approach can be extended to a 3-D analysis, and (4) introduce the use of mean facies lengths and mean transition lengths as adjustable parameters in a geostatistical inversion, thereby allowing the approach to be extended to greater than two category Markov-chain models.