Efficient estimation of hydraulic conductivity heterogeneity with non-redundant measurement information

Solving the inverse problem of identifying groundwater model parameters with measurements is a computationally intensive task. Although model reduction methods provide computational relief, the performance of many inversion methods depends on the amount of often highly correlated measurements. We propose a measurement reduction method that only incorporates essential measurement information in the inversion process. The method decomposes the covariance matrix of the model output and projects both measurements and model response on the eigenvector space corresponding to the largest eigenvalues. We combine this measurement reduction technique with two inversion methods, the Iterated Extended Kalman Filter (IEKF) and the Sequential Monte Carlo (SMC) methods. The IEKF method linearizes the relationship between measurements and parameters, and the cost of the required gradient calculation increases with increase of the number of measurements. SMC is a Bayesian updating approach that samples the posterior distribution through sequentially sampling a set of intermediate measures and the number of sampling steps increases with increase of the information content. We propose modified versions of both algorithms that identify the underlying eigenspace and incorporate the reduced information content in the inversion process. The performance of the modified IEKF and SMC methods with measurement reduction is tested on a numerical example that illustrates the computational benefit of the proposed approach as compared to the standard IEKF and SMC methods with full measurement sets.


Introduction
Aquifer parameters governing groundwater flow, such as the hydraulic conductivity, are highly variable and heterogeneous. Knowledge of the spatial parameter distribution is essential for designing successful remediation systems and developing reliable water supply resources. However, directly measuring these parameters at local scale is practically impossible. Therefore, aquifer parameters are indirectly estimated using hydraulic head measurements under different hydraulic stress conditions, such as pumping tests. The measurements are related to the sought parameters through the numerical solution of a mathematical model. The task of determining the model parameters through comparing the model response with measurements of the underlying physical system requires the solution of an inverse problem.
Groundwater inversion involves high-dimensional parameter spaces and sparse or noisy measurement data, and hence it is typically an ill-posed problem. Conventionally, the parameter space is simplified by assuming the spatial distribution of parameters to be zoned according to an interpretation of geologic information. However, such interpretation is often subjective. Recently, Hydraulic Tomography (HT) techniques have been developed to cast the model inversion problem in a probabilistic framework that treats the model parameters as random fields. The verification of HT has been demonstrated with sandbox experiments (Illman et al. 2010) and at a highly heterogeneous field site with abundance of hydraulic head measurements and validation samples (Berg and Illman 2011). One HT approach is to formulate the model inversion in a probabilistic framework through Bayesian analysis (Ghosh et al. 2007). In this context, the information contained in the measurement data is expressed by the likelihood function, whereas any knowledge on the parameter values prior to conducting the measurements is expressed by the prior distribution. The posterior distribution of the parameter is obtained through Bayes' rule as the normalized product of the prior distribution and likelihood function, whereby the prior distribution acts as a regularizer. The mathematical formulation for Bayesian inverse problems involving spatially variable parameters is given in Stuart (2010). Bayesian analysis is flexible with respect to the data types, it enables quantifying the uncertainty of the obtained parameter values and it facilitates spatial modeling through random fields. For these reasons its application is popular in HT problems (e.g. see Copty et al. 1993;Hachich and Vanmarcke 1983;Kitanidis and Vomvoris 1983;McLaughlin and Townley 1996).
The normalizing constant in Bayesian analysis consists of a high dimensional integral, which generally needs to be solved numerically. This leads to considerable computational cost, especially in cases where the model is computationally intensive. Analytical solutions are available for Gaussian random fields and linear inversion models (i.e., linear relationships between measurements and model parameters). Examples of Bayesian analysis with linear inversion models are given in Woodbury and Ulrych (2000) and Jiang and Woodbury (2006), where the moments of the random variables are treated as hyperparameters. Also, quasi-linear approaches are available, as given in Fienen et al. (2008). An extension of analytical solutions to nonlinear inversion models with Gaussian prior fields is given by the Successive Linear Estimator (SLE) method (Yeh and Zhang 1996). This approach iteratively linearizes the nonlinear relationship between hydraulic pressure head and the spatial distribution of hydraulic conductivity (Yeh and Liu 2000). Apart from hydraulic pressure head, other types of data can be incorporated, such as flux measurements (Zha et al. 2014), or geological data as prior information (Zhao et al. 2016). A method that is closely related to SLE is the Extended Kalman Filter (EKF) (Jazwinski 2007), which is applicable for mildly nonlinear problems. In contrast to the SLE, the measurement error, which is inherent in all field measurements, is an explicit part of the formulation. An improvement for nonlinear problems is the Iterated Extended Kalman Filter (IEKF) (Jazwinski 2007), which also explictly accounts for measurement uncertainty. Similar to SLE, the IEKF approach is based on a successive linearization of the model at a sequence of updated estimates.
The linearization assumption in SLE and IEKF, which renders the methods effective in problems where the hydraulic conductivity variance is small, leads to convergence issues in inversion problems with large variance. This can be circumvented through the application of sampling-based methods, which generate samples that follow the true posterior distribution of the sought parameters and use the actual forward simulation model without linearization. Markov Chain Monte Carlo (MCMC) methods are a popular group of algorithms for solving Bayesian inverse problems. A review of MCMC methods for subsurface problems is given in Yustres et al. (2012). One of the early applications of MCMC to hydraulic tomography is found in Oliver et al. (1997). Efendiev et al. (2005) suggested to accelerate MCMC by using intermediate, coarser-scaled models until the final model scale is evaluated. However, for MCMC methods, it is difficult to assess if the resulting distributions are accurate. Moreover, the sampling process may become computationally expensive as often large burn-in periods are required to reach the stationary distribution of the Markov chain. Alternative approaches are Importance Sampling (IS, Li and Lin 2015), blocking Markov Chain Monte Carlo (BMCMC, Fu and Gómez-Hernández 2009), Hybrid Nested Sampling (HNS, Elsheikh et al. 2014) and Sequential Monte Carlo (SMC, Neal 2001;Chopin 2002;Del Moral et al. 2006;Beskos et al. 2015). SMC generates samples from a sequence of distributions that gradually approach the posterior distribution. The method has found several applications in the subsurface inversion problem (e.g. Chang et al. 2012;Field et al. 2016;Pasetto et al. 2012;Montzka et al. 2011;Schöniger et al. 2012). Other sampling-based inversion methods include the Implicit sampling method (Chorin et al. 2010;Morzfeld et al. 2015), which is a special case of chainless sampling. In addition, the Ensemble Kalman Filter, which is a version of the Kalman Filter where the covariance is evaluated based on samples, has found wide use in different fields. It was first presented by Evensen (1994), and has found numerous applications and modifications (Kovachki and Stuart 2019;Gu et al. 2007;Iglesias et al. 2018).
The availability of measurements has been increasing, due to growing affordability of instrumentation and advancements in storage and transfer technologies. Although a large amount of data is beneficial for inversion problems, since their availability generally reduces the inversion uncertainty, many algorithms are likely to become less efficient with increasing observation data amount. This dilemma, most commonly solved by choosing a subset of data either randomly or subjectively based on selected measurements events, has also been tackled from different angles in a more systematic manner. In surface water model calibration, using a representative short calibration period instead of applying the full available measurement period has been suggested (Razavi and Tolson 2013). This is achieved by determining the period which has a similar mean squared error as the full model period. Another method of reducing the number of measurements with respect to time is limiting the measurements used for inversion by a sliding window, thus ignoring the influence of measurements that range further back in time than a certain threshold (Lamberti et al. 2017). However, the selection of this threshold seems to be arbitrary and could ignore long-term effects. In Mirhoseini et al. (2015), a method of sparsifying measurement correlation matrices has been introduced, in order to decrease computing effort for regression techniques, as well as to facilitate the data transfer between parallel computing units. Regression has been conducted with a subset of principal components of the predictor variables, obtained through a principal component analysis (PCA) (Huang and Antonelli 2001). Within this work, the optimal compression ratio with respect to reconstruction residuals was also investigated. The works of Chen et al. (2017) and Schillings et al. (2019) introduce Hessian information to improve the convergence of Monte Carlo-based method in inversion problems with large measurement sets. This paper proposes a method for compressing measurement data in the context of groundwater inversion. For measurements with high correlation and therefore redundancy, the proposed compression approach allows decreasing the computational effort of model inversion. The proposed approach is similar to the method in Fischer et al. (2017), but instead of applying PCA to the parameter space and reducing the number of parameters, it is applied to the observation space to express the measurements with a reduced number of fixed measurement component terms, as in Rezaie et al. (2012). In meteorological forecasting, PCA has been applied for measurement reduction in data assimilation (Collard et al. 2010). Their work primarily concentrates on the removal of measurement noise by decomposition and reconstruction of measurements using PCA. The focus of this paper is the potential in computational gain of groundwater inversion algorithms through measurement reduction. The method is combined with two Bayesian inversion algorithms, the IEKF and SMC methods, and modified versions of these algorithms are proposed. The performance of the latter is investigated with a numerical example and it is demonstrated that the proposed modifications result in a substantial decrease of the computational cost without significant loss of accuracy.
The paper is organized as follows: First, the theory behind groundwater inversion based on Bayesian Analysis is presented. Then, a method of reducing the number of measurements while keeping the contained information is introduced. The application and the benefit of this measurement reduction method for two inversion algorithms, the IEKF and SMC methods, is described next. Subsequently, we use a numerical example of a theoretical aquifer to illustrate the behavior and accuracy of the modified variants of IEKF and SMC. Finally, we discuss the observations and thoughts we obtained from the numerical studies.

