Coupled inverse modeling of a controlled irrigation experiment using multiple hydro-geophysical data

Highlights

• An irrigation experiment is monitored using time-lapse multi-offset GPR and ERT.

• GPR guided wave data are inverted to produce depth of the infiltration front.

• We implemented a Data Assimilation procedure based on Richards’ equation modelling.

• Both ERT and GPR data are assimilated to give a best estimate of hydraulic parameters.

Abstract

Geophysical surveys can provide useful, albeit indirect, information on vadose zone processes. However, the ability to provide a quantitative description of the subsurface hydrological phenomena requires to fully integrate geophysical data into hydrological modeling. Here, we describe a controlled infiltration experiment that was monitored using both electrical resistivity tomography (ERT) and ground-penetrating radar (GPR). The experimental site has a simple, well-characterized subsoil structure: the vadose zone is composed of aeolic sand with largely homogeneous and isotropic properties. In order to estimate the unknown soil hydraulic conductivity, we apply a data assimilation technique based on a sequential importance resampling (SIR) approach. The SIR approach allows a simple assimilation of either or both geophysical datasets taking into account the associated measurement uncertainties. We demonstrate that, compared to a simpler, uncoupled hydro-geophysical approach, the coupled data assimilation process provides a more reliable parameter estimation and better reproduces the evolution of the infiltrating water plume. The coupled procedure is indeed much superior to the uncoupled approach that suffers from the artifacts of the geophysical inversion step and produces severe mass balance errors. The combined assimilation of GPR and ERT data is then investigated, highlighting strengths and weaknesses of the two datasets. In the case at hand GPR energy propagates in form of a guided wave that, over time, shows different energy distribution between propagation modes as a consequence of the evolving thickness of the wet layer. We found that the GPR inversion procedure may produce estimates on the depth of the infiltrating front that are not as informative as the ERT dataset.

Introduction

Hydrological research increasingly requires detailed information to feed data-hungry numerical models. For this reason, geophysical data are increasingly called into play to fill the lack of spatial and sometimes temporal resolution of traditional hydrological data. This is particularly true for the vadose zone, where the difficulties for obtaining direct measurements, the general lack of knowledge and the uncertainty on the soil parameters and their spatial heterogeneity often lead to develop numerical models that cannot reproduce the behavior of the real systems, unless they are strongly constrained by multiple, extensive and complementary data.

The vadose/unsaturated zone is home to a number of complex key processes that control the mass and energy exchanges in the subsurface (soil water migration) and between the subsurface and the atmosphere (rain infiltration, soil evaporation and plant transpiration). The understanding of vadose zone fluid-dynamics is key to the comprehension of a large number of hydrologically-controlled environmental problems, with strong implications in water resources management and subsurface contaminant hydrology. Unsaturated processes are also key factors in a number of important issues, such as the availability of water for agriculture, slope stability, and floods. The dependence of the hydro-geophysical response on changes in soil moisture content is the key mechanism that allows the monitoring of the vadose zone in time-lapse mode via non-invasive techniques. The use of these techniques can provide high-resolution images of hydro-geological structures in the shallow and deep vadose zones and, in some cases, a detailed assessment of dynamical processes in the subsurface.

The estimation of the time and space variations of water content using non-invasive methodologies has been the focus of intensive research over the past three decades. Among the numerous techniques developed in literature for such a goal, such as electromagnetic induction, off-ground ground-penetrating radar, surface nuclear magnetic resonance, in this work we consider electrical resistivity tomography (ERT) and ground-penetrating radar (GPR). These techniques measure the electrical resistivity ρ (Ωm) and the relative dielectric permittivity ${\epsilon }_r$ (–) of the porous media, respectively. For both methods the determination of soil water content is based upon existing relationships that link water content to the geophysical quantities measured (e.g., [1], [7], [8], [50], [58]).

When used to study hydrological dynamics, GPR surveys are often performed to detect changes in soil moisture content via the variation of dielectric permittivity, generally measured from GPR travel times in a variety of configurations (e.g., [15], [16], [32]), such as borehole-to-borehole (e.g., [49], [51], [52]) or borehole-to-surface (e.g., [64]). However, the most common setup uses GPR antennas from the ground surface, even though only few studies with this configuration have been focused on the understanding of the dynamics of the water front during irrigation (e.g., [26], [36], [40], [45]) or using natural rainfall [10]. When working solely from the ground surface, three approaches are possible to determine soil moisture content: (a) use the velocity of the direct ground wave, (b) estimating velocity from the reflected events, (c) estimating impedance and thus velocity from the reflected GPR signal. Approaches (a) and (b) share in fact the same operational characteristics, needing the two antennas to be separated from each other. Approach (c) does not require antenna separation and exploits the physics of the reflection mechanism, with its own advantages and disadvantages (e.g., [37], [53]), and with more limited applications so far. When the two antennas are separated from each other, the survey can be conducted in wide angle reflection and refraction (WARR) mode (e.g., [63]), where one antenna is kept fixed while the other is moved, or common midpoint (CMP) [25], [28], [56], where both antennas are moved simultaneously to keep the same mid-point. Both sounding techniques allow for a good identification of direct waves through the air and the ground. These methods are also employed for the estimation of velocity from the reflected events, even though for this use the normal move-out approach, typical of seismic processing, may not be ideal (see [3] for a discussion). The estimation of velocity from the direct wave through the ground is the most widely adopted approach for vadose zone applications (e.g. [30], [31], [63]). However, in some cases direct arrivals are not so straightforward to identify and can be confused with other events. This can happen in the presence of critically refracted radar waves [6] or guided waves [2], [57], [61]. A water front that infiltrates from the surface can give rise to such ambiguous situations, as the wet and consequently low velocity layer, lying on top of a faster (drier) media, can give rise to critically refracted waves [6] as well as act as a waveguide confined between two faster layers: the air above and the drier media below [57], the two situations being defined by the ratio between the wavelength and the layer thickness. Therefore, to study infiltrating fronts, maximum care must be given in understanding the nature of the observed, multi-offset GPR signal, possibly exploiting the entire information content of the data (e.g. [9]).

