Efficient geostatistical inversion of transient groundwater flow using preconditioned nonlinear conjugate gradients

Highlights

• We present a preconditioned conjugate gradient method for geostatistical inversion.

• The prior covariance matrix is used as preconditioner at negative computational cost.

• The approach incorporates linearized uncertainty quantification.

• The method is particularly efficient for the inversion of transient data.

• Transient inversion may speed up experiments and improve quality of inversion results.

Abstract

In the geostatistical inverse problem of subsurface hydrology, continuous hydraulic parameter fields, in most cases hydraulic conductivity, are estimated from measurements of dependent variables, such as hydraulic heads, under the assumption that the parameter fields are autocorrelated random space functions. Upon discretization, the continuous fields become large parameter vectors with $\mathcal{O}({10}^4-{10}^7)$ elements. While cokriging-like inversion methods have been shown to be efficient for highly resolved parameter fields when the number of measurements is small, they require the calculation of the sensitivity of each measurement with respect to all parameters, which may become prohibitive with large sets of measured data such as those arising from transient groundwater flow.

We present a Preconditioned Conjugate Gradient method for the geostatistical inverse problem, in which a single adjoint equation needs to be solved to obtain the gradient of the objective function. Using the autocovariance matrix of the parameters as preconditioning matrix, expensive multiplications with its inverse can be avoided, and the number of iterations is significantly reduced. We use a randomized spectral decomposition of the posterior covariance matrix of the parameters to perform a linearized uncertainty quantification of the parameter estimate.

The feasibility of the method is tested by virtual examples of head observations in steady-state and transient groundwater flow. These synthetic tests demonstrate that transient data can reduce both parameter uncertainty and time spent conducting experiments, while the presented methods are able to handle the resulting large number of measurements.

Introduction

The spatially distributed hydraulic conductivity K(x) $\left[L{T}^{-1}\right]$ is widely considered as the most important parameter field for groundwater flow and, indirectly, solute transport. Estimating K(x) is a prerequisite for groundwater management, assessment of solute transport in aquifers, and remediation of contaminated sites. Because direct measurement of hydraulic conductivity can only be done in the laboratory using samples of aquifer material, the estimation of K(x) in the field mostly relies on hydraulic-head measurements, either under conditions of ambient groundwater flow, or during hydraulic tests (Liu, Yeh, Gardiner, 2002, Liu, Illman, Craig, Zhu, Yeh, 2007, Straface, Yeh, Zhu, Troisi, Lee, 2007, Zhu, Yeh, 2005). Additional information about hydraulic conductivity can be obtained from concentration data, e.g., from tracer tests, due to the dependence of the velocity on K(x) (Datta-Gupta, Yoon, Vasco, Pope, 2002, Vasco, Datta-Gupta, 1999, Yeh, Zhu, 2007). Traditional methods of estimating K(x) rely on analytical solutions of the groundwater-flow and potentially solute-transport equations, typically assuming uniform parameter values (e.g., Kruseman, Ridder, et al., 1990, Rumynin, 2011). In hydrogeological and engineering practice, numerical models are commonly calibrated to fit measured hydraulic heads and other observations, conceptualizing the subsurface as a set of space-filling, simple-shaped geometric bodies, such as layers and blocks, with uniform hydraulic properties (e.g., Hill and Tiedemann, 2007). Identification of the hydraulically relevant geobodies, however, remains a challenge, and the hydraulic variability within the geological units may be larger than between them.

In the last three decades, geostatistical methods have become widely accepted in order to characterize the spatial variability of hydraulic subsurface properties across many scales, notwithstanding the identification of dominant geological features (e.g., Rubin, 2003; Kitanidis, 1997; Zhang, 2002). The key assumption is that the hydraulic-conductivity field is the outcome of a spatially autocorrelated random process that is considered stationary unless data strongly supports spatial variation of statistical properties. The geostatistical premise can also be used as prior information in the estimation of hydraulic conductivity as a continuous spatial field from direct data, i.e., conductivity values themselves, and indirect data such as hydraulic-head measurements. Upon discretization, one hydraulic-conductivity value is assigned to each node or element of a numerical grid, leading to $\mathcal{O}\left({10}^4\right)$ parameter values in 2-D, and $\mathcal{O}({10}^6-{10}^7)$ in 3-D applications. In this context, the geostatistical prior information acts as a regularization term, because the number of measurements typically is by orders of magnitude smaller than the number of grid cells, resulting in an ill-posed inverse problem.

Most geostatistical inverse methods assume that the logarithm of hydraulic conductivity Y(x) is an – at least blockwise – second-order stationary multi-Gaussian field. Some approaches replace the measurements of indirect data, such as hydraulic heads, by virtual measurements of log-hydraulic conductivity itself at so-called pilot points (Doherty, 2003, RamaRao, LaVenue, De Marsily, Marietta, 1995), master points (Gómez-Hernánez et al., 1997), or anchors (Rubin et al., 2010). Then the log-hydraulic conductivity field between these points is estimated by interpolation or conditional simulation, and the inversion becomes estimating the virtual log-conductivity measurements. This approach effectively reduces the number of parameters to be estimated, but can lead to bias in the conditional uncertainty analysis, because the spatial patterns of the conditional Y(x)-field and its uncertainty differ whether hydraulic-head or log-conductivity measurements are used for conditioning Kitanidis (1998).

As alternatives to reduce the dimension of the parameter space, the eigen-decomposition-based Karhunen–Loève expansion (e.g., Li and Cirpka, 2006; Li and Zhang, 2007) and truncated singular-value decomposition methods (Kitanidis, Lee, 2014, Lee, Kitanidis, 2014) have been used in geostatistical inversion. These techniques act as low-pass filters. Deciding upon the right number of terms requires balancing the computational effort and the accuracy of the inversion.