Problem setting
We consider an aquifer , being a bounded subdomain of R d , d ∈ N, with hydraulic conductivity field K ∈ X , where X is a suitable space of functions defined over the spatial domain . We are interested in inferring the spatial distribution of K based on measurements of the hydraulic head at certain spatial locations from a number of pumping tests. The relationship between K and hydraulic head h j for the j-th pumping test is given by the steady-state groundwater flow equation: where Q j is the pumping/injection rate at location x Q j and δ is the Dirac function. At the model boundary ∂ , the hydraulic head is set to constant value h 1 .
Equation (1) is typically solved with a numerical method, such as the finite element (FE) method. After discretization into n elements, the FE method results in a system of linear equations that takes the following form: Here, K ∈ R n denotes the hydraulic conductivity values at the finite elements and h j ∈ R n h is the hydraulic head values at the nodal points; A(K) is the stiffness matrix that depends on the element conductivity values K, and b j is the force vector due to the j-th pumping test. The hydraulic head values at the measurement locations can be extracted as follows: where E m = [e 1 , . . . , e n m ], e i is the standard basis vector corresponding to the location of the i-th measurement and n m is the number of measurements. The vector h collects the hydraulic head values at all measurement locations and all pumping tests: with p denoting the number of pumping tests. Through Eqs.
(2)-(4), we obtain a mapping G : R n → R m from the element hydraulic conductivities K to the hydraulic heads at the measurement locations, where m = pn m denotes the total number of measurement points from all pumping tests. The goal is to determine K based on measurements d of h.