ERT measurements [5] have been widely employed to monitor water dynamics, as variations of moisture content [4], [22] and salinity of pore water [47] leads to changes in the electrical properties of the media [17], [35]. However, it is well known that resolution limitations [23] can produce severe mass balance errors [54] even in the most favorable cross-hole configurations. The problem is even more serious when only surface ERT are used to monitor natural or artificial irrigation from the ground surface [13], [20], [21], [43], [60] where resolution dramatically drops with depth and a direct conversion of inverted resistivity values into estimates of soil moisture content may prove elusive.

Geophysical measurements can be informative of the hydrological response of the soil and subsoil if applied in time-lapse monitoring mode: some geophysical quantities (in this case, ρ and ${\epsilon }_r$) are useful indicators of changes in the hydrological state variables, such as moisture content or pore water salinity. However, in order to extract this hydrological information, the assimilation of measurements in a hydrological model is needed. Two different approaches may be applied, named respectively “uncoupled” and “coupled” hydro-geophysical inversions [24], [29]. The procedure for an uncoupled inversion can be summarized by the following steps:

• The spatial distribution of the geophysical quantity of interest (e.g. electrical resistivity for ERT) is derived from the inversion of geophysical field data.

• The application of a petro-physical relationship leads to obtaining, from the geophysical quantity, an estimation of moisture content distribution.

• The estimated hydrologic state variable, in its spatio-temporal distribution, is used to calibrate and constrain a hydrological model, thus identifying the corresponding governing parameters.

The inversion of geophysical measurements is usually an ill-posed inversion problem that can be tackled introducing prior information. If no solid independent information is available, the most common approach is the introduction of a regularizing functional, commonly a smoothness constraint [42]. As a consequence of ill-posedness and regularization, the inversion procedure can lead to artifacts, misinterpretations and unphysical results, especially in the subsurface regions where the sensitivity of the measurements is low (consider e.g. [23]). To overcome these problems, a coupled hydro-geophysical modeling can be applied:

• A hydrological model is used to predict the evolution of hydrological state variables – e.g. moisture content – on the basis of a set of hydrological governing parameters, the identification of which is the final aim of the inversion.

• A suitable petrophysical relationship (same as for point (2) above) translates hydrological state variables into geophysical quantities, such as resistivity or dielectric permittivity.

• The simulated geophysical quantities are used to predict the geophysical field measurements.

• A comparison between predicted and measured geophysical field measurements allows a calibration of the complex of hydrological and geophysical models (thus the name “coupled inversion”), leading to the identification of the hydrological parameters, that is the key objective of the study.

In this work we follow a coupled approach within the framework of data assimilation (DA). DA schemes are mathematical tools of common use in hydrological applications. The main idea behind DA is using the field measurements to correct numerical simulations obtained with a hydrological model, thus modifying their governing parameters. This is possible by the recursion of forecast steps, which simulate the time-evolution of the probability density function (pdf) of the hydrological process, and analysis (or update) steps, which compute a posterior pdfs of the model parameters and state variables by assimilating the measurements (e.g., [41], [44]). A few examples of coupled hydro-geophysical inversion exist in the literature (e.g., [10]) but the use of DA techniques is less widespread [48], [59].

The present work focuses on a field experiment where artificial irrigation is monitored in time-lapse mode from the surface via both ERT and GPR. The homogeneous nature of the site, made of aeolic sand deposits, provides a simplified case study suitable to evaluate the performance of coupled hydro-geophysical inversion and test the information content of different geophysical data. Both GPR and ERT geophysical measurements are assimilated into the hydrological model CATHY [11], that is employed for the numerical simulation of the experiment. We elected to use the iterative sequential importance resampling (SIR) proposed by Manoli et al. [39] as a DA technique to estimate the model saturated hydraulic conductivity. This technique is particularly designed to assimilate geophysical measurements in a coupled hydro-geophysical model: the geophysical measurements are blended in the simulation to update the state of the system, estimate the model parameters and quantify the model uncertainties.

The specific goals of this work are:

• To analyze in detail the nature of the WARR GPR data collected during the irrigation experiment, verifying whether or not complex refraction and waveguide phenomena occur during the progression of the wetting front, and how and to what extent this type of data can be processed and inverted.

• To assess the effectiveness of incorporating ERT and GPR data in a coupled hydro-geophysical inversion procedure that, using the unsaturated flow equations, point directly at the estimation of the saturated hydraulic conductivity, and to compare this approach with the results of a classical uncoupled inversion approach.

• To evaluate to what extent the information that can be obtained from GPR and ERT data corroborate each other, how the independent assimilation of each data type performs, if the assimilation of both geophysical techniques adds information with respect to separate procedures, and finally what is the value of using both techniques to monitor the infiltration process.

The paper is organized as follows: Section 2 is dedicated to the description of the hydrological model and the DA procedure used for the coupled inversion of the geophysical data. After presenting the hydrological experiment taken into consideration (Section 3), in Sections 4 GPR data analysis, 5 ERT data analysis we analyze the GPR and ERT data, respectively. In Section 6 we describe the setup for the DA procedure in this experiment. The benefits of the coupled inversion are presented in Section 7. The major conclusions of this work are summarized in Section 8.