Image-guided inversion in steady-state hydraulic tomography
Introduction
Hydraulic conductivity (or transmissivity) distribution is of course the holy grail of subsurface hydrogeology. Adequate knowledge of this parameter is required to model groundwater flow and the migrations of contaminants in aquifers. Hydraulic conductivity or transmissivity are usually determined by performing pumping tests [2], [7], [39], [49], [54]. These tests consist in pumping water in one well (or one well interval) and monitoring the water level changes in a set of piezometers or well intervals. This analysis gives an averaged and localized estimate of the transmissivity in the vicinity of the wells. This estimate can be however erroneous in the case of heterogeneous aquifers (e.g., [33], [57]).
Hydraulic tomography has been proposed as a method to investigate the distributions of the hydraulic parameters of heterogeneous aquifers. Hydraulic tomography consists in stimulating an aquifer by performing a sequence of pumping tests or injection/extraction tests in a set of wells, possibly at different depths between packers, and to record the hydraulic head response at the remaining wells. The response is then analyzed in terms of finding the transmissivity distribution on a grid (tomogram) as done typically in geophysics or medical imaging. Over the past two decades, several numerical and field studies have demonstrated how steady-state hydraulic tests can be used to image hydraulic conductivity or transmissivity (e.g., [12], [21], [40], [59]). Transient pumping tests can also be used to estimate both of the specific storage and hydraulic conductivity fields (e.g., [3], [27], [64]). Hydraulic conductivity tomography can also be combined with the temporal moments method to reconstruct the hydraulic parameters heterogeneities (e.g.: [38], [39], [62], [65]). This method has found a broad number of applications in the recent literature (e.g., [5], [48]). In steady-state conditions, such problem is similar to electrical resistivity tomography in geophysics (see recently [63]). The inverse problem is generally non-unique and the resolution depends on the number of hydraulic data and prior information.
To incorporate additional information in the inverse problem, various geostatistical frameworks have been introduced in order to use any known prior knowledge of the spatial variability of the transmissivity. Two approaches have been extensively used in hydraulic tomography: the quasi-linear geostatistical approach by [35] and the Sequential Linear Estimator (SLE) by Yeh et al. [61]. Despite the differences in their mathematical formulations, these two approaches are similar in their spirit of using spatial statistics (mean, variance, and correlation length scales) of hydraulic properties as prior information. The quasi-linear geostatistical approach of Kitanidis [35] is formally equivalent to the Tikhonov regularization family of methods. These approaches are also used in the realm of geophysical inversion with a generalized covariance functions included in a regularization term added to the model misfit term in the objective function to minimize. This approach is based on the linearization of the functional relationship between parameters (hydraulic conductivities) and measured quantities (hydraulic head) and it comprises a step of estimating the geostatistical parameters, such as the variance and correlation lengths, used for the prior covariance function [13], [43], [45].
The Sequential Linear Estimator (SLE) approach of [61] treats the model parameter as a stochastic processes characterized by their spatial statistics (mean, variance, and correlation lengths). Afterwards, the SLE uses kriging to obtain the most likely parameters values given the data. It then updates (i) the estimated model parameters using the differences between the observed heads and (ii) the covariances to obtain residual covariances. In addition, the SLE goes beyond cokriging because it assumes a non-linear relationship between parameters fields and hydraulic data [61]. The SLE approach has proven to be an efficient tool in hydraulic tomography. Its success has been demonstrated in saturated flow synthetic tests [59], [64]; [60], [58]; [51], in variably-saturated flow synthetic tests [42], in pneumatic tomography [44], in sandbox experiments [3], [32], [40], [41], and in field experiments [4], [6], [30], [31], [49].
Another approach is the pilot point method, which is applied as a parameterization technique for reducing the number of unknown parameters to estimate by evaluating the model parameters only at few key-locations in the aquifer (called pilot points) where the values of the transmissivity can be assessed to fit the hydraulic head measurements via deterministic or stochastic optimization processes [16], [36], [37], [55]. The full spatial field is obtained by geostatistical interpolation (kriging) of the values of the transmissivity at the pilot points by assuming a cross-covariance function. However, one limitation of the pilot point method is due to the fact that the covariance function used for the interpolation of the model parameter values from those obtained at the pilot points is fixed [30]. This means that the pilot point method does not optimally exploit the information provided by the hydraulic head data in the sense that it does not fully take advantage of the changes in the flow field due to changes of the pumping locations.
The previous approaches provide some smoothness in the estimated hydrological parameters because they estimate the conditional means of such parameters. However, these smooth parameters successfully predict and reproduce the best unbiased mean data. These approaches are all adequate to characterize smoothly heterogeneous aquifers. However the distribution of the hydraulic parameters in the subsurface is not necessarily smooth. The subsurface is generally characterized by a set of facies, each of them with their degree of heterogeneities. Inside each facies, the distribution of the material properties may be described with classical geostatistical indicators such as semi-variograms. That said, they may exist drops in the distribution of the transmissivity at the boundaries between the facies. Note that under some particular conditions (especially fractured media), one can imagine that the drop in the transmissivity is not necessarily determined by the geology alone and tectonic aspects need to be accounted for as well. Several approaches have been suggested to overcome the existence of discontinuities in the distribution of the hydraulic properties. Hu and Chugunova [29] provided for instance a review of the multiple-point geostatistical approach. They studied the efficiency of using multiple-point simulations on multiple grids to recover the patterns observed in training images. Jafarpour et al. [34] developed a new regularization approach, which has proven to be efficient for the reconstruction of geologic features, facies, and porosity distributions. Cardiff and Kitanidis [11] used the level set approach to identify the shape and location of the boundaries between separate hydrofacies in some simple geometries.
The strategy we follow in this paper is to use additional information such as those provided by seismic or GPR cross-sections, or geological expertise in the hydraulic tomography problem. At the opposite of the methods derived by Brauchler et al. [1], [10], our approach does not require a statistical relationship between seismic/electromagnetic properties (e.g., P-wave velocity) and hydraulic properties. It only requires that GPR or seismic reflectors are associated with hydraulic conductivity drops that are not necessarily correlated in magnitude but are collocated with hydraulic conductivity or transmissivity drops. To take advantage of this colocation, we use an approach derived by Hale [24], who introduced two new techniques called image-guided inversion and interpolation methods. The first approach is to extract and use structural features from an image to guide the inversion of geophysical data. The second approach is similar to focusing techniques in image processing starting with the fundamental question: how can we refocus a blurred image? There are several new ideas behind these methods. First, in contrast to classical interpolation methods requiring the definition of Euclidian distances from locations where data are unknown to locations where data are known, the image-guided interpolation uses orientations and shapes of features in images to predict unknown data values. This approach has proven to be very useful in geophysical tomography to obtain inverted models that honor geological structures and give more accurate results in terms of the value of the material properties. For instance, Ma et al. [42] used the image-guided interpolation in waveform inversion of seismic data using as guiding image a migrated seismic section. They reduced the number of model parameters by discretizing the space parameters in a sparser space on which they apply the image-guided interpolation. Zhou et al. [63] used the image-guided inversion to invert apparent resistivity data using guiding images derived from georadar (GPR) or geological cross-sections. In addition, they applied the image-guided interpolation on their inverted model in order to focus smooth resistivity tomograms using the structural constraints.
As mentioned above, the steady-state transmissivity inverse problem is mathematically similar to the electrical resistivity inverse problem. The main difference is that usually much more data are available in electrical resistivity inversion because of the easiness in collecting resistivity data. In this paper, we adapt the methodology used recently in Zhou et al. [63] in order to perform image-guided inversion for transmissivity tomography. As far as we know, such method has never been used in hydrogeology. We solve the hydraulic inverse problem for four synthetic 2D-case studies in steady state conditions with the main goal to introduce these new techniques to simple case studies. The structural constraints will be derived from an image corresponding to the hydraulic pattern (assumed to be perfectly known in this paper for Case Studies 1, 2, and 4) and imposed via a four direction smoothing matrix. Finally, we will study the effect of incorporating erroneous hydraulic pattern on the quality of the inverted model (Case Study 3).
Section snippets
Forward problem
The groundwater flow in a saturated, heterogeneous, but locally isotropic confined aquifer under steady state regime is governed by the following partial differential equationwhere denotes the transmissivity (in ), denotes the hydraulic head (in m) and (in ) is the pumping rate. Eq. (1) is subject to the following boundary conditionswhere is the prescribed head at Dirichlet boundary . Note that other boundary conditions can be eventually employed
Numerical case studies
In this section, we explore the efficiency of the image-guided inversion for estimating the transmissivity field from the hydraulic head data of a shallow unconfined aquifer. Four synthetic cases are studied. In order to demonstrate the advantage of our guided inversion approach, we compare its result to the conventional Gauss–Newton method (using simple isotropic smoothing along the x- and y-directions) for each of the four test cases. For the first two cases and the last one, we suppose that
Conclusion
We have investigated the advantage of using the image-guided inversion to assess the transmissivity field of heterogeneous aquifers. We use the structural information contained in a guiding image (either obtained from seismic data, GPR data or from some geological expertise) to guide the inversion of the transmissivity field through image-processing. The image-guided approach offers a promising alternative to incorporate geophysical or geological information to reconstruct the hydraulic
Acknowledgments
AJ thanks “Agence de l’Eau of Haute-Normandie” (France), “La Maison de l’Estuaire”, and GPMH for funding the project “Hydromar”. AR thanks NSF for the SmartGeo Educational Program (Project IGERT: Intelligent Geosystems; DGE-0801692), the NSF PIRE project and DOE (Geothermal Technology Advancement for Rapid Development of Resources in the U.S., GEODE, Award #DEE0005513). We thank the Editor and four referees for their constructive comments.
References (65)
- et al.
Improved predictions of saturated and unsaturated zone drawdowns in a heterogeneous unconfined aquifer via transient hydraulic tomography: laboratory sandbox experiments
J Hydrol
(2012) - et al.
2-D joint structural inversion of cross-hole electrical resistance and ground penetrating radar data
J Appl Geophys
(2012) - et al.
Efficient solution of nonlinear, underdetermined inverse problems with a generalized PDE model
Comput Geosci
(2008) - et al.
Stochastic inversion of pneumatic cross-hole tests and barometric pressure fluctuations in heterogeneous unsaturated formations
Adv. Water Resour.
(2008) - et al.
Derivation of site-specific relationships between hydraulic parameters and p-wave velocities based on hydraulic and seismic tomography
Water Resour Res
(2012) Dynamics of fluids in porous media
(1988)- et al.
Capturing aquifer heterogeneity: comparison of approaches through controlled sandbox experiments
Water Resour Res
(2011) - et al.
Three-dimensional transient hydraulic tomography in a highly heterogeneous glaciofluvial aquifer–aquitard system
Water Resour Res
(2011) - et al.
Comparison of hydraulic tomography with traditional methods at a highly heterogeneous site
Groundwater
(2015) - et al.
A field assessment of the value of steady shape hydraulic tomography for characterization of aquifer heterogeneities
Water Resour Res
(2007)
The need to adapt the exploration model from the oil patch to contaminated-site characterization: a case from Hill AFB, Utah, USA
Leading Edge
A field assessment of high-resolution aquifer characterization based on hydraulic travel time and hydraulic attenuation tomography
Water Resour Res
Bayesian inversion for facies detection: an extensible level set framework
Water Resour Res
3-D transient hydraulic tomography in unconfined aquifers with fast drainage response
Water Resour Res
A field proof-of-concept of aquifer imaging using 3-D transient hydraulic tomography with modular, temporarily-emplaced equipment
Water Resour Res
Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model
Geostat Nat Resour Charact Part
Constructing piecewise-constant models in multidimensional minimum-structure inversions
Geophysics
Fast structural interpretation with structure-oriented filtering
Geophysics
Generalized cross-validation for large-scale problems
J Comput Graphical Stat
Identification of the permeability distribution in soil by hydraulic tomography
Inverse Prob
Three-dimensional modelling and inversion of DC resistivity data incorporating topography—II. Inversion
Geophys J Int
Analysis of discrete ill-posed problems by means of the L-curve
SIAM Rev
Hydraulic tomography for detecting fracture zone connectivity
Ground Water
Estimation of effective hydrogeological parameters in heterogeneous and anisotropic aquifers
J Hydrol
Multiple-point geostatistics for modeling subsurface heterogeneity: A comprehensive review
Water Resour. Res.
Robustness of joint interpretation of sequential pumping tests: numerical and field experiments
Water Resour Res
Hydraulic tomography in fractured granite: Mizunami under- ground research site, Japan
Water Resour Res
Cited by (25)
Image-guided structure-constrained inversion of electrical resistivity data for improving contaminations characterization
2024, Computers and GeosciencesReconstruction of hydrogeological parameter distributions by exploiting structural similarities
2023, Advances in Water ResourcesCharacterization of the non-Gaussian hydraulic conductivity field via deep learning-based inversion of hydraulic-head and self-potential data
2022, Journal of HydrologyCitation Excerpt :This is because the information on K heterogeneity contained in hydraulic heads and SP signals is complementary (Soueid Ahmed et al., 2016b; Kang et al., 2020). To tackle the non-Gaussian K inversion problem, previous studies have adapted the traditional geostatistical approaches by utilizing the knowledge of geological structure for parameter zonation (e.g., Soueid Ahmed et al., 2015; Tso et al., 2016; Zhao et al., 2016; Zhao and Illman, 2017; Zha et al., 2017). Consequently, the estimation accuracy of the adapted geostatistical approaches relies primarily on the quantity and quality of prior geological information.
Mapping of hydraulic transmissivity field from inversion of tracer test data using convolutional neural networks. CNN-2T
2022, Journal of HydrologyCitation Excerpt :To overcome this issue, constraints and prior information must be included which narrow the possibilities and promote convergence to a physically meaningful model (Kitanidis, 1997; Zha et al., 2017). This optimization can be performed with an algorithm that belongs to one of three types: deterministic, stochastic, and global (Carrera et al., 2005; Reuschen et al., 2020; Sen et al., 1995; Soueid Ahmed et al., 2015; Xu and Gómez-Hernández, 2018). In deterministic algorithms, such as Gauss-Newton and Conjugate Gradient, optimization aims at a local minimum by iteratively computing the gradient of the objective function.