Random fields and Bayesian formulation
To approach the problem introduced in Sect. 2.1, i.e. determine K by means of measurements of h, we move into a Bayesian probabilistic framework. The Bayesian framework allows the construction of a stable solution of the inverse problem, through introducing prior information on K. This prior information acts as a regularizer to the generally ill-posed inversion problem. Additionally, the Bayesian approach enables obtaining a probabilistic estimate of K, rather than a deterministic one, which allows for quantifying the uncertainty of the estimate. The prior information is given in terms of a prior probability distribution that expresses the observer's belief on the distribution of K prior to collecting the measurements. In this context, K represents discrete instances of a spatially variable uncertain quantity, i.e. the hydraulic conductivity field K , which is modeled probabilistically by a random field. Here, we model the logarithm of K , log(K ), by a homogeneous Gaussian random field. This implies that the marginal distribution of K , i.e. its distribution at every point in space, is lognormal and hence it cannot take negative values. The field log(K ) can be completely defined by its mean and auto-covariance function, which are chosen based on past geologic investigations and information from literature. Taking log(K) to be the random variables corresponding to the midpoints of each finite element, their prior distribution will be jointly Gaussian, with mean values equal to the mean of the random field and covariance matrix evaluated based on the auto-covariance function of the random field and the coordinates of the element midpoints.
To decrease the dimensionality of the parameter space and thus further facilitate the solution of the inverse problem, the random vector log(K) can be approximated by a truncated Karhunen-Loève (KL) Expansion (Ghosal and Van der Vaart 2017), expressed as: where μ log K is the mean vector of log(K) and the pairs {λ KL,i , v KL,i } are the k largest eigenvalues and corresponding eigenvectors of the covariance matrix of log(K). The truncation order k is selected such that the vector log(K) is represented without a significant loss of accuracy, in practice through requiring that the approximation captures a significant portion of the total variance of log(K). The vector θ follows the k-dimensional independent standard Gaussian distribution N (0, I).
The parametrization of Eq. (5) implies that the identification problem is now shifted to the identification of the random variables in θ . This is done through Bayesian updating, which applies Bayes' rule to update the prior standard Gaussian density of θ , denoted f (θ ), to a posterior density: Here, L(θ) is the likelihood function that describes the measurement information and c E is a normalizing constant that ensures that the posterior density f (θ ) integrates to one: The likelihood function L(θ ) is proportional to the probability of the measurements d given a parameter state θ . The model outcome can be expressed as a function of the vector θ as h = G(θ ), where G : R k → R m is the composition of operator G and Eq. (5). The data vector d can then be expressed as a realization of G(θ ) superimposed by the outcome of an error term , i.e.
The error models the discrepancy between measurements and model outcomes due to measurement noise and/or model uncertainty. We assume that follows a zero-mean Gaussian distribution with covariance matrix C . As a consequence, the likelihood function takes the following form: where (θ) is the so-called potential (or negative log-likelihood), defined as The closer the model outputs to the measurements, the higher is the likelihood that this model input θ is close to the true value of the parameter. The covariance matrix C of the error is chosen to reflect the measurement and/or model uncertainty. For the case where models only measurement noise, it is often reasonable to assume that errors at different locations are statistically independent. In such case, the matrix C is a diagonal matrix whose diagonal elements are the variances of the errors at the measurement locations. The Bayesian framework serves as a basis for a number of inversion algorithms. These include approximation methods that attempt to approximate the posterior distribution and sampling-based approaches that generate samples from the posterior usually based on MCMC methods. The performance of most inversion algorithms deteriorates with increase of the dimension of the data d. This is related to the fact that an increase in the data dimension usually implies a higher information content. This results in a highly peaked likelihood, which hinders exploration of the high probability mass area of the posterior distribution. Next, we discuss an approach to condense the information contained in high-dimensional data vectors d. In the subsequent sections, we demonstrate how to incorporate this approach within an approximation and a sampling method, namely the IEKF and SMC methods.

