Development and application of a coupled-process parameter inversion model based on the maximum likelihood estimation method

The coupled flow-mass transport inverse problem is formulated using the maximum likelihood estimation concept. An evolutionary computational algorithm, the genetic algorithm, is applied to search for a global or near-global solution. The resulting inverse model allows for flow and transport parameter estimation, based on inversion of spatial and temporal distributions of head and concentration measurements. Numerical experiments using a subset of the three-dimensional tracer tests conducted at the Columbus, Mississippi site are presented to test the model's ability to identify a wide range of parameters and parametrization schemes. The results indicate that the model can be applied to identify zoned parameters of hydraulic conductivity, geostatistical parameters of the hydraulic conductivity field, angle of hydraulic conductivity anisotropy, solute hydrodynamic dispersivity, and sorption parameters. The identification criterion, or objective function residual, is shown to decrease significantly as the complexity of the hydraulic conductivity parametrization is increased. Predictive modeling using the estimated parameters indicated that the geostatistical hydraulic conductivity distribution scheme produced good agreement between simulated and observed heads and concentrations. The genetic algorithm, while providing apparently robust solutions, is found to be considerably less efficient computationally than a quasi-Newton algorithm.     

A combinatorial optimization scheme for parameter structure identification in ground water modeling

This research develops a methodology for parameter structure identification in ground water modeling. For a given set of observations, parameter structure identification seeks to identify the parameter dimension, its corresponding parameter pattern and values. Voronoi tessellation is used to parameterize the unknown distributed parameter into a number of zones. Accordingly, the parameter structure identification problem is equivalent to finding the number and locations as well as the values of the basis points associated with the Voronoi tessellation. A genetic algorithm (GA) is allied with a grid search method and a quasi-Newton algorithm to solve the inverse problem. GA is first used to search for the near-optimal parameter pattern and values. Next, a grid search method and a quasi-Newton algorithm iteratively improve the GA's estimates. Sensitivities of state variables to parameters are calculated by the sensitivity-equation method. MODFLOW and MT3DMS are employed to solve the coupled flow and transport model as well as the derived sensitivity equations. The optimal parameter dimension is determined using criteria based on parameter uncertainty and parameter structure discrimination. Numerical experiments are conducted to demonstrate the proposed methodology, in which the true transmissivity field is characterized by either a continuous distribution or a distribution that can be characterized by zones. We conclude that the optimized transmissivity zones capture the trend and distribution of the true transmissivity field.     

Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging

The assessment of hydraulic conductivity of heterogeneous aquifers is a difficult task using traditional hydrogeological methods (e.g., steady state or transient pumping tests) due to their low spatial resolution. Geophysical measurements performed at the ground surface and in boreholes provide additional information for increasing the resolution and accuracy of the inverted hydraulic conductivity field. We used a stochastic joint inversion of Direct Current (DC) resistivity and self-potential (SP) data plus in situ measurement of the salinity in a downstream well during a synthetic salt tracer experiment to reconstruct the hydraulic conductivity field between two wells. The pilot point parameterization was used to avoid over-parameterization of the inverse problem. Bounds on the model parameters were used to promote a consistent Markov chain Monte Carlo sampling of the model parameters. To evaluate the effectiveness of the joint inversion process, we compared eight cases in which the geophysical data are coupled or not to the in situ sampling of the salinity to map the hydraulic conductivity. We first tested the effectiveness of the inversion of each type of data alone (concentration sampling, self-potential, and DC resistivity), and then we combined the data two by two. We finally combined all the data together to show the value of each type of geophysical data in the joint inversion process because of their different sensitivity map. We also investigated a case in which the data were contaminated with noise and the variogram unknown and inverted stochastically. The results of the inversion revealed that incorporating the self-potential data improves the estimate of hydraulic conductivity field especially when the self-potential data were combined to the salt concentration measurement in the second well or to the time-lapse cross-well electrical resistivity data. Various tests were also performed to quantify the uncertainty in the inverted hydraulic conductivity field.     

Identification of hydraulic conductivity distributions in density dependent flow fields of submarine groundwater discharge modeling using adjoint-state sensitivities

A mathematical optimal control method is developed to identify a hydraulic conductivity distribution in a density dependent flow field. Using a variational method, the adjoint partial differential equations are obtained for the density- dependent state equations used for the saline aquifer water flow. The adjoint equations are numerically solved in through a finite difference method. The developed method is applied to identify the hydraulic conductivity distribution through the numerical solution of an optimal control problem. To demonstrate the effectiveness of the optimal control method, three numerical experiments are conducted with artificial observation data. The results indicate that the developed method has the potential to accurately identify the hydraulic conductivity distribution in a saline water aquifer flow system.     

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.     

Parameter estimation in ensemble based snow data assimilation: A synthetic study

