An iterative particle filter approach for coupled hydro-geophysical inversion of a controlled infiltration experiment

Abstract

The modeling of unsaturated groundwater flow is affected by a high degree of uncertainty related to both measurement and model errors. Geophysical methods such as Electrical Resistivity Tomography (ERT) can provide useful indirect information on the hydrological processes occurring in the vadose zone. In this paper, we propose and test an iterated particle filter method to solve the coupled hydrogeophysical inverse problem. We focus on an infiltration test monitored by time-lapse ERT and modeled using Richards equation. The goal is to identify hydrological model parameters from ERT electrical potential measurements. Traditional uncoupled inversion relies on the solution of two sequential inverse problems, the first one applied to the ERT measurements, the second one to Richards equation. This approach does not ensure an accurate quantitative description of the physical state, typically violating mass balance. To avoid one of these two inversions and incorporate in the process more physical simulation constraints, we cast the problem within the framework of a SIR (Sequential Importance Resampling) data assimilation approach that uses a Richards equation solver to model the hydrological dynamics and a forward ERT simulator combined with Archie's law to serve as measurement model. ERT observations are then used to update the state of the system as well as to estimate the model parameters and their posterior distribution. The limitations of the traditional sequential Bayesian approach are investigated and an innovative iterative approach is proposed to estimate the model parameters with high accuracy. The numerical properties of the developed algorithm are verified on both homogeneous and heterogeneous synthetic test cases based on a real-world field experiment.

Introduction

Electrical Resistivity Tomography (ERT) is a practical, cost-effective, indirect tool for collecting soil and moisture content data in subsurface environments [1], [2], [3], [4], [5]. When applied to the simulation of the dynamics of the vadose zone, ERT relies on the inversion of the Direct Current (DC) flow equation providing an image of the electrical resistivity [4], with the soil moisture pattern reconstructed from petrophysical relations, such as, e.g., Archie's Law [6]. A second inverse problem is finally used to estimate hydrological model parameters. It is well known that inverse modeling of a parabolic diffusion equation is generally an ill-posed problem and regularization techniques are often employed to achieve well-posedness [2], [7], [8], [9]. Traditional geophysical inversion is at the same time an over- and under-constrained problem, in the sense that the problem character can change in space, and benefits from the use of prior information embedded in the regularization procedure [10]. However, imposing smoothness via regularization may introduce inaccuracies or even unphysical constraints into the estimates of the hydrological properties [11]. ERT has been widely used to monitor vadose zone processes (e.g. [1], [12], [13]) but it is well known that the inversion procedure can produce mass balance errors [14] especially when surface ERT is used to monitor water infiltration into soil [15], [5], [16] due a rapid decrease of ERT resolution with depth. To cope with this limitation coupled hydro-geophysical approaches seem highly promising [17]. By these procedures, the spatial distribution and the temporal dynamics of the geophysical properties are enforced by a physically based hydrologic model combined with petrophysical relations, and explicit assumptions for spatial and temporal regularization are no longer needed.

Even though the coupled approach avoids an independent geophysical inversion, estimation of the hydrologic properties (e.g. soil hydraulic parameters) is still a highly non-linear, mixed-determined inversion problem. For these reasons, although parameter estimation can be made theoretically well-posed, the physical interpretation of the estimated parameters is still not well understood [18]. The presence of structural model errors (model approximations, uncertain initial conditions, etc.), as well as measurement uncertainties, suggests that a deterministic search for the best parameters is not likely to converge to a single set of “true” values. A stochastic approach based on ensemble forecasting seems therefore the most appropriate solution procedure [18], [19].

Sequential Data Assimilation (S-DA) methods (typically called filters) have been successfully applied to improve model predictions by incorporating real system observations onto the dynamical model and have been already employed to correct the hydrological states of groundwater infiltration models [20]. Their ability to include structural and parametric error distributions make them particularly attractive for application to the problem of dynamic parameter estimation [18]. Because of the high nonlinearity of porous media infiltration models, the typical filtering method used in hydrological applications is the Ensemble Kalman Filter (EnKF) [21]. Notwithstanding the linear optimality properties of the Kalman Gain [22], the main limitation of EnKF is that it is based on the Gaussian approximation of the filtering probability distribution, possibly leading to inaccurate results or even divergence of the posterior pdfs in presence of a strongly nonlinear relation between observations and state variables [23], [24], [25]. To cope with arbitrary non-Gaussian prior distributions, the family of particle filters is a highly attractive alternative, as it is directly based on the Bayesian filtering rule [26], [27]. Particle filters have been recently introduced into hydrology [28], [29], [30], [31], [25] and used also for estimation of hydrological model parameters [32], [33], [34]. All these latter studies focus on the assimilation of direct hydrological information (e.g. discharge [25] or soil moisture data [35], [36], [37]). A coupled hydro-geophysical parameter estimation procedure by S-DA has been presented by [38], but its ability to provide accurate estimates of unknown model parameters remains to be proven, as shown by the consistent underestimation of saturated hydraulic conductivity in the results of [38]. As a matter of fact, the structural uncertainties of both the hydrologic evolution and geophysical observation models strongly affect the estimated parameters. Sequential filters correct both model parameters and state variables at each assimilation time, yielding identified parameter values that vary in time [18]. Compared to smoothers or other more sophisticated inversion methods (e.g., Markov Chain Monte Carlo methods [39], [40]) the filtering approach is computationally more efficient when dealing with a detailed and spatially resolved simulation model such as the coupled Richards equation-ERT solver here employed.

In this paper we propose an iterative procedure to overcome the problem of the sensitivity to the initial guess and provide accurate identification of unknown model parameters from indirect state information. The method is grounded on a Sequential Importance Resampling (SIR) particle filter, already tested in similar hydrological applications [25], [38], whereby an ERT forward simulation model is embedded into the observation equation and both parameter and state distributions are updated at each assimilation step. Iteration is introduced by sequentially repeating until convergence the same simulation period, using as initial guess the state values and parameter pdfs evaluated from the results of the previous iteration. Compared to more sophisticated statistical updates, the use of iterations allows the inclusion of a less accurate but computationally more efficient inversion scheme able to cope with large dimensional problems.

We validate the methodology on synthetic test cases and apply the methods to a field experiment comparing the results of our procedure with traditional uncoupled inversion of ERT data. We focus on both homogeneous and heterogeneous systems with parameters distributed by zones. The proposed procedure displays convergence of the posterior distribution towards the correct value of the hydraulic conductivity in both the homogeneous and heterogeneous scenarios independently from the initial guess. The numerical results obtained from the synthetic test cases show that the iterative approach yields faster convergence with respect to standard DA methods, using consistently smaller ensemble sizes and a drastic reduction of the number of forward model runs, in particular for the heterogeneous test case. The results obtained in the application to the real world problem are consistent with the desired physical constraints at relatively low computational costs, thus improving significantly on existing coupled flow-ERT procedures.