Coping with correlated measurements
The measurements collected for an aquifer characterization effort are often numerous. This is either owed to the small temporal interval necessary for determining the specific storage of the aquifer or due to the small spatial interval, which can be necessary to identify medium-scale heterogeneity. However, in practice, it is computationally impossible to incorporate large numbers of measurements in the inversion process. We now discuss an approach for reducing the number of measurements, without neglecting potentially important information.
We aim at reducing the number of measurements for the updating process by determining linear combinations thereof with high information content. We do this through projecting the measurements on the space defined by the eigenvectors of the covariance matrix C hh = Cov(h, h) of the model response corresponding to the measurements h = G(θ). The matrix C hh can be decomposed as: is an m×m matrix whose columns are the eigenvectors of C hh and D = diag(λ i ) is a diagonal matrix whose entries are the corresponding eigenvalues of C hh . The matrix V is orthogonal, i.e. V T V = I where I is the m×m identity matrix. Therefore, the columns of V imply orthogonal transformations in the Euclidean space R m . We define the transformation: The components of a are uncorrelated, i.e. it holds The vector h can be expressed in the transformed coordinate system as: The sum of the variances of each component of h can be expressed as the sum of the eigenvalues, i.e.
We now sort the eigenvector matrix V and eigenvalue matrix D in order of decreasing eigenvalue. h can be approximated by keeping only the first m r ordered components in Eq. (14), which leads toĥ The total variance of the approximation in Eq. (16) reads: Through Eqs. (15) and (17), we can choose the number of eigenvectors that reflect a target percentage α of the variability in h. That is, we select m r such that We now formulate the Bayesian inversion problem in the transformed space defined by the truncated eigenvector matrix V r of dimensions m × m r . To this end, we project the data vector d on the truncated space, which gives The relationship between the reduced measurement components d r and the truncated model a r reads The modified error r is a linear transformation of a Gaussian error and is therefore itself Gaussian with zero mean and covariance matrix C r r = V T r C V r . The likelihood function expressing the reduced information is given as where the reduced potential r (θ ) reads We remark that the measurement reduction approach presented in this section requires the evaluation of the covariance matrix C hh . Exact evaluation of C hh is not possible, as h depends on the input parameters θ through the solution of a finite element system. Hence, C hh needs to be approximated. Optimally, the approximation should be based on the posterior distributions of θ and therefore h, in order to accurately reflect the variability of the measurements. The approach adopted to approximate C hh depends on the algorithm used to solve the inversion problem and includes methods such as sampling or first-order Taylor approximations. We note that, if C hh is estimated based on samples from h, the projections a coincide with the principle components of these samples and, hence, the approach is equivalent to a PCA of the model output.
In the following sections, we discuss two inversion algorithms for solving the modified inversion problem of Eq. (20) that are based on extensions of the IEKF and SMC methods.