Estimating erroneous parameters in ensemble based snow data assimilation system has been given little attention in the literature. Little is known about the related methods’ effectiveness, performance, and sensitivity to other error sources such as model structural error. This research tackles these questions by running synthetic one-dimensional snow data assimilation with the ensemble Kalman filter (EnKF), in which both state and parameter are simultaneously updated. The first part of the paper investigates the effectiveness of this parameter estimation approach in a perfect-model-structure scenario, and the second part focuses on its dependence on model structure error. The results from first part research demonstrate the advantages of this parameter estimation approach in reducing the systematic error of snow water equivalent (SWE) estimates, and retrieving the correct parameter value. The second part results indicate that, at least in our experiment, there is an evident dependence of parameter search convergence on model structural error. In the imperfect-model-structure run, the parameter search diverges, although it can simulate the state variable well. This result suggest that, good data assimilation performance in estimating state variables is not a sufficient indicator of reliable parameter retrieval in the presence of model structural error. The generality of this conclusion needs to be tested by data assimilation experiments with more complex structural error configurations.     

Simultaneous determination of aquifer parameters and zone structures with fuzzy c-means clustering and meta-heuristic harmony search algorithm

This study proposes an inverse solution algorithm through which both the aquifer parameters and the zone structure of these parameters can be determined based on a given set of observations on piezometric heads. In the zone structure identification problem fuzzy c-means (FCM) clustering method is used. The association of the zone structure with the transmissivity distribution is accomplished through an optimization model. The meta-heuristic harmony search (HS) algorithm, which is conceptualized using the musical process of searching for a perfect state of harmony, is used as an optimization technique. The optimum parameter zone structure is identified based on three criteria which are the residual error, parameter uncertainty, and structure discrimination. A numerical example given in the literature is solved to demonstrate the performance of the proposed algorithm. Also, a sensitivity analysis is performed to test the performance of the HS algorithm for different sets of solution parameters. Results indicate that the proposed solution algorithm is an effective way in the simultaneous identification of aquifer parameters and their corresponding zone structures.     

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.     

Multi-level parameter structure identification for two-phase porous-media flow problems using flexible representations

We consider identification of absolute permeability (hydraulic conductivity) based on time series of pressure data in sparsely distributed wells for two-phase porous-media flow. For this problem, it is impossible to recover all details of the parameter function. On the other hand, a coarser, approximate recovery may be sufficient for many applications. We propose a novel solution approach, based on reparametrization, for such approximate identification of the parameter function. We use a nonlinear, composite representation, which is detached from the computational grid, allowing for a flexible representation of the parameter function at many resolution levels. This is utilized in a sequential multi-level estimation of the parameter function, starting at a coarse resolution, which is then gradually refined. The composite representation is designed to allow for smooth as well as sharp transitions between regions of nearly constant parameter value. Moreover, it facilitates the estimation also of the structure and smoothness of the parameter function itself. As a limiting case, the chosen representation is reduced to a zonation with implicit representation of the interior boundaries that is equivalent to a level-set representation. A motivation for the selected representation and the multi-level estimation is presented in terms of an analysis of sensitivity and nonlinearity. Numerical examples demonstrate identification of coarse-scale features of reference permeability distributions with varying degree of smoothness. Comparisons show how the multi-level strategy stabilize the identification and avoid local minima of the objective function compared to a single-level strategy.     

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.     

Full-Wave Seismic Data Assimilation: Theoretical Background and Recent Advances

The seismological inverse problem has much in common with the data assimilation problem found in meteorology and oceanography. Using the data assimilation methodology, I will formulate the seismological inverse problem for estimating seismic source and Earth structure parameters in the form of weak-constraint generalized inverse, in which the seismic wave equation and the associated initial and boundary conditions are allowed to contain errors. The resulting Euler?CLagrange equations are closely related to the adjoint method and the scattering-integral method, which have been successfully applied in full-3D, full-wave seismic tomography and earthquake source parameter inversions. I will review some recent applications of the full-wave methodology in seismic tomography and seismic source parameter inversions and discuss some challenging issues related to the computational implementation and the effective exploitation of seismic waveform data.     

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.     

A brief review of kriging and its application to optimal interpolation and observation well selection

ABSTRACT

The current state of kriging in subsurface hydrology is critically reviewed. In an application to a region where boreholes already exist, methods of optimal location of additional observation wells for geophysical parameter investigation and optimal interpolation for the purpose of solving the inverse problem are investigated. The particular case of the location of wells for the measurements of transmissivity and hydraulic head in the Kennet Valley Chalk aquifer, UK, is examined. Results of interpolation of measured hydraulic conductivity values by kriging are compared with results from a standard graphical package for interpolation. Reference is also made to the distribution obtained by the inverse method (in which the conductivity distribution is obtained from the head distribution). On the basis of the application, the conditional simulation (in which the generated data are both consistent with field values and the field statistical structure) is deemed to be the best. It is also found that different methods of interpolation give widely different distributions in the case of hydraulic conductivity. It is suggested that the kriged map or conditional map of the transmissivity should serve as the basis for regional discretization to which corrections via the inverse model may be made.     

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.     