An extended group of mathematically similar geostatistical inverse methods rely on variants of the cokriging equations, which can be derived from the Gauss–Newton method applied to the geostatistical inverse problem (Hoeksema, Kitanidis, 1984, Kitanidis, 1995, Rubin, Dagan, 1987, Valstar, McLaughlin, Te Stroet, van Geer, 2004, Yeh, Gutjahr, Jin, 1995, Yeh, Jin, Hanna, 1996). While in the original Gauss–Newton method, the order of the system of equations to be solved is the number of parameters (McLaughlin and Townley, 1996), it is reduced to the number of measurements plus the number of deterministic drift coefficients in the cokriging equations. A benefit of the approach is that it directly yields a linearized uncertainty estimate. A disadvantage is that it requires the calculation of the full sensitivity matrix, that is the sensitivity of each measurement with respect to each parameter. This can be achieved by analytical linearization about the prior mean (Hoeksema, Kitanidis, 1984, Rubin, Dagan, 1987), thus neglecting the non-linear relationship between hydraulic heads and log-conductivity in the inversion, direct numerical differentiation, or application of adjoint-state methods (Sun and Yeh, 1990). The latter two approaches require solving equations of the same complexity as the forward problem as often as either the number of parameters or the number of measurements. That is, the non-linear cokriging-type inversions become computational prohibitive when the log-conductivity field is highly resolved and many measurements are to be inverted.

A large number of measurements is typically achieved in monitoring of transient state variables, such as hydraulic heads during transient pumping tests or concentrations during tracer tests, which may pose a problem to cokriging-like geostatistcal inversion as discussed above. Some authors recommend picking a selected small number of observations from the continuous data streams (Cardiff, Barrash, 2011, Cardiff, Barrash, Kitanidis, 2012, Liu, Illman, Craig, Zhu, Yeh, 2007), whereas others condensate the information to a few temporal moments (Cirpka, Kitanidis, 2000, Leube, Nowak, Schneider, 2012, Li, Nowak, Cirpka, 2005, Pollock, Cirpka, 2010, Zhu, Yeh, 2006). Both approaches are associated with loosing information.

Recently, Ensemble-Kalman methods have become popular in geostatistical inversion (Chen, Zhang, 2006, Franssen, Kinzelbach, 2008, Franssen, Kinzelbach, 2009, Nowak, 2009, Schöniger, Nowak, Hendricks Franssen, 2012). In these methods, an ensemble of realizations is used to compute the covariance matrices involving the dependent variables. By this, the computation of the sensitivity matrix is circumvented. However, to avoid bias in the uncertainty estimates, denoted filter inbreeding, either a large number of realizations are needed or techniques to inflate the conditional covariance have to be applied after each conditioning step. As an alternative to performing Monte–Carlo simulations, second-order moment-generating equations have been suggested to compute the underlying covariance matrices, but here the number of partial differential equations to be solved equals the number of grid elements (Panzeri, Riva, Guadagnini, Neuman, 2013, Panzeri, Riva, Guadagnini, Neuman, 2014).

In summary, the existing geostatistical inverse methods either rely on parameterizing the continuous fields with a few terms, require the evaluation of as many adjoint problems as the number of measurements, condensate the data, or demand the evaluation of large ensembles of parameter fields. Thus, despite progress made in inversion methods in the last three decades, a choice between accuracy and computational efficiency has to be made.

The geostatistical inverse problem for transient data sets is conceptually similar to large-scale data-assimilation problems, such as those occurring in weather forecast (Caya, Sun, Snyder, 2005, Marécal, Mahfouf, 2002, Županski, Mesinger, 1995) and oceanography (Weaver et al., 2003). Here, very large vectors of state variables (with up to $\mathcal{O}\left({10}^8\right)$ elements) are updated using smaller, but still large vectors of observations (with up to $\mathcal{O}\left({10}^6\right)$ elements) that depend on the state variables in nonlinear ways. In the variational approaches toward data assimilation, an objective function is minimized that resembles that of the geostatistical inverse problem discussed here. The sheer size of the state and observation vectors prohibits methods involving the full sensitivity matrix. Conjugate Gradient methods, by contrast, require only calculating the sensitivity of the single-valued objective function with respect to the vector of values to be estimated. The latter task can efficiently be done by adjoint-state methods. The method can be expedited by choosing an appropriate preconditioning matrix.

Conjugate Gradient methods have already been used in the inversion of groundwater data, among others, in the pilot-point/master-point approaches (Gómez-Hernánez, Sahuquillo, Capilla, 1997, RamaRao, LaVenue, De Marsily, Marietta, 1995); but so far they have not been applied to the full geostatistical problem without additional parameterization. A potential reason for that could be that evaluating the gradient of the prior term requires solving a large system of equations involving the dense prior covariance matrix on the left-hand side. In the present paper, we will show that the latter can be avoided by using the prior covariance matrix for preconditioning of the Conjugate Gradient method.

The remaining article is structured as follows: Section 2 describes the inverse problem that is to be solved, stating the forward model, its governing equations and the derivation of an objective function for the inverse problem through Bayesian statistics. Section 3 discusses the solution of the arising optimization problem, describes the Gauss–Newton algorithm as a reference scheme that is often applied and introduces the Conjugate Gradient method preconditioned with the prior covariance matrix as an alternative. Section 4 lists the software packages used for the numerical implementation and specifies implementation details. Section 5 shows the application of the discussed methods in a synthetic test case, comparing the schemes regarding their computational cost and accuracy. Finally, Section 6 introduces an eigen-decomposition-based algorithm for linearized uncertainty quantification and applies it to the aforementioned synthetic test case.