Iterated extended Kalman filter (IEKF)
The IEKF method is an approximation method for solving Bayesian inverse problems that is based on an iterative linearization of the computational model and application of the Kalman filter on the respective linearized system. This results in successive estimates of the conditional mean and conditional covariance matrix of the parameters given the data. We introduce the IEKF for solving the Bayesian inverse problem of Eq (8). The conditional mean of the vector θ , containing the random variables in the KL expansion of Eq. (5), μ l θ , is updated at iteration step l with the information from the vector of measurements d by the following expression: where h l is the output of the model with parameters μ l θ , i.e. h l = G(μ l θ ), and J l is the sensitivity of the measurements with respect to the model parameters. The initial value of the mean is taken as the prior mean μ 0 θ = 0. The covariance matrix of the linearized model outputs C l hh is calculated at each step l with J l and the initial parameter covariance matrix as: where the matrix C 0 θθ is set to the unconditional parameter covariance matrix, C 0 θθ = I. The cross-covariance between model response and parameters is obtained as: The iterative updating process is continued until there is no significant change in the estimate of μ θ . Thereafter, the conditional covariance matrix of the parameters is estimated as follows: Conventionally, the sensitivity J l is obtained by the Adjoint State Method, e.g. Cao et al. (2003). The (i, j) component of matrix J l can be written as We define the adjoint vector λ T i = e T i A −1 . λ i is obtained by solving the so-called adjoint system, which is defined by: The adjoint system needs to be solved once for each measurement, but only one time for all parameters. This is advantageous for cases where the number of parameters is significantly larger than the number of measurements. However, at large numbers of pumping tests and corresponding measurements, a large number of adjoint systems needs be solved and hence, the sensitivity evaluation becomes computationally expensive. After solving the adjoint system, the sensitivity of model response h l i with respect to parameter θ j is obtained as: The partial derivative of the stiffness matrix ∂A(θ ) ∂θ j is evaluated at the element level using the expression for the element stiffness matrix and the KL expansion of Eq. (5); the element matrices are then assembled to obtain the global derivative matrix. The original IEKF algorithm is shown in Algorithm 1. We remark that the final estimate of μ θ obtained by the IEKF method is also a local mode of the posterior density of θ (Jazwinski 2007

IEKF with reduced measurement components
The computational cost of the IEKF method scales linearly with the number of measurements m. Assuming a direct solver for the sparse matrix in the adjoint system, the cost per iteration is O(n 2 h m), where n h is the number of degrees of freedom in the finite element system (Davis 2006). This is because at each iteration, m adjoint systems need to be solved. To enhance its efficiency, we combine IEKF with the measurement reduction approach, introduced in Sect. 3. That is, we decompose the covariance matrix C hh and project the inversion problem on the space defined by the eigenvectors V r corresponding to the m r largest eigenvalues of C hh . As a first approximation, we compute the covariance matrix through its first-order approximation at the initial point C 0 hh , following Eq. (24). The sensitivity J l r of the projected model responses a l r = V T r h l is evaluated with the adjoint method. Let us consider the case, where the data come from a single pumping test. In such case, the modified adjoint system for evaluating J l r reads: with λ T r,i = v T i E T , v i is the eigenvector corresponding to the i-th projected response a r,i and E = [e 1 , . . . , e m ] collects the basis vectors corresponding to each measurement location. The sensitivity is then obtained as: The cost of the Jacobian evaluation is O(n 2 h m r ), which can be significantly lower than for the case where all measurements are used, since usually m r m. Having evaluated the matrix J l r , the covariance of the projected model responses is calculated through: The cross-covariance between projected response and parameters is obtained as follows: We are then able to perform the updating by means of the reduced measurement components d r as follows: The parameter covariance matrix is updated through: To encode the information relevant for updating the current state of the system, we periodically update the projection based on updated estimates of the full covariance matrix. This requires evaluating the Jacobian matrix of the original response, through application of Eq. (27). To avoid a large number of Jacobian evaluations with the original response variables, we only update the projection if the available eigenvector matrix of the reduced measurement components is sufficiently non-orthogonal with respect to the current estimate of the covariance matrix C l aa . We assess this through examining the off-diagonal terms of the following matrix: where the eigenvector matrix V a is obtained by the eigendecomposition of C aa = V a DV T a . If the ratio of the norm of the off-diagonal terms D − diag(D) divided by the norm of the diagonal terms diag(D) surpasses a threshold value η, we estimate the covariance of the original response and update the projection. On this occasion, an updating step is performed with the full measurement set. The algorithm is terminated, when the change between two solutions is smaller than a threshold, and a final updating step is performed with all measurements. The modified IEKF method is shown in Algorithm 2.
We remark that for the case where p pumping tests are performed, a total of p adjoint systems need to be solved per eigenvector for evaluating the Jacobian matrix at each iteration. Therefore, for large p, the computational cost for estimating the Jacobian becomes considerable, even for cases where the number of measurements are reduced significantly. This can be alleviated through decomposing the covariance matrix of the model responses separately for each pumping test and aggregating the resulting eigenvectors in V r .

Sequential Monte Carlo (SMC)
SMC samplers is a category of methods that generate samples from a target distribution through sequentially sampling a set of intermediate distributions (Neal 2001;Chopin 2002;Del Moral et al. 2006). Consider a density sequence { f l (θ), l = 0, . . . , L}, where f 0 is a density that can be easily sampled from, and f L is the target density. Each density f l , l > 0, is known up to a normalizing constant, i.e. f l (θ ) = c −1 l η l (θ ), and supp( f l ) ⊆ supp( f l−1 ). Assume that at step l −1, a set of samples (or particles) {θ j , j = 1, . . . , N } is available that provide a sample approximation of f l−1 . These samples are used to construct a weighted sample approximation of f l through importance sampling: where δ θ j is the Dirac mass at θ j and the weights w j are evaluated as: To avoid degeneracy of the weights as we move towards the target distribution, the weighted samples from f l are transformed to uniformly weighted samples through a resample-move scheme. First, we resample the samples {θ j , j = 1, . . . , N } with probability w j assigned to each sample. We then move the resulting samples through applying a Markov move with a Markov kernel that is reversible with respect to f l . In the context of the Bayesian inversion problem of Eq. (8), the density f 0 (θ) is the prior density f (θ) and the target density f L (θ ) is the posterior density f (θ ) with unnormalized form η L (θ ) = L(θ) f 0 (θ). The choice of the intermediate densities should be such that each density f l−1 provides a good proposal density for sampling from f l . A standard choice is to introduce a set of tempering parameters {β l , l = 0, . . . , L} such that 0 = β 0 < . . . < β L = 1, and define each density f l as Neal (2001): Following Eq. (39), the unnormalized weights W j take the form: To ensure that each pair of densities f l−1 and f l do not vary significantly from one another, each tempering parameter β l can be chosen adaptively such that the effective sample size ESS takes a target value (Jasra et al. 2011). ESS is defined by: That is, at each iteration l, the following optimization problem is solved: where τ > 0 is a target value. We remark that the solution of the optimization problem of Eq. (42) does not require additional model evaluations. Therefore, the contribution to the overall computational time is negligible. For application of SMC in high-dimensional parameter spaces, it is essential to perform the move step with a Markov kernel whose acceptance probability does not degenerate with increase of the parameter dimension. Here, we apply the preconditioned Crank-Nicolson (pCN)-Metropolis-Hastings algorithm (Cotter et al. 2013). The pCN algorithm for sampling from f l generates each proposed state ν from a proposal that is reversible with respect to the prior f 0 . For the present case, where the prior is independent standard normal, N (0, I), the proposal is given by where θ denotes the current state. The candidate is accepted with probability α(θ , ν) given by: The parameter b ∈ [0, 1] controls the variance of the proposal distribution; (1 − b 2 ) −1 is the correlation coefficient between candidate and current states (Papaioannou et al. 2015). We choose the parameter b adaptively to ensure a target acceptance probability close to 0.4. Details on the implementation of the adaptive version of the pCN algorithm can be found in Papaioannou et al. (2015) and Papaioannou et al. (2016). The SMC algorithm for sampling the from the posterior distribution of the inverse problem of Eq. (8) is summarized in Algorithm 3.

SMC with reduced measurement components
In problems where the variance of the error in Eq. (8) is small or where the number of measurements is high, the computational cost of SMC is high because a large number of intermediate sampling steps is required to reach the posterior distribution. Evaluate the weights W l through Eq. (40) and normalize them → w l 9: Resample {θ l−1 , w l } 10: Move the samples by a pCN move ] to get samples θ l 11: while β l < 1 To improve the performance of SMC for such problems, we combine the method with the measurement reduction of Sect. 3. We first estimate the covariance matrix C 0 hh using the samples from the prior distribution at the initial step of the SMC method. The estimate of the covariance matrix is as follows: where h i is the output of the model with parameters θ i , i.e. h i = G(θ i ), and θ i ∼  N (0, I). The matrix C 0 hh is decomposed to obtain the projection V 0 r . SMC is then applied to obtain samples from the posterior associated with the reduced measurement components: The unnormalized weights with the reduced measurement components are evaluated as follows: where 0 r (θ ) is the reduced potential given in Eq. (22) with V r set to V 0 r . Moreover, the acceptance probability of the pCN algorithm is given by: The SMC algorithm with reduced measurement components is presented in Algorithm 4. We remark that in contrast to the IEKF approach with reduced measurement components, we do not update the estimate of the covariance matrix in the modified SMC approach. That is, we only estimate C hh once with samples from the prior distribution of θ . An alternative approach would be to update the estimate of the covariance matrix with samples from the posterior distribution based on the first projection and apply a bridging step (e.g. see Gelman and Meng 1998) to obtain samples from the next posterior distribution based on the new projection. We tested this procedure and observed that the additional bridging step may result in an increase of the total number of sampling steps and, hence, in an increased computational cost of the method, without a significant improvement of the obtained parameter estimate.

