Image-guided inversion in steady-state hydraulic tomography

Highlights

• Image-guided inversion and interpolation applied to hydraulic tomography.

• Hydrofacies boundaries provided by georadar or some prior geological expertise.

• The inverted hydraulic conductivity distributions are close to the true ones.

• The prior textural information is enforced in the model covariance matrix.

Abstract

In steady-state hydraulic tomography, the head data recorded during a series of pumping or/and injection tests can be inverted to determine the transmissivity distributions of an aquifer. This inverse problem is usually under-determined and ill-posed. We propose to use structural information inferred from a guiding image to constrain the inversion process. The guiding image can be drawn from soft data sets such as seismic and ground penetrating radar sections or from geological cross-sections inferred from the wells and some geological expertise. The structural information is extracted from the guiding image through some digital image analysis techniques. Then, it is introduced into the inversion process of the head data as a weighted four direction smoothing matrix used in the regularizer. Such smoothing matrix allows applying the smoothing along the structural features. This helps preserving eventual drops in the hydraulic properties. In addition, we apply a procedure called image-guided interpolation. This technique starts with the tomogram obtained from the image-guided inversion and focus this tomogram. These new approaches are applied on four synthetic toy problems. The hydraulic distributions estimated from the image-guided inversion are closer to the true transmissivity model and have higher resolution than those computed from a classical Gauss–Newton method with uniform isotropic smoothing.

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).