Inverse methods in hydrogeology: Evolution and recent trends

Parameter identification is an essential step in constructing a groundwater model. The process of recognizing model parameter values by conditioning on observed data of the state variable is referred to as the inverse problem. A series of inverse methods has been proposed to solve the inverse problem, ranging from trial-and-error manual calibration to the current complex automatic data assimilation algorithms. This paper does not attempt to be another overview paper on inverse models, but rather to analyze and track the evolution of the inverse methods over the last decades, mostly within the realm of hydrogeology, revealing their transformation, motivation and recent trends. Issues confronted by the inverse problem, such as dealing with multiGaussianity and whether or not to preserve the prior statistics are discussed.     

Parameter structure identification using tabu search and simulated annealing

In groundwater modeling the identification of an optimal flow or transport parameter that varies spatially should include both the values and structure of the parameter. However, most existing techniques for parameter identification only consider the parameter values. In this study, the problem of identifying optimal parameter structure is treated as a large combinatorial optimization problem. Two recently developed heuristic search techniques, simulated annealing and tabu search, are used to solve the large combinatorial optimization problem. The effectiveness and flexibility of these two techniques are evaluated and compared with simple grid search and descent search, using preliminary results from one-dimensional examples. Among the techniques examined in this paper, tabu search performs extremely well in terms of the total number of function evaluations required.     

Adaptive multi-scale parameterization for one-dimensional flow in unsaturated porous media

In the analysis of the unsaturated zone, one of the most challenging problems is to use inverse theory in the search for an optimal parameterization of the porous media. Adaptative multi-scale parameterization consists in solving the problem through successive approximations by refining the parameter at the next finer scale all over the domain and stopping the process when the refinement does not induce significant decrease of the objective function any more. In this context, the refinement indicators algorithm provides an adaptive parameterization technique that opens the degrees of freedom in an iterative way driven at first order by the model to locate the discontinuities of the sought parameters. We present a refinement indicators algorithm for adaptive multi-scale parameterization that is applicable to the estimation of multi-dimensional hydraulic parameters in unsaturated soil water flow. Numerical examples are presented which show the efficiency of the algorithm in case of noisy data and missing data.     

Simultaneous identification of parameter,initial condition,and boundary condition in groundwater modelling

In this paper, we perform an inverse method to simultaneously estimate aquifer parameters, initial condition, and boundary conditions in groundwater modelling. The parameter estimation is extended to a complete inverse problem that makes the calibrated groundwater flow model more realistic. The adjoint state method, the gradient search method, and the least square error algorithm are combined to build the optimization procedure. Horizontal two‐dimensional groundwater flow in a confined aquifer is exemplified to demonstrate the correlation between unknowns, the contribution of observation, as well as the suitability of applying the inverse method. The correlation analysis shows the connection between storage coefficient and initial condition. Besides, transmissivity and boundary conditions are also highly correlated. More observations at different location and time are necessary to provide sufficient information. A time series of unsteady head is requested for estimation of storage coefficient and initial condition. Observation near boundary is very effective for boundary condition estimation. The observation at pumping well mostly contributes to the estimation of transmissivity. According to all observations, it is possible to identify parameters, initial condition, and boundary condition simultaneously. Furthermore, the results not only illustrate the traditional assumption of known boundary condition but also initial condition, which may cause an incorrect estimation. Copyright © 2009 John Wiley & Sons, Ltd.     

Comparison between gradient based UCODE_2005 and the ensemble Kalman Filter for transient groundwater flow inverse modeling

Gradient based UCODE_2005 and data assimilation based on the Ensemble Kalman Filter(EnKF) are two different inverse methods. A synthetic two-dimensional flow case with four no-flow boundaries is used to compare the UCODE_2005 with the Ensemble Kalman Filter(EnKF) for their efficiency to inversely calculate and calibrate a hydraulic conductivity field based on hydraulic head data. A zonal, random heterogeneous conductivity field is calibrated by assimilating the time series of heads observed in monitoring wells. The study results indicate that the two inverse methods, UCODE_2005 and EnKF, could be used to calibrate the hydraulic conductivity field to a certain degree. More available observations and information about the conductivity field, more accurate inverse results will be obtained for the UCODE_2005. On the other hand, for a realistic zonal heterogeneous hydraulic conductivity field, EnKF can only efficiently determine the hydraulic conductivity field at the first several assimilated time steps. The results obtained by the UCODE_2005 look better than those by the EnKF. This is possibly due to the fact that the UCODE_2005 uses observed head data at every time step, while EnKF can only use observed heads at first several steps due to the filter divergence problem.