Model and parameters
We investigate the performance of the proposed algorithms with a numerical example. A hypothetical aquifer is simulated by a two-dimensional finite element groundwater model. The model specifications are taken from a numerical experiment presented in Huang et al. (2011), which is based on an actual field experiment. The model has a grid with 21 by 21 elements with a cell size of 1 by 1 meter (m).
The prior distribution of the log-hydraulic conductivity log(K) is modeled by a Gaussian random field with mean value of −6.2 and variance of 1.6. The spatial correlation structure is assumed to be exponential, hence the correlation ρ i j between the conductivity in element i and j is given by: where x i j is the Euclidean distance between the midpoints of elements i and j and λ is the correlation length, which is set to 5 m in this example. The conductivities at the elements log(K) are modeled with a truncated KL expansion [cf. Eq. (5)], and we choose the truncation order k to ensure that 99% of the prior parameter variability is represented. This leads to k = 236 instead of n = 441 estimation parameters. 10 monitoring wells are placed in the model, the locations within the model are similar to the real field setup. A total of 7 hydraulic pumping tests is simulated with these wells, and, without consideration of the drawdown at the respective pumping well, a Fig. 1 Model mesh, monitoring and pumping wells total of 70 steady state measurements are generated. The pumping rate for each test is set to 1. An overview of the mesh and the locations of pumping and monitoring wells is shown in Fig. 1. The true parameter values are generated by sampling from the prior distribution of log(K). The observations d are given by the model evaluation of the true parameter superimposed by a realization of a Gaussian measurement noise ∼ N (0, σ I), with noise level of σ = 0.05. The realization of the log(K) field is shown in Fig. 2a.
To assess the quality of the results, we use two metrics. The first is the so-called relative error, related to the estimation parameters, which we define by: where D KL is the diagonal matrix of the eigenvalues in the KL expansion, μ θ is the estimate of the posterior mean of the random variables obtained by the respective method and θ true are the true variables. The second metric we introduce is the relative misfit between model output and measurements: For the SMC, the statistics of the errors are evaluated based on 50 independent simulation runs.

