Bayesian Poroelastic Aquifer Characterization from InSAR Surface Deformation Data Part II: Quantifying the Uncertainty
Amal Alghamdi, Marc Hesse, Jingyi Chen, Umberto Villa, Omar Ghattas
mmanuscript submitted to
Water Resources Research
Bayesian Poroelastic Aquifer Characterization fromInSAR Surface Deformation DataPart II: Quantifying the Uncertainty
Amal Alghamdi , Marc A. Hesse , , Jingyi Chen , , Umberto Villa , OmarGhattas , , University of Texas at Austin, Oden Institute for Computational Engineering and Sciences, Austin, TX,United States University of Texas at Austin, Geological Sciences, Austin, TX, United States University of Texas at Austin, Aerospace Engineering & Engineering Mechanics, Austin, TX, UnitedStates Washington University in St. Louis, Electrical and Systems Engineering, St. Louis, MO, United States University of Texas at Austin, Mechanical Engineering, Austin, TX, United States
Key Points: • Using InSAR data reduces the uncertainty in selected quantities of interestcompared to using prior knowledge only • The preconditioned Crank–Nicolson (pCN) MCMC method is extended to ex-ploit posterior curvature and allow better chain mixing • We demonstrate the intrinsic low dimensionality of the poroelastic inverse prob-lem that is critical for the success of the MCMC method
Corresponding author: Amal Alghamdi, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] F e b anuscript submitted to Water Resources Research
Abstract
Uncertainty quantification of groundwater (GW) aquifer parameters is critical for ef-ficient management and sustainable extraction of GW resources. These uncertaintiesare introduced by the data, model, and prior information on the parameters. Here wedevelop a Bayesian inversion framework that uses Interferometric Synthetic ApertureRadar (InSAR) surface deformation data to infer the laterally heterogeneous perme-ability of a transient linear poroelastic model of a confined GW aquifer. The Bayesiansolution of this inverse problem takes the form of a posterior probability density ofthe permeability. Exploring this posterior using classical Markov chain Monte Carlo(MCMC) methods is computationally prohibitive due to the large dimension of thediscretized permeability field and the expense of solving the poroelastic forward prob-lem. However, in many partial differential equation (PDE)-based Bayesian inversionproblems, the data are only informative in a few directions in parameter space. Forthe poroelasticity problem, we prove this property theoretically for a one-dimensionalproblem and demonstrate it numerically for a three-dimensional aquifer model. Wedesign a generalized preconditioned Crank–Nicolson (gpCN) MCMC method that ex-ploits this intrinsic low dimensionality by using a low-rank based Laplace approxima-tion of the posterior as a proposal, which we build scalably. The feasibility of ourapproach is demonstrated through a real GW aquifer test in Nevada. The inherentlytwo dimensional nature of InSAR surface deformation data informs a sufficient numberof modes of the permeability field to allow detection of major structures within theaquifer, significantly reducing the uncertainty in the pressure and the displacementquantities of interest.
Efficient groundwater (GW) resources management is critical for mitigating theglobal-scale problem of GW depletion (Wada et al., 2010, 2012, 2014; Famiglietti,2014). This is the result of extracting GW at rates that exceed the natural recharge,which then leads to significant drops in GW pressure and (possibly irreversible) com-paction (Wada et al., 2010; Konikow & Kendy, 2005; Xue et al., 2005; Holzer &Galloway, 2005). Making informed management and control decisions requires accu-rate GW aquifer models, typically in the form of a time-dependent partial differentialequation (PDE). These models are used to predict the aquifer response (e.g. GW pres-sure drop and aquifer compaction) to GW extraction activities. The major source ofuncertainty in these models is the uncertainty in the aquifer properties (Eaton, 2006;Bohling & Butler, 2010; Oliver & Chen, 2011). The goal of our work is to devise aframework to systematically quantify the uncertainty in GW aquifer model properties.Here we only consider uncertainty in one parameter field, the aquifer permeability,and one data source, surface deformation data, but the framework we present can begeneralized to multiple parameter fields and data sources. To estimate the perme-ability, we take advantage, in particular, of the evergrowing InSAR data that providefrequent, large-scale, and up to sub-millimeter accurate surface deformation measure-ments globally (Ferretti et al., 2007; Tom´as et al., 2014).Bayes’ theorem—which has been adopted in GW applications since as early as1986 (Carrera & Neuman, 1986a, 1986b, 1986c; McLaughlin & Townley, 1996; Linde etal., 2017)—provides a natural means of incorporating observational data, e.g., pressureor surface deformation data at the aquifer site, to quantify the uncertainty in theaquifer model parameters (such as physical properties, source terms, and/or boundaryor initial conditions). The Bayesian framework updates our assumed “prior” knowledgeabout the aquifer parameters through what is known as the “likelihood” distribution,which determines how likely it is for the observational data to result from the modeledsystem with a particular parameter realization. The updated probability distributionis called the “posterior” distribution and is regarded as the solution of the inverse –2–anuscript submitted to
Water Resources Research problem. Bayesian inversion, therefore, provides a characterization of the uncertaintyin the parameters, as opposed to finding just a point estimate as is with deterministicinversion.The Markov chain Monte Carlo (MCMC) method is often the method of choicefor characterizing the posterior distribution for PDE-based Bayesian inverse problems.MCMC generates samples of the posterior probability density function (PDF)—eachof which requires solution of at least the forward PDEs—from which sample statisticscan be computed. Conventional MCMC methods view the parameter-to-observablemap as a black box and thus are not capable of exploiting the structure of this map toaccelerate convergence of MCMC chains. As such, most work on Bayesian inversiongoverned by PDEs has been restricted to simple PDEs or low parameter dimensions(or both).Over the past decade, several powerful MCMC methods have emerged that ex-ploit the geometry and smoothness—and in some cases the underlying low-dimensionality—of the parameter-to-observable map to enable efficient MCMC sampling from high-dimensional PDE-based posteriors (Cui et al., 2016; Flath et al., 2011; Martin et al.,2012; Bui-Thanh et al., 2013; Petra et al., 2014; Beskos et al., 2017). These methodshave been applied to several large-scale Bayesian inverse problems governed by com-plex forward problems, including seismic wave propagation (Martin et al., 2012) andice sheet flow (Isaac et al., 2015). What these methods all have in common is the useof the Hessian of the negative log posterior to exploit the geometry of the posterior.Moreover, the low rank structure of the Hessian can be used to identify the (often)low dimensional manifold (of dimension r much smaller than the nominal dimensionof the discretized parameters) on which the parameters are informed by the data. Tothat end, randomized methods (e.g., randomized singular value decomposition and ran-domized generalized eigensolvers) are devised to identify this low dimensional manifoldat a cost of O ( r ) PDE solves only. Characterizing high-dimensional posterior distri-butions governed by three-dimensional transient poroelastic GW aquifer models hasbeen explored by Hesse and Stadler (2014) via a Gaussian approximation of the pos-terior distribution for model problems. Full MCMC sampling of posteriors governedby high-dimensional GW models of such complexity has yet to be explored.In Alghamdi et al. (2020), Part I of this two-part series of articles, we built aBayesian framework to characterize the lateral heterogeneity in a three-dimensionaltime-dependent poroelastic aquifer model and applied it to a test case in Nevadausing InSAR surface deformation data. There, we addressed the choice of the priordistribution and the data noise model. We also took a first step in exploring theposterior by inverting for the maximum a posteriori (MAP) point, which is the mostlikely permeability realization with respect to the posterior distribution. We usedindependent GPS measurements (Burbey et al., 2006) for validation. A brief summaryof the pertinent results from Part I is provided in section 3.We summarize the contribution of Part II as follows: (1) We demonstrate the de-caying eigenvalues of the data misfit part of the Hessian analytically, through Fourieranalysis, for a simplified 1D poroelasticity problem. We additionally provide numericalevidence of the data misfit Hessian compactness for the three-dimensional fully-coupledquasi-static linear poroelastic aquifer model considered here. This property determinesthe applicability of the aforementioned class of scalable MCMC methods. (2) We pro-pose a generalization of the preconditioned Crank–Nicolson (pCN) MCMC method(Cotter et al., 2012), referred to as gpCN thereafter, that identifies and exploits thislow dimensional manifold using randomized eigensolvers to construct a low-rank basedLaplace approximation of the posterior. Besides the initial cost of creating the low-rank based Laplace approximation, which is negligible in comparison to the cost ofgenerating MCMC chains, gpCN has no additional cost per MCMC sample in com-parison to pCN. Furthermore, we compare convergence properties of pCN and gpCN. –3–anuscript submitted to Water Resources Research (3) We apply gpCN to characterize the high-dimensional subsurface permeability forthe Nevada test case aquifer, and quantify the associated uncertainty (both in theparameters and the selected state-variable-based quantities of interest (QoIs)). Wediscuss how information from the data reduces the prior uncertainty in the parametersand the state variables in terms of the location in the physical domain and the dom-inant parameters modes inferred from the data. The three integrated contributionsdemonstrate a successful design and application of a scalable MCMC method to areal-world GW Bayesian inverse problem governed by a coupled forward problem withhigh dimensional parameters and using InSAR data.The reminder of the paper is structured as follows. In section 2, we review thequasi-static linear poroelasticity model and the Bayesian framework that we estab-lished in Part I. We discuss low-rank based Laplace approximation of the posteriorand present gpCN, our generalization of pCN MCMC method. We describe the testcase to which we apply our framework in section 3. We demonstrate gpCN perfor-mance and present results on characterizing the uncertainty in the permeability andselected state-variable-based QoI for the Nevada test case in section 4. We summarizeour conclusions in section 5.
In Part I of this two-part series of articles, we formulate a Bayesian inversionframework governed by quasi-static linear poroelasticity. The solution of the Bayesianinverse problem takes the form of a posterior distribution that characterizes the hetero-geneity of groundwater aquifer permeability. Here, we review the forward model, thequasi-static linear poroelasticity model, and the formulation of the inverse problem,the Bayesian framework, in sections 2.1 and 2.2, respectively. We also demonstrate thecompactness of the data misfit Hessian for a simplified one dimensional poroelasticityproblem. In sections 2.3 and 2.4, we build a low-rank based Laplace approximationand present the gpCN MCMC method that exploits this intrinsic low dimensionalityin the poroelasticity inverse problem.
In the Bayesian inversion framework, we model the aquifer as a poroelasticmedium governed by the linear poroelasticity theory developed by Biot (1941). Biottheory describes the coupling between the fluid pressure field in a saturated porousmedium and the accompanying elastic deformation of the solid skeleton. The classicalformulation of Biot system in a space-time domain Ω × (0 , T ] can be written as( S (cid:15) p + α ∇ · u ) t − ∇ · (cid:18) κµ ∇ p (cid:19) = f p −∇ · ( σσσ ( u ) − αp I ) = f u , (1)where ( · ) t denotes time derivative, p = p ( x , t ) is the deviation from hydrostatic pres-sure, u = u ( x , t ) is the displacement of the solid skeleton, f p = f p ( x , t ) is a fluidsource per unit volume, f u = f u ( x , t ) is the body force per unit volume, S (cid:15) is the spe-cific storage, κ ( x ) is the medium permeability field, µ is the dynamic viscosity of thepore fluid, and α is the Biot–Willis coupling parameter. The tensor σσσ ( u ) is the stresstensor for isotropic linear elasticity, and it is parameterized by the drained shear mod-ulus G and the Poisson’s ratio ν . We write the unknown medium permeability κ ( x ) interms of the log-permeability field, which we denote m ( x ), so that κ ( x ) = e m ( x ) . Thisparameterization ensures that the inferred permeability field is positive when solvingthe inverse problem for m ( x ). The system (1) together with well-defined boundaryand initial conditions, given in Alghamdi et al. (2020), form the forward problem. –4–anuscript submitted to Water Resources Research
We use a three-field formulation of the poroelasticity system (1) in which Darcyflux is introduced as a state variable q ( x , t ) = − κµ ∇ p , adding a third equation tothe system (1). Following Ferronato et al. (2010), we employ a mixed finite elementmethod (MFEM) to discretize the three-field formulation in space and use backwardEuler method for discretization in time. In this MFEM, the pressure p is approx-imated by piecewise constant functions, the displacement u is approximated in thefirst-order Lagrange polynomial space, and the flux q is approximated in the lowest-order Raviart–Thomas space. This MFEM discretization has the advantage of con-serving mass discretely and reducing the non-physical pressure oscillations that canform in some other discretizations of the poroelasticity system (Phillips & Wheeler,2009; Ferronato et al., 2010; Haga et al., 2012). To simplify notation, we define thediscretized “global” space-time system, which combines solving for all time steps si-multaneously, as ¯S ( m ) ¯X = ¯F , (2)where the block lower triangular matrix ¯S combines the discretized differential op-erators of the three-field formulation of the poroelasticity system (1). ¯S depends onthe discretized log permeability m ∈ R n n , which we approximate in the first orderLagrange polynomial space, where n n is the number of degrees of freedom (DOFs) ofthe discretized parameter field. The vector ¯F combines the source terms and the nat-ural boundary conditions at all time steps, and the vector ¯X consists of the pressure,displacement, and fluid flux DOFs at all time steps as well. The discretization detailscan be found in Part I (Alghamdi et al., 2020). Our goal is to characterize the aquifer heterogeneous permeability κ = e m in theporoelasticity model (1) from InSAR Line-of-Sight (LOS) surface deformation data d obs ∈ R n obs , where n obs is the number of observational data (total number of pixelsin deformation maps in this case). The LOS surface deformation u LOS is given by u LOS = α u + α u + α u , (3)where u , u and u are the east, north and vertical surface displacements, respec-tively, and α , α and α are the eastward, northward and vertical components of theLOS unit look vector α LOS = [ α , α , α ]. In this study, we only use one LOS netdeformation map available at t = t InSAR . In general, multiple deformation maps ordeformation time series can be available.Assuming modeling errors are negligible compared to noise in the data, for agiven log permeability field m ( x ), the relation between the observational data d obs and the model output is given by d obs = F ( m ) + ηηη. (4)The map F ( · ) : R n n → R n obs is the parameter-to-observable map which, for a given m , gives the observables, i.e the modeled LOS surface deformation values at the ob-servational data locations (pixels centers) and times. Evaluating F ( m ) requires solv-ing model (1) using the given log permeability m . In particular, F ( m ) := ¯B ¯X = ¯B¯S ( m ) − ¯F , where we define the pointwise linear observation operator ¯B which ex-tracts the pointwise observables from the modeled displacement DOFs included in theglobal vector ¯X . The solution ¯X is obtained by solving the “global” system (2). Therandom vector ηηη ∼ N ( , Γ noise ) is measurement noise which we assume is additive andGaussian with mean and covariance matrix Γ noise ∈ R n obs × n obs .We use relation (4) to define the likelihood probability distribution that deter-mines how likely it is for the observational data d obs to be measured from an aquifer –5–anuscript submitted to Water Resources Research with permeability m —in the absence of model error. In other words, the likelihooddistribution is the probability distribution of the discrepancy between the observeddata and the model outcome which, based on the relation (4), is given by π likelihood ( d obs | m ) ∝ e − (cid:96) ( d obs ; m ) , (5)where the symbol ∝ denotes equality up to a multiplicative constant and (cid:96) is thenegative log of the likelihood, or the data misfit, which is given by (cid:96) ( d obs ; m ) = 12 ( F ( m ) − d obs ) T Γ − ( F ( m ) − d obs ) . (6)Bayes theory gives the posterior probability distribution π post ( m | d obs ), whichquantifies the uncertainty in the parameter m given the observational data d obs . Theposterior distribution incorporates information about the model and the data, via thelikelihood distribution. Statistical prior assumptions about the parameters are encodedinto the prior distribution π prior ( m ). Bayesian formulation of inverse problems wasintroduced in (Tarantola et al., 1982) and is given by π post ( m | d obs ) ∝ π likelihood ( d obs | m ) π prior ( m ) ∝ exp (cid:18) −
12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs ) −
12 ( m − ¯ m ) T Γ − ( m − ¯ m ) (cid:19) . (7)We assume the prior distribution is a discretized Mat´ern Gaussian field with mean¯ m ∈ R n n and covariance matrix Γ prior ∈ R n n × n n . The latter dictates our assumptionsabout the pointwise variance and the spatial correlation of the log permeability fieldfeatures. For this class of priors, the precision operator Γ − admits a representation Γ − = AM − m A in terms of an elliptic differential operator A and a mass matrix M m in the discretized parameter space (Lindgren et al., 2011; Alghamdi et al., 2020).The posterior distribution is a distribution in a very large dimensional space R n n ,and evaluating the posterior at a single realization m requires solving the model (1)forward in time. Computing moments and other expected values of interest from theposterior distribution is computationally prohibitive using traditional MCMC samplingmethods. Fortunately, in many PDE-based Bayesian inversion problems, the dataare informative in only relatively few directions in parameter space r (cid:28) n n . Thisproperty can be exploited to obtain low-rank based approximations of the posterior,which additionally can be used to build effective MCMC proposals, as we present inthe remaining of this section. In Part I we focused on solution of the optimization problem of maximizingthe posterior distribution to find the MAP point m MAP . Besides being the mostlikely realization of the permeability field, the MAP point is also needed to form theLaplace approximation of the posterior, which is a Gaussian approximation obtainedby making a quadratic approximation of the negative log posterior at the MAP point(MacKay, 2003). This approximation is valid if the parameter-to-observable map isapproximately linear over the support of the prior distribution. Additionally, it canbe used in building effective proposals for dimension invariant MCMC algorithms aswe discuss in section 2.4.The negative log of the posterior (7) can be written as − log ( π post ( m ))= 12 (cid:0) F ( m ) − d obs (cid:1) T Γ − (cid:0) F ( m ) − d obs (cid:1) + 12 ( m − ¯ m ) T Γ − ( m − ¯ m ) − c , (8) –6–anuscript submitted to Water Resources Research where c is the log of the posterior normalization constant. Using a Taylor expansionat the point m MAP in parameter space, a quadratic approximation of (8) is given by − log ( π post ( m MAP + ˆ m )) ≈ c + 12 ˆ m T H ( m MAP ) ˆ m , (9)where ˆ m = m − m MAP . To obtain (9), we use the fact that the gradient of the negativelog posterior with respect to the parameter m is zero at the MAP point. The constant c is defined as c = 12 (cid:0) F ( m MAP ) − d obs (cid:1) T Γ − (cid:0) F ( m MAP ) − d obs (cid:1) + 12 ( m MAP − ¯ m ) T Γ − ( m MAP − ¯ m ) − c . (10)The matrix H ( m MAP ) ∈ R n n × n n is the Hessian operator of the negative log of theposterior (7) evaluated at the MAP point and it is given by H ( m ) = H misfit ( m ) + Γ − , (11)where the data misfit Hessian H misfit ( m ) is given by H misfit ( m ) = ∇ F ( m ) T Γ − ∇ F ( m ) + n obs (cid:88) i (cid:0) Γ − ( F ( m ) − d obs ) (cid:1) i ∇ F i ( m ) . (12)The resulting Laplace approximation is the Gaussian distribution: π Laplace ( m ) ∼N (cid:0) m MAP , H ( m MAP ) − (cid:1) . (13) In many PDE-based inverse problems, it is often the case that the data misfitHessian operator is compact and only a few parameter dimensions, O ( r ), are informedby the data (Flath et al., 2011; Bui-Thanh et al., 2012; Isaac et al., 2015; Villa et al.,2020). In Alghamdi (2020), we provide an analytical formula for the eigenvalues andthe eigenfunctions of the Hessian of an idealized one-dimensional steady-state poroe-lasticity problem. We assume this problem is defined on an interval [0 , L ], and that wehave continuous deformation observations everywhere in this interval. Specifically, theeigenvalues λ i and eigenfunctions φ i of the data misfit Hessian operator in this caseare given by λ i = f L µ e κ (cid:18) Liπ (cid:19) , φ i ( x ) = cos (cid:18) iπx L (cid:19) , i = 1 , , , ..., (14)where f L is the Darcy flux at the boundary x = L , e is Young’s modulus, and κ isthe constant permeability value at which the data misfit Hessian operator is evalu-ated. We note from the expressions in (14) that the eigenvalues have a rapid decayof O ( i ), and the eigenfunctions corresponding to dominant eigenvalues are smoother.This means that the more oscillatory a permeability mode is, the more difficult it isto infer from the deformation data. Note that the expected information gain (fromprior to posterior) from the data in each mode is given by log(1 + λ i ) under a Gaussianassumption (Alexanderian et al., 2016). The rapid decay of eigenvalues in (14) impliesthat the data inform the permeability field in a low dimensional manifold. Althoughthis analysis is carried out for an idealized one-dimensional case, numerical evidencesuggests rapid eigenvalue decay—independent of the discretization dimension—of thedata misfit Hessian of the 3D poroelasticity problem we study; see section 4. More-over, analogous to the 1D analysis result, larger eigenvalues correspond to smoothereigenvectors. –7–anuscript submitted to Water Resources Research
The gpCN method described in the next section requires the ability to gener-ate samples ηηη
Laplace from the distribution π Laplace ( m ). To sample from this Gaussiandistribution, H − is applied to a white noise vector. However, constructing the Hes-sian matrix explicitly and computing its Cholesky factorization is computationallyprohibitive in high parameter dimensions. The Hessian construction costs n n pairsof incremental forward (A8) and incremental adjoint (A9) poroelasticity solves; seeAppendix A for details.Instead of building the Hessian explicitly, we take advantage of the fact thatapplying the Hessian to a vector can be carried out in a matrix-free manner and costsonly a pair of incremental forward and incremental adjoint solves, and the data misfitHessian is inherently low-rank. This allows us to build a low-rank based representationof the Hessian at m MAP by solving a generalized eigenvalue problem using a random-ized generalized eigensolver—at a cost of just O ( r ) incremental forward and adjointporoelasticity solves. We use this approximation to create a low-rank based Laplaceapproximation, π Laplace,r ( m ), of the posterior.The following approach of approximating the Hessian at the MAP point hasbeen applied to the inverse problem for the basal friction coefficient of the Antarcticice sheet (Isaac et al., 2015), and is integrated in the hIPPYlib library (Villa et al.,2018). This approach is based on solving the generalized eigenvalue problem: H misfit ( m MAP ) w i = λ i Γ − w i for i = 1 , .., r. (15)The eigenvalues λ i of this generalized eigenvalue problem indicate whether the cor-responding eigendirection w i in the high dimensional parameter space R n n is mostlyinformed by the data (when λ i >
1) or by the prior (when λ i < λ i reflects the level of certainty of inferring the corresponding eigendirectionfrom the data. The diagonal matrix of the first r dominant eigenvalues, λ i , and thematrix of the corresponding eigenvectors, w i , of this generalized eigenvalue problemare denoted as Λ r = diag([ λ , λ , ..., λ r ]) and W r = [ w , .., w r ], respectively. The ma-trix W r has Γ − -orthonormal columns, that is W Tr Γ − W r = I r , where I r denotesthe identity matrix in R r .The eigenvalue problem (15) along with the Sherman–Woodbury formula are thenused to obtain the following low-rank approximation of the Laplace approximationcovariance matrix (the inverse of the Hessian at the MAP point); see (Villa et al.,2020): H ( m MAP ) − = ( H misfit ( m MAP ) + Γ − ) − ≈ Γ prior − W r D r W Tr . (16)The matrix D r is given by D r = diag([ λ λ , λ λ , ..., λ r λ r ]). The value λ i λ i is neg-ligible for λ i (cid:28)
1, which provides a basis for truncation. We denote the low-rankbased covariance matrix given by the inverse of the Hessian in (16) by Γ Laplace,r := Γ prior − W r D r W Tr . A sample ηηη Laplace,r from this low-rank based Laplace approxima-tion π Laplace,r ( m ) = N ( m MAP , Γ Laplace,r ) can be generated scalably using (Villa et al.,2020): ηηη
Laplace,r = ( I − W r S r W Tr Γ − ) ηηη + m MAP (17)where S r = diag([1 − √ λ , − √ λ , ..., − √ λ r ]) and ηηη ∼ N (0 , Γ prior ). The costof sampling from this approximation is dominated by the cost of generating the priorsample ηηη (mainly applying A − via a scalable multigrid solve as discussed in Part I),and the matrix–vector multiplications by sparse and low-rank matrices. –8–anuscript submitted to Water Resources Research
To characterize the posterior distribution, MCMC methods generate samples ofthe posterior PDF (each of which requires at least a forward PDE solve) from whichsample statistics can be computed. As stated earlier, for PDE-governed posteriors in ahigh dimensional parameter space, it is often the case that the low rank structure of theHessian can be used to identify the low dimensional manifold on which the parametersare informed by the data, which limits the expensive PDE solves to exploration of justa low dimensional subspace of the parameter space.In this work, we design a generalization of the preconditioned Crank–Nicolson(gpCN) MCMC method. Preconditioned Crank–Nicolson (pCN) is a dimension-invariantMCMC method that uses prior information to build the proposal distribution (Cotteret al., 2012; Pinski et al., 2015). The gpCN algorithm is also dimension-invariant,yet has the advantage of capturing posterior curvature information in a relatively in-expensive manner. It uses an approximation of the gradient and the low-rank basedapproximation of the Hessian at the MAP point discussed in section 2.3, rather thancomputing exact gradient directions and Hessian matrices at numerous samples in pa-rameter space. The gpCN proposed sample ηηη ( k )prop in the k th Metropolis Hastings (MH)iteration (described in Algorithm 1) is given by ηηη ( k )prop = m MAP + (cid:112) − β ( ηηη ( k − − m MAP ) + βηηη ( k )0 (18)where ηηη ( k )0 ∼ N (0 , Γ Laplace,r ) and ηηη ( k − is the ( k − th MH iteration posterior sample.The proposed sample ηηη ( k )prop is then subject to acceptance or rejection based on the MHalgorithm. The parameter 0 < β < β , the higher the acceptance rate but thelonger the chain takes to explore the posterior distribution and generate independentsamples. Larger β values, however, allow the chain to visit farther points in parameterspace but lead to lower acceptance rates. The cost of each MH iteration is dominatedby the cost of evaluating π post ( ηηη ( k )prop ), which amounts to a single solution of the forwardporoelasticity problem (1) for the log permeability realization ηηη ( k )prop . The gpCN methodis integrated in hIPPYlib (Villa et al., 2018). Algorithm 1
MH algorithm for generalized preconditioned Crank–Nicolson (gpCN)method
Input β , Γ Laplace,r , m MAP , number-of-samples for k = 1 , , , ..., number-of-samples do draw ηηη ( k )0 ∼ N (0 , Γ Laplace,r )calculate ηηη ( k )prop using formula (18) α gpCN ← π post ( ηηη ( k )prop ) π Laplace,r ( ηηη ( k − ) π post ( ηηη ( k − ) π Laplace,r ( ηηη ( k )prop ) draw scalar s ∼ U ([0 , if s < min(1 . , α gpCN ) then ηηη ( k )post ← ηηη ( k )prop else ηηη ( k )post ← ηηη ( k − end ifend for –9–anuscript submitted to Water Resources Research
Figure 1: Overview of the Nevada test case. (a) A conceptual illustration of the aquifermodel. The well is screened in the lowest half of the aquifer (blue segment). (b) InSAR-observed LOS deformation (east, north, and vertical displacements projected onto theradar line-of-sight) between May 4, 2003 and October 26, 2003 over the study area. (c)The computational mesh we use to model the aquifer. The mesh has 4,081 nodes. (d)Two-dimensional horizontal slice of the inferred permeability, e m MAP , at mid-aquifer depth(in decimal logarithm) using InSAR data in Figure 1 b and the mesh in Figure 1 c . In pan-els (b)–(d), x is the east direction and y is the north direction relative to the well location( x, y ) = (0 ,
0) (marked by the red star).
Our study site is located southeast of Mesquite, Nevada. The aquifer is ∼ ∼
85 m thick brittle layer at the top and a claylayer that is at least 400 m thick at the bottom (Figure 1 a ). A controlled aquifertest was performed at a newly installed municipal well WX31 between May 7, 2003and July 9, 2003 (Burbey, 2006). This well pumped 9,028 m of GW per day duringthe study period, and it was relatively isolated from other wells in the area. Theregional GW recharge is due mainly to winter precipitation in the surrounding moun-tains, which was negligible during the summer pumping test. Burbey (2006) analyzedthe aquifer geomechanical response to the pumping test using InSAR and GPS data.They found significant deviations between the observed deformation pattern and theexpected deformation pattern from a homogeneous aquifer. This finding motivates usto investigate the lateral permeability heterogeneity of the aquifer using the Bayesianframework.We processed an Envisat interferogram (05/04/2003-10/26/2003) that capturesthe surface deformation associated with the pumping test along the radar LOS direc-tion (Figure 1 b ). Here the LOS unit look vector is α LOS = [0 . , − . , . d obs is mostly sensitive to vertical motion ( α ≈ b ) shows a subsidence bowl that is shifted toward the southeast of the welland an uplift signal in the northwest of the domain. The latter is likely due to a faultexistence northwest the well (Burbey, 2008; Alghamdi et al., 2020). We found thatthe noise in the Envisat interferogram is mainly due to phase decorrelation that is un-correlated in space ( σ InSAR = 3 . Γ noise to be a diagonal matrix with diagonal values σ .Our aquifer model is described in detail in Alghamdi et al. (2020), we assumea 3D poroelastic medium governed by the Biot model (1) and represent the pumpingwell as a volumetric sink term f p . The parameter values of the forward model (1) areprovided in Table B1 and based on previously estimated values by Burbey (2006). We –10–anuscript submitted to Water Resources Research discretized the model (Figure 1 a ) using a 4,081-node unstructured tetrahedral layeredmesh (Figure 1 c ). The discretization on this mesh results in a total of n state = 72 , n state = 320 ,
824 and n state = 1 , , Γ prior H misfit .We imposed zero normal displacement at the bottom boundary, zero displace-ment at the lateral boundaries, no traction at the top boundary, and no GW flux atall boundaries. Under the assumption of zero initial deviation from the hydrostaticpressure everywhere in the domain at t = 0, we ran the model for T = 175 days using154 variable-length time steps. Initial time steps are 1.2 hours to capture rapid pres-sure decay that occurs after the onset of pumping and increased gradually to reachfive days toward the end of the simulation. Each time step involves solving a linearsystem of size n state . Based on a heuristic approach (Alghamdi et al., 2020), we chose the weakestbilaplacian-like Mat´ern class prior that yet sufficiently regularizes the inverse problemto allow the (Gauss–)Newton method to converge (defined as reducing the gradientnorm by four orders of magnitude). The average-over-domain pointwise standarddeviation (SD) of this prior is ∼ . ρ = 2 km in thelateral direction and an order of magnitude larger in the vertical direction. The latterconstraint inhibits permeability variations in the vertical direction. We impose thisconstraint because (1) vertical variations are likely not informed by deformation dataobtained on the surface; (2) the large aspect ratio of the aquifer horizontal-to-verticaldimensions leads to predominantly horizontal flow.In the Part I paper (Alghamdi et al., 2020), we used the InSAR surface deforma-tion data (Figure 1 b ) to solve for the permeability realization that is most likely to bethe true permeability field. This is known as the maximum a posteriori (MAP) point ofthe posterior distribution (Equation (7)). We found that the MAP point (Figure 1 d )reveals distinct features in the lateral permeability field: high permeability channelextending from south the well to the southwest and low permeability barrier to thenorthwest of the well. We validated the InSAR-based permeability MAP point solutionusing an independent GPS data set as well as a synthesized data set. We computedthe MAP point using (Gauss–)Newton method and observed a mesh-independent con-vergence. We also demonstrated the consistency of the solution with respect to a widerange of prior assumptions, and InSAR data multi-look choices, given a consistentdata noise treatment. The Bayesian inversion framework allows us to quantify the uncertainty in thesolution by sampling the posterior distribution. To estimate the uncertainty in theaquifer permeability, we first build the low-rank based Laplace approximation π Laplace,r of the posterior, a process that entails solving the generalized eigenvalue problem (15)using a randomized generalized eigensolver (Villa et al., 2020), to obtain the low-rankbased approximation of posterior covariance matrix (16). We study two approxima-tions in which we use truncation values r = 80 and r = 150, respectively. This low-rankbased Laplace approximation is feasible to compute—it mainly costs O ( r ) incremental –11–anuscript submitted to Water Resources Research forward (A8) and incremental adjoint (A9) solves—and can serve (depending on howclose the parameter to observable map is to a linear one) as a good approximation ofthe posterior distribution π post . Additionally, we use this approximation in the gpCNproposal (see Algorithm 1) to sample the true posterior.To that end, we generate five MCMC chains, using the gpCN MCMC method(Algorithm 1), each of length 50,000, with additional 1,000 burn-in samples. Eachchain start point ηηη (0)post is an independent sample from the low-rank based Laplaceapproximation π Laplace,r , with r = 80. We similarly generate two pCN chains (Cotteret al., 2012). Generating a single gpCN chain (and similarly a pCN chain) costs50,000 forward solves of the system (1), discretized on the 4,081-node mesh. Thesesolves are required to evaluate the ratio α gpCN in Algorithm 1. We store pressure p and displacement u solutions evaluated at time t InSAR ≈ . b ) is observed relative to the start of thesimulation, for different realizations of m . We set the step size β = 0 .
25 in gpCNproposal formula (18) and set β = 0 .
01 in pCN proposal. These values were foundwith a trial-and-error approach to balance between high acceptance rate and large stepsize. To compare several posterior-based QoI vs prior-based QoI, we generate 1,000samples from the prior distribution and the low-rank based Laplace approximation(with r = 80) directly and solve the system (1) at each sample. We carry out theimplementation using the FEniCS library for finite element discretization in space (Logget al., 2012) and the hIPPYlib library for state-of-the-art Bayesian and deterministicPDE-constrained inversion algorithms (Villa et al., 2018).
Here we first establish the low rank property of the data misfit Hessian for theNevada test-case poroelasticity problem introduced in section 4.1. This property iscritical for the applicability of the gpCN algorithm. Then we study the convergence andmixing of the gpCN algorithm and compare it to the performance of the pCN algorithmin section 4.2. Finally, we use the gpCN algorithm to quantify the uncertainty in theaquifer permeability and state variables QoI, given the InSAR data information, insections 4.3 and 4.4.
The plausibility of the low-rank based Laplace approximation and the computa-tional feasibility of the gpCN method, discussed in sections 2.3 and 2.4, rely on theinherent low rank structure of the data misfit Hessian. This property allows for trun-cating the low-rank based Laplace approximation at r (cid:28) n n . In section 2.3 we provedthat this property holds for an idealized 1D poroelasticity problem. To demonstratethis property for the 3D Nevada test problem introduced in section 3 we computed theeigenvalues and eigendirections of the prior-preconditioned data misfit Hessian evalu-ated at the MAP point (Figure 2). The bottom right panel shows the decay of theeigenvalues on increasingly refined meshes. The decay patterns in the 16,896-nodeand 67,133-node meshes are almost identical, which demonstrates a mesh-independenteigenvalue decay on sufficiently resolved meshes. On all three meshes, only about r = 70–80 (cid:28) n n eigenvalues are larger than one, and therefore, only the r eigendirec-tions that correspond to these eigenvalues can be inferred from InSAR data with highconfidence. Selected eigendirections are shown in Figure 2. Similar to the 1D problemin section 2.3, the larger the eigenvalue the smoother the corresponding eigendirection,indicating that it is easier to infer smooth modes. The distinctive features of the firstfew eigendirections, which are inferred strongly from the data, are localized near thewell and the surrounding subsidence bowl, where surface deformation is largest. –12–anuscript submitted to Water Resources Research
Figure 2: Eigenvectors w i for i = 1 , , , , ,
76 of the generalized eigenproblem (15)using InSAR data and the 4 , λ i = 1. The eigenvalues corresponding tothe eigenvectors w i for i = 1 , , , , ,
76 in the 4 , ∼ λ i > For uncertainty quantification in realistic 3D transient multi-physics problems itis essential to improve the sampling of the parameter space. Incorporating log posteriorHessian information in gpCN proposal, through the low-rank based Laplace approxi-mation, leads to better convergence properties, compared to the pCN algorithm, seeTable 1 for methods performance. The gpCN algorithm achieves an acceptance rateof ∼
22% when using a step size β = 0 .
25. To achieve a similar acceptance rate forthe pCN algorithm, ∼ β = 0 .
01 is required, which leadsto smaller MCMC steps and slower exploration of the feasible parameter space. Con-sequently, the sample autocorrelation in gpCN chains decays much faster than in pCNchains, as shown in the bottom right panel in Figure 3. The large step-size in gpCNchains increases the number of independent MCMC samples and leads to a larger ef-fective sample size compared to pCN chains (Table 1). Visualizing chains mixing atselected points in space (Figure 3) highlights the performance difference between thetwo methods. The pCN chains vary slowly, fail to converge and do not explore thefeasible parameter space while gpCN chains move more rapidly and exhibit more ofa “fuzzy worm” pattern (Bui-Thanh, 2012) in which correlation between samples isweaker. –13–anuscript submitted to
Water Resources Research i − . − . − . − . − . − . − . − . − . l og ( κ i ( x )) x i − . − . − . − . − . − . l og ( κ i ( x )) x i − − − − − − − l og ( κ i ( x )) x i − − − − − − − − l og ( κ i ( x )) x i − . − . − . − . − . − . − . − . − . l og ( κ i ( x )) x lag . . . . . . . au t o c o rr e l a t i on i − . − . − . − . − . − . l og ( κ i ( x )) x x x x x x x gpCN1 gpCN2 gpCN3 gpCN4 gpCN5 pCN1 pCN2 Figure 3: pCN and gpCN chains mixing at selected points in space x = ( − . , − . , . x = ( − . , − . , . x = ( − . , . , . x = (2 . , . , . x = (0 , − . , . x = (0 , − . , .
6) in km. (Bottom right panel) autocorrelation versus lag. Autocorre-lation values are reported as average over the domain. The two pCN chains and the fivegpCN chains are labeled pCN1, pCN2, and gpCN1–gpCN5, respectively. –14–anuscript submitted to
Water Resources Research
Table 1: pCN and gpCN performanceMethod a PDE solves per chain β Acceptance rate IAT b ESS/chain c pCN 2 50,000 0 . ∼
24% 9754 5gpCN 5 50,000 0 . ∼
22% 494 101 a Each chain is of length 50,000 samples. b IAT is the pointwise integrated autocorrelation time reported as average over thedomain and over the number of chains. c ESS =
Chain lengthIAT is the effective sample size.
In this section we present the solution of the Bayesian inverse problem throughthe variance, samples, and selected marginal distributions of the permeability posteriordistribution. In particular we study both the low-rank based Laplace approximationand the MCMC-based posterior distributions.
Figure 4: Two-dimensional horizontal slices of the permeability pointwise standard de-viation (in decimal logarithm) at mid-aquifer depth of (a) the prior distribution, (b) thelow-rank based Laplace approximation ( r = 80), and (c) the MCMC-based true posterior.In all the panels, x is the east direction and y is the north direction relative to the welllocation ( x, y ) = (0 ,
0) (marked by the red star).The posterior distribution incorporates information from InSAR data and theaquifer model. This significantly reduces the pointwise standard deviation (SD) ofthe permeability posterior distributions compared to the prior pointwise SD in mostof the domain. Figure 4 shows this is the case for both the low-rank based Laplaceapproximation and the MCMC-based posterior. The spatial variation of the SD in-ferred from the low-rank based Laplace approximation in panel 4 b is very similar to –15–anuscript submitted to Water Resources Research the SD of the MCMC-based true posterior distribution in panel 4 c . The MCMC-basedtrue posterior, however, infers a slightly larger area where permeability can be inferredwith high confidence. This is also reflected in the samples obtained from these threedistributions. Samples from the prior distribution exhibit varying random patterns(Figure 5 a ) whereas samples from the low-rank based Laplace approximation and theMCMC-based true posterior (Figure 5 b and 5 c ) show consistent features, e.g. a lowpermeability area northwest the well and a high permeability channel extending to-ward the south/southwest. In the latter two distributions, permeability patterns atthe southern part of the domain display more randomness than in the center and thenorth. This randomness is expected because the SD at the southern part is relativelylarger (Figures 4 b and 4 c ) and is likely a result of two factors, the southern part ofthe domain is far from the major part of the subsidence bowl, and the InSAR data weuse in the inversion does not cover that area of the domain.Figure 5: Two-dimensional horizontal slices of permeability samples (in decimal log-arithm) at mid-aquifer depth from (a) the prior distribution, (b) the low-rank basedposterior Laplace approximation ( r = 80), and (c) the true posterior distribution (MCMCsamples). –16–anuscript submitted to Water Resources Research
To illustrate the broad range of results we compare 1D marginals of the priordistribution to posterior distributions inferred using both the low-rank based Laplaceapproximation and MCMC at selected points in the domain ( x i for i = 1 , , ..., r = 80 and r = 150 modes are very similar. This demonstrates that most of the infor-mation from the data is contained in the first 80 modes and adding the next 70 modesin the low-rank based Laplace approximation provides only limited additional informa-tion. The marginals of the MCMC-based true posterior distribution at points x i areobtained by applying 1D Gaussian kernel density estimation (KDE) to the pointwisevalues of the MCMC samples at x i . These marginals confirm the observation from thepointwise SD values (Figure 4) that the posterior distribution is significantly narrowerthan the prior distribution. Moreover, the pointwise marginals reveal that the MCMC-based true posterior exhibits more certainty in the permeability inference compared tothe low-rank based Laplace approximation in most of the selected pointwise locations( x to x ).The certainty in the inferred permeability varies with location in the domain asshown in Figure 6. We notice, for example, inference with high certainty at the location x , which is located in the high permeability channel extending south from the welland within the subsidence bowl. On the other hand, the existence of a low permeabilityregion northwest of the well reduces the GW flux and increases the uncertainty in theparameters in the region beyond this flow barrier, for example at location x . Thiseffect was also observed by Hesse and Stadler (2014) in synthetic model problems andwe refer to it here as the shadowed region. The uncertainty in the permeability of theshadowed region is large, because in the absence of flow there is no surface deformationthat informs the model permeability in this region.Slight deviation of the MCMC-based true posterior from a Gaussian distributioncan be observed from the posterior two-dimensional marginals (see Figure 7). Themean of the marginals at { x , x } and { x , x } is different from the MAP point andthe distribution displays right-skewness of the marginal at { x , x } . On the wholethough, the deviations of the MCMC-based true posterior from Gaussian distributionare minor in this particular case.Figure 8 shows marginals over a line segment parallel to the x axes that ex-tends from the west to the east of the domain and passes by the well at point( x, y, z ) = (0 , , . b ), the posterior distributions characterize the permeability with much more cer-tainty especially in the middle part of the line segment near the well from ∼ − ∼ c and 8 d ). This is mainly because this interval coincides with significantdeformation signal or deformation gradient (Figure 8 a ). In particular, the sharp dropin the permeability just west to the well location is inferred with very high certaintyas evident by the narrow marginal distribution at that location. The location of thisfeature coincides with large gradient in the subsidence from almost no deformation justwest to the well to significant subsidence at and east to the well. This demonstratesthat large gradients in the deformation data inform the permeability characterizationstrongly. w i In the preceding discussion, we studied the permeability characterization location-wise in the physical domain. Uncertainty in the permeability characterization can beviewed from a different perspective, in the directions of the parameter space that aredetermined by the generalized eigenvalue problem (15) (Petra et al., 2014). To thatend, we compute the marginals of the MCMC-based posterior distribution in the first –17–anuscript submitted to
Water Resources Research − − − − − − − − − . . . . . . P D F − − − − − − − − − . . . . . − − − − − − − − − . . . . . . . . . P D F − − − − − − − . . . . . . . − − − − − − − − − ( κ )0 . . . . . . P D F − − − − − − − − − ( κ )0 . . . . . x x x x x x PriorLaplace r=80 Laplace r=150MCMC Prior meanMAP MCMC mean
Figure 6: One-dimensional marginals at selected points x i for i = 1 , , ..., r = 80 and r = 150), and the MCMC-based true posterior distribution. –18–anuscript submitted to Water Resources Research
Figure 7: Two-dimensional marginals at selected pairs of points in space { x i , x j } for i, j = 1 , , , , i (cid:54) = j (coordinates are provided in Figure 3) of the prior, the low-rank based Laplace approximation ( r = 80), and the MCMC-based true posterior. Thevalues of the contour lines of each distribution are c cl × max(Π) for c cl = 0 . , . , . , . –19–anuscript submitted to Water Resources Research − u I nS A R ( mm ) (a) LOS deformation − − − − − l og ( κ ) (b) Prior − − − − − l og ( κ ) (c) Laplace − − x (km) − − − − − l og ( κ ) (d) MCMC LOS deformation Posterior MAP MCMC samples mean Prior mean
Figure 8: Prior, low-rank based Laplace approximation, and true posterior marginalsover a line. (a) InSAR LOS deformation signal over a line segment l parallel to the x axes, extends from the west to the east of the domain, and passes by the well location atthe surface. Marginal distributions, over a line segment l parallel to the x axes, extendsfrom the west to the east of the domain, and passes by the well location at z = . l to generate the approximate marginal PDFs. The grayshade is proportional to the PDF value. Values of the estimated PDF beyond the shadedarea are negligible. The well location is marked with a gray dotted vertical line and theprior mean, the MAP point, and the MCMC samples mean are shown in blue, green anddotted red, respectively (in decimal logarithm). –20–anuscript submitted to Water Resources Research
Figure 9: Marginal distributions of the gpCN MCMC samples components a i in the first r = 80 eigendirections w i , for i = 1 , , ..., r . The PDFs are shifted by their means ¯ a i .(a) PDFs of a i (for better visualization, we display only every fourth PDF and truncatethe PDF plot at PDF = 4). (b) PDFs of a i are plotted vertically in gray scale. Thedarker the gray shade, the higher the PDF value. Values of the estimated PDF beyondthe shaded area are negligible (less than ∼ . r = 80 eigendirections w i (we show selected directions in Figure 2). We denote with a i , for i = 1 , , ..., r , the components of the j th MCMC sample ηηη ( j )post in the eigendi-rections w i , j = 0 , , ..., n s where n s = 250 ,
000 is the total number of gpCN MCMCsamples. We can write the sample ηηη ( j )post as a linear combination of the eigendirections w i as follows: ηηη ( j )post = n n (cid:88) i =1 a i w i (19)Using the Γ − -orthogonality of the directions w i , the components a i are given by a i = w Ti Γ − ηηη ( j )post . We use 1D Gaussian KDE to approximate the distribution ofeach component a i over the samples ηηη ( j )post for j = 0 , , ..., n s . We show a i distributionsin Figure 9, shifted by their averages ¯ a i for ease of comparison. We notice that, asexpected, the more dominant the eigendirection in the generalized eigenvalue prob-lem (15), the more certain the inference of the component in that direction. The firstfew directions, in particular, are inferred with significantly high confidence. We alsoobserve that the distributions of some a i ’s deviate slightly from a Gaussian distributionshowing a mild skewness. Quantifying the propagation of uncertainty in the parameters to the state vari-ables QoIs, e.g. pore pressure and total subsidence, is of interest in GW managementapplications (Qin et al., 2018). We show the prior and posterior predictive distribu-tions of the pressure p ( x i , t InSAR ), at the well location x = (0 , , . x , x , x and x (coordinates are provided in Figure 3) in Figure 10. We also showthese predictive distributions for total subsidence volume − (cid:82) Γ top u ( t InSAR ) ds , whereΓ top is the domain top surface (bottom panel in Figure 10). We use KDE to esti-mate the PDF of each of those QoIs. We note that these distributions can be highlynon-Gaussian and that the prior-based distributions vary across a much wider rangecompared to the narrow posterior-based distributions. In some locations, x , x and x for example, relying only on prior knowledge significantly underestimates the expected –21–anuscript submitted to Water Resources Research pressure drop; the posterior-based QoI mean at these locations is considerably largerin magnitude than that of the prior-based QoI. We also note appreciable differencesbetween the low-rank based Laplace approximation and the MCMC based posteriorpredictives. The mixing behaviour of the chains for these QoIs distributions (shownin Figure 11) is similar to what was observed for the parameter chains in section 4.2.
The use of InSAR surface deformation data to quantify the uncertainty in GWaquifer models via solution of poroelasticity-governed high dimensional Bayesian in-verse problems is important for data-driven, model-predictive management of GW re-sources. The intrinsic low dimensionality of this inverse problem, evident by the rapiddecay in data misfit Hessian eigenfunctions, makes solving this problem computation-ally tractable when the right computational tools are used, e.g. our proposed gpCNMCMC method employing low-rank-based Laplace approximation proposals, scalablepriors based on inverses of differential operators, and an inexact Newton method withadjoint-based gradients and Hessian actions to find the MAP point. The use of thesetools guarantees that the number of forward poroelasticity PDE solves required tosolve the inverse problem is independent of the parameter dimension. In the Nevadatest case, for example, only r ≈
80 (out of n n = 4 , –22–anuscript submitted to Water Resources Research − . − . . Pressure (Pa) × . . . . . P D F PriorLaplaceMCMC Prior meanLaplace meanMCMC mean − − Pressure (Pa) . . . . . P D F − − Pressure (Pa) . . . . . P D F − − − Pressure (Pa) . . . . . P D F − − − Pressure (Pa) . . . . . P D F − − − − Subsidence volume ( m ) . . . . P D F . . . . P r i o r P D F × − . . . . . P r i o r P D F . . . . P r i o r P D F . . . . P r i o r P D F . . . . . P r i o r P D F . . . . . P r i o r P D F x x x x x Figure 10: Prior, low-rank based Laplace approximation, and MCMC-based posteriorpredictives for pointwise pressure values at selected locations ( x , x , x , x , and x )and total subsidence volume (bottom panel). In all panels, the values of the Laplaceapproximation-based and MCMC-based PDFs are measured by the left y -axis scale whilethe right y -axis scale measures the prior PDF. –23–anuscript submitted to Water Resources Research i − − − − − − − − P r e ss u r ea t x ( P a ) x i − − − − − − P r e ss u r ea t x ( P a ) x i − − − − − − P r e ss u r ea t x ( P a ) x i − − − − − − − P r e ss u r ea t x ( P a ) x i − − − − − − − − P r e ss u r ea t x ( P a ) x i − − − − − − S ub s i den c e v o l u m e ( m ) gpCN1 gpCN2 gpCN3 gpCN4 gpCN5 pCN1 pCN2 Figure 11: pCN and gpCN MCMC-based mixing of the pointwise pressure values at se-lected locations ( x , x , x , x , and x ) and total subsidence volume (bottom right panel).The two pCN chains and the five gpCN chains are labeled pCN1, pCN2, and gpCN1–gpCN5, respectively. –24–anuscript submitted to Water Resources Research
Appendix A Adjoint-Based Derivation of the Hessian Operator
We derive the Hessian operator in the Newton iteration that we form for findingthe MAP point. Maximizing the posterior distribution (7) with respect to the param-eter m to find the MAP point is equivalent to minimizing the negative log posterior,i.e. min m ∈ IR nn J ( m ) = 12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs ) + 12 ( m − ¯ m ) T Γ − ( m − ¯ m ) , (A1)where ¯X depends on m through the solution of the forward problem ¯S ( m ) ¯X = ¯F .To derive the gradient using the adjoint method, we form the Lagrangian for theconstrained optimization problem (A1) to be: L ( m , ¯X , ¯Y ) = 12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs )+ 12 ( m − ¯ m ) T Γ − ( m − ¯ m ) + ¯Y T ( ¯S ( m ) ¯X − ¯F ) . (A2)where ¯Y is the Lagrange multiplier (also called the adjoint variable).At a minimum of (A2), the partial derivatives of the Lagrangian with respect tothe state variable, ∂ L /∂ ¯X , the parameter, ∂ L /∂ m , and the adjoint variable, ∂ L /∂ ¯Y ,vanish. We set ∂ L /∂ ¯X and ∂ L /∂ ¯Y to zero, which gives the state equation and theadjoint equation respectively: ∂ L ∂ ¯Y ( ¯X , ¯Y , m ) = ¯S ( m ) ¯X − ¯F = 0 (state) , (A3) ∂ L ∂ ¯X ( ¯X , ¯Y , m ) = ¯B T Γ − ( ¯B ¯X − d obs ) + ¯S T ( m ) ¯Y = 0 (adjoint) . (A4)We satisfy these two conditions by solving the discretized state and adjoint equationsexactly for the given value of m . Thus we seek the parameter m that ensures that thegradient ∂ L ∂ m ( ¯X , ¯Y , m ) = Γ − ( m − ¯ m ) + ¯C T ( m ) ¯Y (A5)vanishes. The operator ¯C is the partial derivative (sensitivity) of the residual R = ¯S ¯X − ¯F with respect to the parameter m¯C = ∂∂ m ( ¯S ¯X − ¯F ) = ∂∂ m ( ¯S ¯X ) . (A6)We solve the equation ∂ L /∂ m = 0 using Newton’s method.To derive the Newton iteration we form the Lagrangian L H : L H := δ ¯X T (cid:0) ¯B T Γ − ( ¯B ¯X − d obs ) + ¯S T ( m ) ¯Y (cid:1) + δ ¯Y T (cid:0) ¯S ( m ) ¯X − ¯F (cid:1) + δ m T (cid:0) Γ − ( m − ¯ m ) + ¯C T ( m ) ¯Y (cid:1) , (A7)where δ ¯X and δ ¯Y are Lagrange multipliers for the adjoint (A4) and state (A3) equa-tions, and δ m is the direction in which the Hessian acts. The so-called incrementalforward and incremental adjoint problems are defined as follows: ∂ L H ∂ ¯Y = ¯S ( m ) δ ¯X + ¯C ( m ) δ m = 0 . (incremental forward) , (A8) ∂ L H ∂ ¯X = ¯B T Γ − ¯B δ ¯X + ¯S T ( m ) δ ¯Y + (cid:16) ∂∂ ¯X ( ¯C T ( m ) ¯Y ) (cid:17) T δ m = W ¯X ¯X δ ¯X + ¯S T ( m ) δ ¯Y + W ¯Xm ( m ) δ m = 0 . (incremental adjoint) (A9) –25–anuscript submitted to Water Resources Research
The Hessian action in the direction δ m is given by: ∂ L H ∂ m = ∂∂ m ( δ ¯X T ¯S T ( m ) ¯Y ) + ∂∂ m ( δ ¯Y T ¯S ( m ) ¯X ) + Γ − δ m + ∂∂ m ( δ m T ¯C T ( m ) ¯Y )= W m ¯X ( m ) δ ¯X + ¯C T ( m ) δ ¯Y + Γ − δ m + R mm ( m ) δ m . (A10)Substituting the solutions for δ ¯X and δ ¯Y from the systems (A8) and (A9) into theexpression (A10) gives the Hessian expression: H ( m ) = Γ − + H misfit ( m ) , (A11)where H misfit ( m ) := R mm ( m ) + ¯C T ( m ) ¯S − T ( m )( W ¯X ¯X ¯S − ( m ) ¯C ( m ) − W ¯Xm ( m )) − W m ¯X ( m ) ¯S − ( m ) ¯C ( m ) . (A12)The matrices W ¯X ¯X , W ¯Xm , W m ¯X , and R mm are defined as follows W ¯X ¯X = ¯B T Γ − ¯B , W ¯Xm = ∂∂ m ( ¯S T ¯Y ) , W m ¯X = ∂∂ ¯X ( ¯C T ¯Y ) , R mm = ∂∂ m ( ¯C T ¯Y ) . Forming the Hessian, evaluated at a given m , explicitly, requires applying theoperator ¯S − to each column of ¯C to form the matrices product ¯S − ¯C . Applyingthe operator ¯S − to a vector amounts to an incremental forward poroelasticity solve,equation (A8). Therefore, forming the GN Hessian explicitly costs n n poroelasticityPDEs solves (note that no additional solves are required for forming the product ¯C T ¯S − T since it is the transpose of the product ¯S − ¯C ). However, applying H to adirection mainly costs applying the operators ¯S − and ¯S − T each to a vector whichis the cost of two poroelasticity PDEs solves (one incremental forward (A8) and oneincremental adjoint (A9)). Appendix B Forward Model Parameters –26–anuscript submitted to
Water Resources Research
Table B1: Forward model parameters (a) Model parameters
Volumetric force ( f u ) 0Pumping rate ( − (cid:82) Ω f p d x ) 9 ,
028 m / dayFluid viscosity ( µ ) 0 .
001 Pa · sBiot-Willis coefficient ( α ) 0 . h ) 885 mWater density ( ρ ) 997 .
97 Kg / m Gravitational acceleration ( g ) 9 . / s Water compressibility ( β w ) 4 . × − Pa − (b) Layer parameters Units Aquifer Confining layer Lower clayPoisson’s ratio ( ν ) - 0 .
25 0 .
25 0 . G ) Pa 3 . × . × × Specific storage ( S (cid:15) ) Pa − . × − . × − . × − Prior mean value (log ( κ )) m (for κ ) -11.2 -14.2 -14.2 Acknowledgments
This work was supported by National Science Foundation (NSF) Grants CBET–1508713and ACI-1550593, and DOE grant DE-SC0019303. A.A. would like to acknowledgefunding from the Ministry of Education in Saudi Arabia. The Envisat SAR imageryused in this study can be downloaded through the UNAVCO Data Center SAR archive.
References
Alexanderian, A., Gloor, P. J., & Ghattas, O. (2016). On Bayesian A-and D-optimalexperimental designs in infinite dimensions.
Bayesian Analysis , (3), 671–695. doi: 10.1214/15-BA969Alghamdi, A. (2020). Bayesian inverse problems for quasi-static poroelasticitywith application to ground water aquifer characterization from geodetic data (Unpublished doctoral dissertation).Alghamdi, A., Hesse, M. A., Chen, J., & Ghattas, O. (2020). Bayesian poroe-lastic aquifer characterization from InSAR surface deformation data. PartI: Maximum a posteriori estimate.
Water Resources Research , (10),e2020WR027391. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020WR027391 (e2020WR02739110.1029/2020WR027391) doi: 10.1029/2020WR027391Beskos, A., Girolami, M., Lan, S., Farrell, P. E., & Stuart, A. M. (2017). Geomet-ric MCMC for infinite-dimensional inverse problems. Journal of ComputationalPhysics , , 327-351.Biot, M. A. (1941). General Theory of Three- Dimensional Consolidation. Journal ofApplied Physics . doi: 10.1063/1.1712886 –27–anuscript submitted to
Water Resources Research
Bohling, G. C., & Butler, J. J. (2010). Inherent Limitations of Hydraulic Tomogra-phy.
Ground Water , (6), 809–824. doi: 10.1111/j.1745-6584.2010.00757.xBui-Thanh, T. (2012, May). A gentle tutorial on statistical inversion using thebayesian paradigm (Technical report No. ICES-12-18). Institute for Computa-tional Engineering and Sciences. (ICES Report)Bui-Thanh, T., Burstedde, C., Ghattas, O., Martin, J., Stadler, G., & Wilcox, L. C.(2012). Extreme-scale UQ for Bayesian inverse problems governed by PDEs.In
Sc12: Proceedings of the international conference for high performancecomputing, networking, storage and analysis.
Bui-Thanh, T., Ghattas, O., Martin, J., & Stadler, G. (2013). A computationalframework for infinite-dimensional Bayesian inverse problems Part I: The lin-earized case, with application to global seismic inversion.
SIAM Journal onScientific Computing , (6), A2494-A2523. doi: 10.1137/12089586XBurbey, T. J. (2006). Three-dimensional deformation and strain induced by munici-pal pumping, part 2: Numerical analysis. Journal of Hydrology , (3-4), 422 -434.Burbey, T. J. (2008). The influence of geologic structures on deformation due toground water withdrawal. Ground water , (2), 202–211.Burbey, T. J., Warner, S. M., Blewitt, G., Bell, J. W., & Hill, E. (2006). Three-dimensional deformation and strain induced by municipal pumping, part 1:Analysis of field data. Journal of Hydrology , (1-4), 123 - 142.Carrera, J., & Neuman, S. P. (1986a). Estimation of Aquifer Parameters UnderTransient and Steady State Conditions: 1. Maximum Likelihood Method Incor-porating Prior Information. Water Resources Research , (2), 199–210. doi:10.1029/WR022i002p00199Carrera, J., & Neuman, S. P. (1986b, 2). Estimation of Aquifer Parameters Un-der Transient and Steady State Conditions: 2. Uniqueness, Stability, andSolution Algorithms. Water Resources Research , (2), 211–227. Re-trieved from http://doi.wiley.com/10.1029/WR022i002p00211 doi:10.1029/WR022i002p00211Carrera, J., & Neuman, S. P. (1986c). Estimation of Aquifer ParametersUnder Transient and Steady State Conditions: 3. Application to Syn-thetic and Field Data. Water Resources Research , (2), 228–242. doi:10.1029/WR022i002p00228Cotter, S. L., Roberts, G. O., Stuart, A. M., & White, D. (2012). MCMC methodsfor functions: modifying old algorithms to make them faster.(submitted)Cui, T., Law, K., & Marzouk, Y. (2016). Dimension-independent likelihood-informedMCMC. Journal of Computational Physics , , 109-137.Eaton, T. T. (2006, 2). On the importance of geological heterogeneity for flowsimulation. Sedimentary Geology , (3-4), 187–201. Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0037073805003714 doi: 10.1016/j.sedgeo.2005.11.002Famiglietti, J. S. (2014). The global groundwater crisis. Nature Climate Change , (11), 945.Ferretti, A., Savio, G., Barzaghi, R., Borghi, A., Musazzi, S., Novali, F., . . . Rocca,F. (2007). Submillimeter accuracy of insar time series: Experimental vali-dation. IEEE Transactions on Geoscience and Remote Sensing , (5), 1142-1153.Ferronato, M., Castelletto, N., & Gambolati, G. (2010). A fully coupled 3-D mixedfinite element model of Biot consolidation. Journal of Computational Physics , (12), 4813 - 4830. Retrieved from doi: 10.1016/j.jcp.2010.03.018Flath, P. H., Wilcox, L. C., Ak¸celik, V., Hill, J., van Bloemen Waanders, B., &Ghattas, O. (2011). Fast algorithms for Bayesian uncertainty quantification –28–anuscript submitted to Water Resources Research in large-scale linear inverse problems based on low-rank partial Hessian ap-proximations.
SIAM Journal on Scientific Computing , (1), 407-432. doi:10.1137/090780717Haga, J. B., Osnes, H., & Langtangen, H. P. (2012). Biot’s consolidation, pressureoscillations, elastic locking, low-permeable media, finite elements. Interna-tional Journal for Numerical and Analytical Methods in Geomechanics , (12),1507–1522.Hesse, M., & Stadler, G. (2014). Joint inversion in coupled quasistatic poroelasticity. Journal of Geophysical Research: Solid Earth , (2), 1425–1445.Holzer, T. L., & Galloway, D. L. (2005). Impacts of land subsidence caused by with-drawal of underground fluids in the united states. Humans as geologic agents , , 87.Isaac, T., Petra, N., Stadler, G., & Ghattas, O. (2015, September). Scalable andefficient algorithms for the propagation of uncertainty from data through in-ference to prediction for large-scale problems, with application to flow of theAntarctic ice sheet. Journal of Computational Physics , , 348-368. doi:10.1016/j.jcp.2015.04.047Konikow, L. F., & Kendy, E. (2005). Groundwater depletion: A global problem. Hy-drogeology Journal , (1), 317–320.Linde, N., Ginsbourger, D., Irving, J., Nobile, F., & Doucet, A. (2017). On uncer-tainty quantification in hydrogeology and hydrogeophysics. Advances in WaterResources , , 166–181.Lindgren, F., Rue, H., & Lindstr¨om, J. (2011). An explicit link between Gaussianfields and Gaussian Markov random fields: the stochastic partial differentialequation approach. Journal of the Royal Statistical Society: Series B (Sta-tistical Methodology) , (4), 423–498. Retrieved from http://dx.doi.org/10.1111/j.1467-9868.2011.00777.x doi: 10.1111/j.1467-9868.2011.00777.xLogg, A., Mardal, K.-A., & Wells, G. N. (Eds.). (2012). Automated solution of dif-ferential equations by the finite element method (Vol. 84). Springer. doi: 10.1007/978-3-642-23099-8MacKay, D. J. (2003).
Information theory, inference and learning algorithms . Cam-bridge university press.Martin, J., Wilcox, L. C., Burstedde, C., & Ghattas, O. (2012). A stochastic New-ton MCMC method for large-scale statistical inverse problems with applicationto seismic inversion.
SIAM Journal on Scientific Computing , (3), A1460-A1487. doi: 10.1137/110845598McLaughlin, D., & Townley, L. R. (1996). A reassessment of the groundwater in-verse problem. Water Resources Research , (5), 1131–1161. doi: 10.1029/96WR00160Oliver, D. S., & Chen, Y. (2011, 1). Recent progress on reservoir history matching:a review. Computational Geosciences , (1), 185–221. Retrieved from http://link.springer.com/10.1007/s10596-010-9194-2 doi: 10.1007/s10596-010-9194-2Petra, N., Martin, J., Stadler, G., & Ghattas, O. (2014). A computational frame-work for infinite-dimensional Bayesian inverse problems: Part II. StochasticNewton MCMC with application to ice sheet inverse problems. SIAM Journalon Scientific Computing , (4), A1525–A1555.Phillips, P. J., & Wheeler, M. F. (2009). Overcoming the problem of locking inlinear elasticity and poroelasticity: an heuristic approach. Computational Geo-sciences , (1), 5–12.Pinski, F. J., Simpson, G., Stuart, A. M., & Weber, H. (2015). Algorithms forKullback–Leibler approximation of probability measures in infinite dimensions. SIAM Journal on Scientific Computing , (6), A2733–A2757.Qin, H., Andrews, C. B., Tian, F., Cao, G., Luo, Y., Liu, J., & Zheng, C. (2018).Groundwater-pumping optimization for land-subsidence control in beijing –29–anuscript submitted to Water Resources Research plain, china.
Hydrogeology Journal , (4), 1061–1081.Tarantola, A., Valette, B., et al. (1982). Inverse problems = quest for information. Journal of geophysics , (1), 159–170.Tom´as, R., Romero, R., Mulas, J., Marturi`a, J. J., Mallorqu´ı, J. J., Lopez-Sanchez, J. M., . . . Blanco, P. (2014). Radar interferometry techniquesfor the study of ground subsidence phenomena: a review of practical is-sues through cases in spain. Environmental Earth Sciences , (1), 163–181. Retrieved from https://doi.org/10.1007/s12665-013-2422-z doi:10.1007/s12665-013-2422-zVilla, U., Petra, N., & Ghattas, O. (2018). hIPPYlib: An extensible software frame-work for large-scale deterministic and Bayesian inversion. Journal of OpenSource Software , (30), 940.Villa, U., Petra, N., & Ghattas, O. (2020). hIPPYlib: An extensible software frame-work for large-scale inverse problems governed by PDEs; Part I: Deterministicinversion and linearized bayesian inference. Transactions on MathematicalSoftware . (in print)Wada, Y., Van Beek, L., & Bierkens, M. F. (2012). Nonsustainable groundwater sus-taining irrigation: A global assessment.
Water Resources Research , (6).Wada, Y., Van Beek, L. P., Van Kempen, C. M., Reckman, J. W., Vasak, S., &Bierkens, M. F. (2010). Global depletion of groundwater resources. Geophysi-cal research letters , (20).Wada, Y., Wisser, D., & Bierkens, M. F. (2014). Global modeling of withdrawal,allocation and consumptive use of surface water and groundwater resources. Earth System Dynamics Discussions , (1), 15–40.Xue, Y.-Q., Zhang, Y., Ye, S.-J., Wu, J.-C., & Li, Q.-F. (2005). Land subsidence inchina. Environmental geology , (6), 713–720.(6), 713–720.