IEKF with all measurements
As a base case, the complete set of 70 measurements from the numerical model are used for inversion. The convergence threshold δ is set to 0.001. The posterior of log(K) is shown in Fig. 2b and the match of the model outputs to the measurements is shown in Fig. 3a. By observation, IEKF is clearly capable of identifying the high and low K regions in-between the monitoring wells. As shown in Table 1, the minimum relative error achieved is 0.8731 and the final relative misfit is 3.6 × 10 −5 . For obtaining this result, one gradient calculation per iteration and the solution of one adjoint system per measurement including the subsequent vector-matrix calculations on the element level were required. This results in 70 measurements × 8 iterations = 560 adjoint system evaluations.

IEKF with reduced measurement components
We now conduct the inversion with Algorithm 2. To obtain the projection, we separate the calculation and eigendecomposition of the observation covariance matrix for each of the 7 pumping tests. This is computationally favorable, since it decouples the calculation of the Jacobian matrix with respect to the eigenmodes and therefore avoids the solution of several adjoint systems per output sensitivity. The estimated mean of log(K), μ θ , is shown in Fig. 2c and the match of the model outputs G(μ θ ) to the measurements d is shown in Fig. 3b. The accuracy of the result decreases slightly compared to the solution in Fig. 2b, but the general features of the low and high K regions match well with the true solution. The input parameters defined for this algorithm and the obtained outputs are given in Table 1 adjoint system evaluations was required, which is about 38 % of the effort of the case presented in Sect. 6.2.1, while losing accuracy with an increase of the relative error by about 1 %. The relative misfit increases from 3.6 × 10 −5 to 8.0 × 10 −4 , which is also apparent when comparing Figs. 3a, b. We remark that the IEKF approach with reduced measurement components has an additional computational cost due to the eigenmode decomposition of up to O(m 3 )  (Parlett 2000) per decomposition step, whereas the Jacobian calculation has complexity of O(n 2 h m r ). That is, the additional cost of the eigenmode decomposition becomes important in cases where m n h , (at least > n 2 h ). This will hardly occur in problems that involve computationally intensive models (i.e. where n h is large). In the problem at hand, it is m = 70 and n h = 484 and, hence, this additional cost is negligible.
We conduct a parameter study to explore the influence of α through varying α between 0.75 and 0.99. The results, as presented in Fig. 4, show that the number of model evaluations increases with increasing α, since the number of final eigenmodes increases from 7 to 32, but, as expected, the model error is reduced from 3% to a minimum of 0.1% compared to the base case with all measurements.

SMC
The performance of the SMC algorithm is tested with the numerical example introduced in Sect. 6.1. Besides the relative error and the relative misfit, the number of sampling levels (tempering steps) is reported as a measure of the computational effort of the method. For all the study cases, except the study of the influence of number

SMC with all measurements
As a base case, the SMC was run for all 70 measurements. Figure 5a shows the average SMC result of all 50 independent simulation runs and Fig. 6a the model output using the average result. The average relative error obtained is 1.194, which is significantly larger than the corresponding IEKF result using all measurements. The relative misfit of the SMC solution is two magnitudes larger than the relative misfit of the IEKF solution. The number of levels averages to 39.6.

SMC with reduced measurement components
Figures 5b and 6b show the solution applying Algorithm 4 with α set to 0.95. The relative error is 1% higher, and the relative misfit is increases from 1.4 × 10 − 3 to 5.8 × 10 − 3. The inversion is about 10 % faster.

Parameter study
We conducted a parameter study to evaluate the influence of the most important parameters in Algorithm 4. It was performed with a Monte Carlo realization with smaller parameter variance (0.1 instead of 1.6) to limit the computation time for the different cases. A reference solution with IEKF was generated and is presented in Fig. 7a. The corresponding SMC solution is presented in Fig. 7b.
In order to focus the investigation on the influence of the number of samples, instead of obtaining the observation covariance through sampling, SMC was run with constant model output eigenmodes obtained from the IEKF solution. Hence, the sampling error in the estimation of the reduced measurement components is eliminated. The model output truncation was chosen such that 99% of the model output variability was preserved (i.e., α = 0.99). As presented in Fig. 8, increasing the number of samples Influence of truncation threshold α results in a better estimate, since the relative error decreases on average. However, one observes that the variability of the error remains constant as the samples per level increase further than 1500. The average relative misfit and its variability decrease until the number of samples per level reach 1500. Using larger sample sets, the estimate and variability of the relative misfit remain constant.
Furthermore, we vary the threshold α which is used to determine the number of eigenmodes in the projection. The relative error, the relative misfit and the number of levels over the different α are presented in Fig. 9. Using more eigenmodes decreases the relative misfit, but increases the variability of the relative error. This is related to the MCMC algorithm employed within each sampling step of the SMC approach. As the number of levels grows, the subsequent MCMC steps result in progressively correlated samples, which increases the variance of the posterior estimates based on the samples in the final level. Therefore, using less eigenmodes improves the solution, since it results in less tempering steps, and hence in a reduced uncertainty of the parameter estimates.

Concluding remarks
This paper presents a method that compresses measurement information from pumping tests and therefore accelerates groundwater parameter inversions. The method decomposes the model output covariance matrix and projects the measurements and model outputs on the eigenvector space corresponding to the largest eigenvalues. Therefore, it incorporates only essential information in the inversion process. The method is combined with the IEKF and SMC algorithms and modified versions of the algorithms are proposed that include reduced measurement information. The IEKF algorithm, which linearizes the relationship between measurement perturbation and parameters, relies on repeated evaluations of the sensitivity (Jacobian) matrix. The computational effort for evaluating the sensitivity matrix increases linearly with increase of the number of measurements. SMC, a sampling-based inversion method with importance sampling and move steps, yields samples from the posterior distribution of the parameters. The width of the likelihood function, describing the measurement information, decreases with increase of the measurement dimension, which leads to an increase of the number of sampling steps. Hence, a measurement reduction has potential for improvement of both inversion methods. As demonstrated by a numerical example, the inversion is generally faster than using all measurements with negligible loss of solution accuracy, provided that the parameters of the proposed algorithms are carefully chosen.
The benefit of the proposed approach increases with increase of both the size of the computational model and the number of measurements. For the IEKF method, this is because the computational cost per adjoint system evaluation is O(n 2 h m), where n h is the number of degrees of freedom of the computational model and m is the number of measurements. We note that the computational cost per eigenmode decomposition is up to O(m 3 ), so the computational benefit of the measurement reduction is more pronounced in cases where m n h . For large-scale inversion problems (large n h ) we expect that the additional computational cost of the eigendecomposition will be small compared to the computational gain. Another crucial aspect for an improved inversion process is the structure that governs the eigendecomposition. The standard approach is to collect the full model response dataset, form a covariance matrix, and then perform the decomposition. This implies that each eigenvector contains information on all pumping tests together. When the IEKF method is applied with relatively large numbers of pumping tests and small numbers of data per pumping test, the eigendecomposition should be conducted separately for each pumping test, as this decreases the cost required for evaluating the sensitivity matrix of the projected responses.
Another aspect we address in the paper is the frequency of the covariance calculations for obtaining an optimal measurement reduction. Since the prior distribution of the model parameters can differ considerably from their posterior distribution, the measurement reduction based on the prior distribution might not reflect the true behavior of the measurement correlation. Therefore, in the modified IEKF algorithm, we propose to update the response covariance matrix frequently throughout the inversion process, whenever the projection becomes too inaccurate (see Algorithm 2). Regarding the modified SMC algorithm, we observed that updating the covariance matrix leads to additional sampling steps and, hence, increased computational costs compared to the standard SMC algorithm. Hence, we recommend to form the observation covariance matrix with samples from the prior conductivity distribution and employ this first estimate in the updating process (see Algorithm 4).
Comparing the computational cost of the two inversion methods in the worked numerical example, we observe that the IEKF, although being a linearization method, performed significantly better than the SMC approach. In particular, the high variance of the estimated K field did not seem to impose a large problem for the IEKF approach, which obtained an accurate parameter estimate in a small number of iterations. In contrast, the SMC method required more sampling iterations and resulted in parameter estimates with significant associated uncertainty. The modified SMC method could be improved by changing the move step through the implementation of a more efficient MCMC algorithm, e.g. which includes gradient information (Roberts et al. 1996;Cui et al. 2016;Rudolf and Sprungk 2018). This could potentially decrease the required samples per level to obtain an accurate approximation of the posterior. Another sampling-based inversion method worth investigating in this context is the ensemble Kalman filter. This method has recently found popularity in Bayesian inverse problems and its performance can potentially profit from a reduced measurement set.