Bayesian Poroelastic Aquifer Characterization from InSAR Surface Deformation Data. Part I: Maximum A Posteriori Estimate
BBayesian Poroelastic Aquifer Characterization fromInSAR Surface Deformation Data.Part I: Maximum A Posteriori Estimate
Amal Alghamdi , Marc A. Hesse , , Jingyi Chen , , Omar Ghattas , , 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 University of Texas at Austin, Mechanical Engineering, Austin, TX, United States
Key Points: • A Bayesian framework was developed for characterizing heterogeneous aquifer prop-erties from surface deformation data. • Scalable algorithms allow the inference of high-dimensional discretized parame-ters in a transient fully coupled 3D poroelastic aquifer model. • InSAR data significantly improves characterization of lateral aquifer heterogene-ity at a test site in Nevada.
Corresponding author: Marc Hesse, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] F e b bstract Characterizing the properties of groundwater aquifers is essential for predicting aquiferresponse and managing groundwater resources. In this work, we develop a high-dimensionalscalable Bayesian inversion framework governed by a three-dimensional quasi-static lin-ear poroelastic model to characterize lateral permeability variations in groundwater aquifers.We determine the maximum a posteriori (MAP) point of the posterior permeability dis-tribution from centimeter-level surface deformation measurements obtained from Inter-ferometric Synthetic Aperture Radar (InSAR). The scalability of our method to high pa-rameter dimension is achieved through the use of adjoint-based derivatives, inexact New-ton methods to determine the MAP point, and a Mat´ern class sparse prior precision op-erator. Together, these guarantee that the MAP point is found at a cost, measured innumber of forward/adjoint poroelasticity solves, that is independent of the parameterdimension. We apply our methodology to a test case for a municipal well in Mesquite,Nevada, in which InSAR and GPS surface deformation data are available. We solve prob-lems with up to 320 ,
824 state variable degrees of freedom (DOFs) and 16 ,
896 param-eter DOFs. A consistent treatment of noise level is employed so that the aquifer char-acterization result does not depend on the pixel spacing of surface deformation data. Ourresults show that the use of InSAR data significantly improves characterization of lat-eral aquifer heterogeneity, and the InSAR-based aquifer characterization recovers com-plex lateral displacement trends observed by independent daily GPS measurements.
The sustainable management of groundwater (GW) resources is quickly emergingas a critical issue as irrigated agriculture continues to grow and rely more and more onGW (Hoekstra & Mekonnen, 2012). Estimates suggest that half the GW extracted forirrigation exceeds aquifer recharge and is hence unsustainable (Rost et al., 2008; Voros-marty, 2000). In particular, areas with persistent water stress and large aquifer systemshave increased their reliance on GW and now experience sustained GW depletion (Hanasakiet al., 2008; Wada et al., 2010; Scanlon et al., 2012), which poses a threat to agriculturalsystems, food security, and livelihoods (UNESCO, 2012).Sustainable GW management under changing climatic and societal conditions re-quires models that can predict the aquifer response to these forcings. The dominant un-certainty in the predictions of any GW model is the insufficient characterization of theaquifer properties (Eaton, 2006; Bohling & Butler, 2010; Oliver & Chen, 2011). Theseproperties can change by several orders of magnitude and display variability on all scales.In principle, the uncertainty in these parameters could be quantified with Bayesian in-ference (Carrera & Neuman, 1986a, 1986b; Yeh, 1986; McLaughlin & Townley, 1996),if a sufficiently rich data set is available. However, standard aquifer characterization basedon sparse measurements at wells is often not sufficient to significantly reduce the uncer-tainty in model predictions (Bohling & Butler, 2010). Hence there is an urgent need toimprove GW monitoring and integrate new data into GW models to reduce the uncer-tainty in predictions and to facilitate decision making.Over the past several decades, advances in satellite remote sensing techniques haverevolutionized capabilities for observing freshwater resources. In particular, spaceborneInterferometric Synthetic Aperture Radar (InSAR) has been measuring surface defor-mation history since 1992 with 10–100 meter spatial resolution and up to millimeter-levelmeasurement accuracy (Hanssen, 2001; Rosen et al., 2000). The potential use of surfacedeformations derived from InSAR in aquifer characterization has been recognized sinceas early as 2000; see D. Galloway and Hoffmann (2007) for a review. Since then space-borne InSAR missions such as ERS-1/2, Envisat, ALOS, and Sentinel-1 have generatedlarge amounts of freely available data (1992–present) for monitoring the state of confined –2– roundwater aquifers. Unlike data types that require field work and instrumentation, In-SAR data can be obtained and processed at relatively low cost and effort.However, despite significant work on InSAR based aquifer monitoring (Amelunget al., 1999; Hoffmann et al., 2001; Du & Olson, 2001; Lu & Danskin, 2001; D. Vasco etal., 2001; Schmidt & Br¨ugmann, 2003; D. Galloway & Hoffmann, 2007; D. W. Vasco etal., 2008; Burbey, 2008; Bell, 2008; D. L. Galloway & Burbey, 2011; Chaussard et al.,2014; J. Chen et al., 2016, 2017; Miller et al., 2017), the potential for this data for in-forming lateral aquifer heterogeneity has not yet been realized. This requires scalablecomputational approaches that allow the inference of high-dimensional parameter fieldsdescribing the aquifer properties. It also requires the use of geomechanical models thatlink the observed surface deformation to GW flow in the aquifer. Such poroelastic in-version methods have only recently been developed by Iglesias and McLaughlin (2012)and Hesse and Stadler (2014) and have not yet been applied to real field data.Here we build on the work of Hesse and Stadler (2014) and develop a high-dimensionalscalable poroelastic Bayesian inversion framework to characterize lateral permeabilityvariations in groundwater aquifers based on InSAR surface deformation data. The Bayesianframework provides a systematic and rational framework for quantifying uncertaintiesin the inverse solution, stemming from uncertainties in the data, model, and any priorinformation on the parameters, along with insensitivity of observables to parameters (Tarantola,2005; Kaipio & Somersalo, 2005; Stuart, 2010). We apply this framework to a well-characterizedfield site in Nevada (Burbey, 2006) and focus on the inference of lateral permeability vari-ations. In particular, we find the maximum a posteriori (MAP) point of the Bayesianposterior distribution. This application is computationally challenging and requires theinference of 16 ,
896 parameters in a three-dimensional geomechanical model with 320 , The solution of realistic aquifer characterization problems, such as the test case insection 3, requires robust, efficient, and scalable methods for both the forward and in-verse problems. The Bayesian inversion framework provides a means of incorporatingprior assumptions about the aquifer properties and the data noise probability distribu-tion to infer model parameters that comply with these assumptions and with a forwardmodel. Section 2.1 describes the formulation and a robust discretization of the poroe-lastic forward problem and section 2.2 outlines the application of the Bayesian inversionframework to this forward problem with InSAR surface deformation data.
Here we use quasi-static linear poroelasticity to model coupled groundwater flowand elastic deformation in a confined aquifer (Biot, 1941; Showalter, 2000; Wang, 2000;Segall et al., 2010). We adopt a three-field formulation and a mixed discretization to con-serve mass discretely and reduce numerical oscillations (Phillips, 2005; Phillips & Wheeler, –3– × (0 , T ] can be written as:( S (cid:15) p + α ∇ · u ) t + ∇ · q = f p (1) −∇ · ( σσσ ( u ) − αp I ) = f u (2) q + e m µ ∇ p =0 , (3)where ( · ) t denotes time derivative, p = p ( x , t ) is the deviation from hydrostatic pres-sure, q ( x , t ) is the volumetric fluid flux, u = u ( x , t ) is the displacement of the solid skele-ton, f p ( x , t ) is a fluid source per unit volume, f u ( x , t ) is the body force per unit volume, S (cid:15) is the specific storage, κ ( x ) = e m ( x ) is the medium permeability field, µ is the porefluid dynamic viscosity, α is the Biot–Willis coupling parameter, and σσσ ( u ) is the stresstensor. The medium permeability is written in terms of the log-permeability field m toensure that the inferred permeability field is positive. The initial and boundary condi-tions for the state variables u , p , and q are given by: p ( x ,
0) = p in Ω q · n = g on (0 , T ] × ∂ Ω np p = p d on (0 , T ] × ∂ Ω dp ( σσσ ( u ) − αp I ) n = g on (0 , T ] × ∂ Ω n u u = u d on (0 , T ] × ∂ Ω d u , (4)where p = p ( x ) is the initial pressure at time t = 0, g ( x ) is the value of the fluidflux across the pressure Neumann boundaries (0 , T ] × ∂ Ω np , p d ( x ) is the prescribed pres-sure on the pressure Dirichlet boundaries (0 , T ] × ∂ Ω dp , g ( x ) is the prescribed tractionon the displacement Neumann boundaries (0 , T ] × ∂ Ω n u , and u d ( x ) is the prescribed dis-placement on the displacement Dirichlet boundaries (0 , T ] × ∂ Ω d u . The stress tensor interms of displacement is given for isotropic linear elasticity by σσσ ( u ) = G ( ∇ u + ∇ u T ) + 2 Gν − ν ( ∇ · u ) I , (5)where G is the drained shear modulus and ν is the Poisson’s ratio.We discretize the system (1)–(3) with the Mixed Finite Element Method (MFEM)proposed by Ferronato et al. (2010), which approximates the fluid flux in the lowest-orderRaviart–Thomas space, the pressure in the space of piecewise constant functions, andthe displacement in the first-order Lagrange polynomial space. We use the implicit Eu-ler method for integration in time. In Appendix A, we provide the details of the weakform and discretization. The resulting linear system of equations that needs to be solvedat each time step is given by: MP k + α DU k + ∆ t k D q Q k =∆ t k F pk + MP k − + α DU k − − α GP k − EU k = − F u k − F u ,bcsk ∆ t k G p P k − ∆ t k KQ k =∆ t k F p,bcsk , (6)where P k , U k and Q k are the discretized pressure, displacement and fluid flux degreesof freedom (DOFs) vectors, respectively. The subscript k denotes the evaluation of thevariable (or the operator) at time t k and ∆ t k = t k − t k − is the k -th time step size.The discrete differential operators are denoted D and D q for the divergence of the dis-placement and the flux, G and G p for the gradient of the pressure in the pressure andthe mixed flow spaces, respectively, and E for the linear elastic operator. The mass ma-trices M and K are weighted by the specific storage and mobility (permeability/viscosity),respectively. Right hand side contributions include the discrete source terms for fluid andbody force, F pk and F u k , as well as boundary tractions, F u ,bcsk , and pressure natural bound-ary conditions, F p,bcsk . –4– o facilitate the derivation of the time-dependent inverse problem in section 2.2.4,we follow Hesse and Stadler (2014) and define the “global” space-time system that rep-resents the solution for all time steps simultaneously and denote it as ¯S ¯X = ¯F , (7)where the block lower triangular matrix ¯S combines the differential operators of system (6), ¯F combines the source terms and the natural boundary conditions, and ¯X consists of thepressure, displacement, and fluid flux DOFs at all time steps (see Appendix A for theexplicit expressions). The bar notation is used to distinguish the space-time discretizedoperators and vectors from the space-only discretized operators and vectors. The Bayesian inverse problem solution is a probabilistic characterization of unknownphysical system parameters, such as material properties, source terms, boundary con-ditions, and initial conditions. Bayes’ theory provides a framework for updating “prior”statistical knowledge about these parameters by using observational data and a math-ematical model of the system via a “likelihood distribution”. This likelihood distribu-tion determines how likely it is that the observed data result from a particular param-eter realization. The updated probability distribution is known as the posterior distri-bution, and is regarded as the solution of the inverse problem. Bayesian inversion, there-fore, provides a characterization of the uncertainty in the parameters due to uncertaintyin the prior information, data, and model, as opposed to finding a “point estimate”, asis the case with deterministic inversion. In the present article (Part I), we formulate theBayesian inversion problem, paying special attention to the prior distribution of the modelparameters and the observational noise model. We also take the first step in exploringthe posterior distribution by introducing fast Hessian-based methods for scalably deter-mining the maximum a posteriori (MAP) point. In a subsequent article (Part II), we willbuild on the tools developed here to construct a Hessian-driven Markov chain Monte Carlomethod (MCMC) to sample the posterior distribution and compute statistics of inter-est as well as posterior predictives.In section 2.2.1, we describe the formulation and discretization of the Bayesian in-verse problem for the three-field formulation of the linear quasi-static poroelasticity modeldescribed in section 2.1. Our methodology follows closely the framework presented byHesse and Stadler (2014), with the following differences or extensions: (1) we base theframework on the three-field Biot system formulation rather than the classic two-field,(2) we use a Mat´ern class prior for the Bayesian inverse problem (section 2.2.2), and (3)we define the likelihood distribution based on a noise model for InSAR surface deforma-tion data (section 2.2.3).
Given the observational data, d obs ∈ IR n obs , and prior statistical knowledge abouta parameter field m , we seek the posterior probability distribution of m . The parame-ter m is related to the observational data through a Gaussian additive noise model, d obs = F ( m ) + η. (8)The parameter-to-observable map , F ( · ) : Z → IR n obs , maps a realization of the pa-rameter field from the infinite dimensional space Z to the finite dimensional space of ob-servables IR n obs . Evaluating F ( m ) involves solution of a system of PDEs in which theparameter m is a coefficient, source term, boundary condition, or initial condition of thesystem, or a combination thereof. In our case, we will invert for the log-permeability field –5– appearing in equation (3). The parameter-to-observable map F ( m ) involves solutionof the poroelasticity model (1)–(3), and extraction of surface displacements.The random variable η is observational noise with known statistical properties. Here, η is assumed to follow a Gaussian distribution with mean zero and covariance matrix Γ noise ∈ IR n obs × n obs , reflecting the estimated noise level and correlation. The likelihood proba-bility measure, which is the probability measure of the difference between the observa-tional data and the simulated observations, d obs − F ( m ), follows the normal distribu-tion of the observational data noise.A general framework for infinite dimensional Bayesian inverse problems governedby PDEs is presented in Stuart (2010). In infinite dimensions, Bayes’ formula reads dµ post dµ prior = 1 Z π likelihood ( d obs | m ) , (9)where π likelihood is the likelihood probability density function, Z = (cid:82) Ω π likelihood ( d obs | m ) dµ prior is the normalization constant required for the posterior measure, µ post , to be a proba-bility measure, and dµ post dµ prior is the Radon–Nikodym derivative of the posterior with respectto the prior.As a first step to characterising the posterior µ post , we wish to find the point thatmaximizes the posterior, the so-called MAP point. We follow the discretize-then-optimize approach because it guarantees a consistent discretized gradient (Gunzburger, 2003). Thecorresponding Bayes formula in finite dimensions is given by: π post ( m | d obs ) ∝ π likelihood ( d obs | m ) π prior ( m ) , (10)where the symbol ∝ indicates equality up to a constant, m ∈ IR n n denotes the discretizedparameters, π prior ( m ) = N ( ¯ m , Γ prior ) , and π likelihood ( d obs | m ) ∝ e − (cid:96) ( d obs ; m ) . (11)Here N is the normal distribution of the prior with mean ¯ m ∈ IR n n and covariance ma-trix Γ prior ∈ IR n n × n n . The function (cid:96) , the negative log of the likelihood, is obtained fromequation (8) and the fact that η ∼ N (0 , Γ noise ) as follows: (cid:96) ( d obs ; m ) = 12 ( F ( m ) − d obs ) T Γ − ( F ( m ) − d obs ) , (12)where F is the discretized parameter-to-observable map. The negative log of the like-lihood, (cid:96) , is also known as the “data misfit term” in the deterministic inversion setting.A typical choice of Γ noise is a diagonal matrix with entries σ i , the noise-level varianceof the i th observation d obs i . This choice of noise is valid for cases in which the noise pol-luting one measurement is uncorrelated with the noise in other measurements. As a prior, we choose a Mat´ern class that can represent several commonly used co-variance models used in geostatistics (Goovaerts, 1997; Deutsch & Journel, 1998; Wack-ernagel, 2010). Our representation of the Mat´ern prior is based on the link between theMat´ern class of Gaussian fields and solutions of elliptic stochastic partial differential equa-tions (SPDE) (Lindgren et al., 2011). We define the parameter m ( x ) to be a Mat´ern classGaussian random field in the domain Ω ∈ IR d , d = 1, 2, or 3. Lindgren et al. (2011)show that its covariance operator is the square of the solution operator of the SPDE( δ − γ ∆) α m ( x ) = s ( x ) in Ω . (13)For d = 3, α should satisfy α > /
2, we choose α = 2. The right-hand side of theSPDE (13), s ( x ), is a white-noise Gaussian random field with unit marginal variance. –6– he parameters γ , δ , and α are chosen to achieve the desired marginal variance and cor-relation length of the features of the parameter field.The parameter m is discretized over the same triangularization defined for the lin-ear poroelasticity PDE problem. We use first-order Lagrange polynomials to approxi-mate m . For our choice α = 2, the Mat´ern class precision operator, R , is constructedas follows: R = AM − m A . The operator A is the finite element discretization of thedifferential operator A = ( δ − γ ∆) and M m is the mass matrix in the finite elementspace where m is approximated. M / m can be efficiently approximated (Villa et al., 2019),in which case the cost of sampling from the prior is mainly the cost of solving the ellip-tic PDE (13) for a realization of the white noise right hand side, which can be carriedout with multigrid solvers which have linear complexity. The covariance of the prior, Γ prior ,is given by R − , the inverse of the precision operator.For an isotropic Gaussian field of Mat´ern class, the correlation length (or range)is ρ = (cid:113)(cid:0) γ ( α − d ) (cid:1) /δ ; in particular, m ( x ) and m ( x ), for any x , x ∈ Ω of dis-tance ρ from each other, have approximately a correlation near 0 . d = 3) and we choosethe elliptic operator power α/ ρ reduces to 2 (cid:112) γ/δ . Choosing γδ = 10 ,for example, dictates a correlation near 0.1 for points that are about 2 ,
000 length unitsapart.
Interferometric SAR (InSAR) computes the phase difference between two syntheticaperture radar (SAR) images. The resulting interferogram can be used to infer a 2D mapof surface deformation between two SAR acquisition times along the radar line-of-sight(LOS) direction. The LOS deformation, u LOS , can be described by the formula: u LOS = α u + α u + α u , (14)where u , u and u are the east, north and vertical displacements, accordingly. The vec-tor α = [ α , α , α ] is the radar LOS direction unit vector that can be calculated basedon the known radar imaging geometry.We define the pointwise InSAR data misfit function (i.e., the negative log likeli-hood), (cid:96) InSAR , as follows: (cid:96)
InSAR ( m ) = n interfero (cid:88) i =1 σ i InSAR ) n i pixels (cid:88) j =1 (cid:16) u i LOS ,j − u i InSAR ,j (cid:17) , (15)where n interfero is the number of interferograms used in the inversion and n i pixels is thenumber of pixels in the i th interferogram. In this data misfit expression, we assume thecovariance operator of the noise is diagonal, where σ i InSAR is the noise level of the i th in-terferogram. This assumption is valid when InSAR measurement noise mostly stems fromphase decorrelation as in the Nevada test case, presented in section 3. In cases where In-SAR measurement noise contains a non-negligible spatially coherent component (e.g. noisefrom atmospheric delays), the covariance operator can be constructed to dictate that noisecorrelation. The value u i InSAR ,j is the LOS deformation measured from the j th pixel inthe i th interferogram, and u i LOS ,j is the simulated LOS deformation of the j th pixel inthe i th interferogram as computed from the forward poroelasticity model. The discretized posterior (10) in general is a distribution in high dimensional pa-rameter space—thousands or millions. It is difficult to explore such distributions due to –7– he high dimensionality and the need to solve the 3D poroelasticity system (1)–(3) toevaluate the posterior at each point in parameter space. Here in Part I, we seek to de-termine the MAP point, which as we will see, amounts to a deterministic inverse prob-lem with specially-chosen weighting matrices. We employ adjoint-based derivatives, theGauss–Newton method, and early termination of conjugate gradient (CG) iterations ateach Gauss–Newton step to ensure the efficiency and scalability to high parameter di-mensions for solving this deterministic inverse problem.The MAP point, m MAP , is the point that maximizes the posterior distribution (10).We can define the parameter-to-observable map, F ( m ), as follows: F ( m ) := ¯B ¯X , wherethe state variables ¯X depend on the permeability parameter vector m through the for-ward problem ¯S ¯X = ¯F (7), and ¯B is the linear observation operator that extracts theobservations from ¯X , i.e. ¯B ¯X is the surface displacements at times and locations of theobserved data. Thus we can rewrite the the posterior distribution (10) as: π post ( m | d obs ) ∝ exp (cid:18) −
12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs ) −
12 ( m − ¯ m ) T Γ − ( m − ¯ m ) (cid:19) . (16)Maximizing the posterior distribution with respect to the parameter m is equivalent tominimizing the negative log posterior. Therefore, our goal is to minimize the objectivefunction:min m ∈ IR nn J ( m ) = 12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs ) + 12 ( m − ¯ m ) T R ( m − ¯ m ) , (17)where ¯X depends on m through the solution of the forward problem ¯S ( m ) ¯X = ¯F . Itcan be seen that Γ − = R acts as a regularization operator in the deterministic prob-lem (17)(Tarantola, 2005). To derive the gradient using the adjoint method, we constructa Lagrangian by adding to the objective function J ( m ) the inner product of the left-hand side of the state equation ¯S ( m ) ¯X − ¯F = 0 with a Lagrange multiplier (also calledthe adjoint variable) ¯Y . In the Lagrangian formulation, ¯X is considered an independentvariable, and the dependence of ¯X on the parameters m is enforced through the Lagrangemultiplier term.At a minimum of (17), the partial derivatives of the Lagrangian with respect to thestate variable, ∂ L /∂ ¯X , the parameter, ∂ L /∂ m , and the adjoint variable, ∂ L /∂ ¯Y , van-ish. We set ∂ L /∂ ¯X and ∂ L /∂ ¯Y to zero, which gives the state equation and the adjointequation respectively. We satisfy these two conditions by solving the discretized stateand adjoint equations exactly for any value of m . Thus we seek the parameter m thatensures that the gradient ∂ L /∂ m vanishes. We solve the equation ∂ L /∂ m = 0 usingNewton method. The details of forming the Lagrangian and constructing the Newtoniteration are left to Appendix B.The Newton system we need to solve at each iteration l is: H l δ m l = − G l , (18)where G l = G ( m l ) = ( ∂ L /∂ ¯X ) l is the gradient and δ m l is the Newton direction. Thesuperscript l we use for the matrices, and the forward and the adjoint solutions indicatesthat these matrices and vectors are evaluated at the l th iteration MAP point approxi-mation, m l . The Hessian H l = H ( m l ) ∈ IR n n × n n is given by H = R + R mm + ¯C T ¯S − T ( W ¯X ¯X ¯S − ¯C − W ¯Xm ) − W m ¯X ¯S − ¯C . (19)The operator ¯C is the partial derivative (sensitivity) of the residual R = ¯S ¯X − ¯F withrespect to the parameter m ¯C = ∂∂ m ( ¯S ¯X − ¯F ) = ∂∂ m ( ¯S ¯X ) , (20) –8– nd the other matrices 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 ) . The direction δ m l is used to update the approximation of the MAP point m l asfollows: m l +1 = m l + θδ m l , (21)where θ is a step length parameter chosen to satisfy the Armijo rule, to guarantee a suf-ficient reduction in the objective function (Nocedal & Wright, 2006). Backtracking linesearch is used to reduce the step size until the objective function is sufficiently reduced.The conjugate gradient method with early termination is used to solve the linear sys-tem (18) at each Newton iteration. This CG iteration is also terminated when a direc-tion of negative curvature is detected, which ensures that the direction δ m l is a descentdirection (Nocedal & Wright, 2006).The full Hessian (19) can be approximated by the Gauss–Newton Hessian H GN ,in which all matrices depending on the adjoint variable ¯Y , i.e. W ¯Xm , W m ¯X , and R mm ,are neglected. The reasoning behind this is that, when the data misfit ¯B ¯X − d obs is zero(i.e. when the model fits the data), the right hand side of the adjoint equation (B3) van-ishes, and since the equation is linear in the adjoint variable ¯Y , it vanishes. Thus whenthe data misfit is small, the Gauss–Newton approximation can be expected to be accu-rate. The resulting approximation is guaranteed to be positive definite for appropriatelychosen regularization, is computationally less expensive, and is expected to converge su-per linearly near the optimal point, m MAP , in the zero-noise case. However, if the noiselevel is high, we can switch to the full Newton Hessian after taking some initial itera-tions with H l GN . The GN Hessian is given by H l GN = R + ¯C l T ¯S l − T ( W ¯X ¯X ¯S l − ¯C l ) . (22)Forming the GN Hessian (and similarly the full Hessian) explicitly each Newtoniteration l (18) is computationally intractable. It requires applying the operator ¯S l − toeach column of ¯C l to form the matrices product ¯S l − ¯C l , each Newton iteration l . Ap-plying the operator ¯S l − to a vector amounts to an incremental forward poroelasticitysolve, equation (B9). Therefore, forming the GN Hessian explicitly costs n n poroelas-ticity PDEs solves (note that no additional solves are required for forming the product ¯C l T ¯S l − T since it is the transpose of the product ¯S l − ¯C l ). When using CG, there is noneed to form the GN Hessian in each Newton iteration l explicitly. What is required, ineach CG iteration i , is computing the matrix vector product R − H l GN d i , where d i is the i th conjugate gradient direction and R − is the prior covariance matrix which is com-monly a favorable choice as a preconditioner for the Hessian. The main cost of comput-ing this product is the cost of applying the operators ¯S l − and ¯S l − T each to a vector wichis the cost of two poroelasticity PDEs solves (one incremental forward (B9) and one in-cremental adjoint (B10)) each CG iteration. We emphasize that it is often the case thatthe number of CG iterations n CG required in each (Gauss–)Newton iteration l is muchsmaller than, and independent of, the number of parameters n n ( n CG (cid:28) n n ) in manyPDE-based optimization problems including those governed by poroelasticity, as has beenproven for some particular problems and observed numerically for a number of others(Bashir et al., 2008; Flath et al., 2011; Bui-Thanh & Ghattas, 2012b, 2013, 2012a; Bui-Thanh et al., 2012, 2013; P. Chen et al., 2017; Alexanderian et al., 2016, 2017, 2014; Cres-tel et al., 2017; Petra et al., 2014; Isaac et al., 2015; Martin et al., 2012; Bui-Thanh &Ghattas, 2015; Hesse & Stadler, 2014). As a consequence of ill-posedness, a typically smalland parameter dimension-independent number of CG iterations is required, correspond-ing to the dominant eigenvalues of the prior preconditioned (Gauss–)Newton Hessian that –9– re distinct from those that cluster around one. For sufficiently large CG tolerance, theeigenvector components of the error corresponding to the remaining eigenvalues clusteredaround one can be eliminated in O (1) CG iteration (Nocedal & Wright, 2006). Further-more, in early Newton iterations when the MAP point approximation m l is away fromthe true solution, solving the linear system arising in the Newton iteration l requires evenfewer CG iterations than the number of the dominant eigenvalues ∼ n CG due to impos-ing Eisenstat-Walker tolerance that prevents oversolving the linear system (Eisenstat &Walker, 1996; Nocedal & Wright, 2006; Villa et al., 2019). In this section, we first review an aquifer pumping test conducted in 2003 by Burbeyet al. (2006) near Mesquite, Nevada. We then introduce the InSAR and GPS data usedto measure the aquifer deformation induced by this test. We finally describe the aquifermodel, and how we set up the Bayesian inversion framework to infer the lateral perme-ability variations of this aquifer.
The Nevada study site is located southeast of Mesquite, Nevada, in the Mesquitebasin of the lower Virgin River valley (Figure 1 a and 1 b ). The aquifer is ∼
400 meterthick, which is confined by a brittle unsaturated layer at the top and a clay layer at thebottom (Figure 1 c ). Hydraulic head records suggest that the regional GW recharge isdue mainly to winter precipitation in the surrounding mountains, and hence is negligi-ble during the duration of the well test (Burbey et al., 2006; Burbey, 2008).A controlled aquifer test was performed at a newly installed municipal well WX31between May 7, 2003 and July 9, 2003 (Burbey, 2006). The well was screened over a 237 . ,
028 cubic meters of water daily. The pumping continued until ∼ October, 2003. Notethat the aquifer was isolated from poromechanical stresses exerted by pumping activ-ities in adjacent wells. Over the test period and within a radius of ∼ ∼ We generate an interferogram spanning May 4, 2003 and October 26, 2003 to mea-sure the LOS deformation caused by the 2003 aquifer test. The LOS unit look vectoris [0 . , − . , . d ). This asymmetrical pattern suggests thepossible existence of a fault as reported in Burbey (2008). Our interferogram shows sim-ilar deformation patterns as the one that spans May 4, 2003 and September 5, 2004 (Burbey,2008). This suggests that little deformation occurred after the end of the well test be- –10– ° . ′ N ° . ′ N ° . ′ N ° . ′ N Well WX31 GPS station (a)
Well (b)
Mesquite
Confining layer (85 m) (c)
Aquifer (400 m) (e) (d)
Lower clay (400 m)
Figure 1: (a) The study site near Mesquite, NV. The red polygon shows the lateralboundaries of our numerical simulation domain. (b) Zoom-in view of the area outlinedin cyan in Figure 1 a . The red star marks the location of well WX31, where Burbey et al.(2006) performed the aquifer test in 2003. Surface deformation was monitored daily usinga network of 10 high-precision GPS stations (white dots). (c) A conceptual illustrationof the aquifer model. The screened segment of the well is marked in blue. (d) InSAR-observed LOS deformation between May 4, 2003 and October 26, 2003 over the studyarea. The the gray polygon marks the extent of the data used in the inversion. (e) A mapview of the observed cumulative lateral displacements at 6 GPS stations that are closer tothe well from day 1 to day 22 of the pumping test. Each black dot shows the daily cumu-lative lateral displacement. The gray shaded circles enclosing the GPS data are the 99%confidence ellipses of the measured data. In both Figure 1 d and 1 e , x is the east directionand y is the north direction relative to the well location (x,y)=(0,0). –11– ween October 26, 2003 and September 5, 2004, confirming that the regional GW rechargeis minimal.We select a high phase-correlation ( > .
3) block bounded by the gray outline inFigure 1 d as the input for the inversion. In this 8-by-7 km region, there are no visibleorbital errors, ionospheric artifacts, or phase unwrapping errors that impact the defor-mation estimates. DEM errors and tropospheric errors are not substantial, because thestudy site is relatively flat with a dry desert climate. The dominant error source hereis due to the change of surface scattering properties between the two SAR acquisitions.This phase noise (known as phase decorrelation) is not correlated in time or space (Zebker& Villasenor, 1992). We estimate the phase decorrelation noise level ( ∼ σ InSAR in equation (15).We use the lateral displacements as recorded at 6 GPS stations during the first 22days of the 2003 pumping test (Figure 1 e ) as an independent validation of the InSAR-based permeability estimates. These GPS stations are closer to the well, and little de-formation was recorded at the other 4 stations that are further away. The joint inver-sion of GPS and InSAR data sets as well as their relative merits will be discussed in afuture article. We model the aquifer site as an 885 m deep layered poroelastic medium (as shownin Figure 1 c ) over the region bounded by the red polygon in Figure 1 a . We assume noGW flux at all boundaries, zero displacement at the lateral boundaries, zero normal dis-placement at the bottom boundary, and a traction-free top surface. We set the initialdeviation from the hydrostatic pressure before the start of the pumping, p , to zero. Wemodel the groundwater extraction as a volumetric sink term, f p in equation (1), locatedin the segment of the well coinciding with the lower half of the aquifer (the blue segmentof the well in Figure 1 c ). Following Burbey (2008), the aquifer parameters are listed inTable 1.We set up the Mat´ern class prior Gaussian field defined in equation (13), which yieldspermeability samples with the following properties: (1) on average, at any given pointin space, the permeability values vary from the prior mean listed in Table 1 b by a ∼ . ρ = 2 km in the lateral direction and an or-der of magnitude larger correlation length in the vertical direction. The large verticalcorrelation length suppresses vertical variations that are not informed by the surface de-formation data. Therefore, the estimated permeability should be interpreted as the ver-tical average of the permeability field. This average is a useful up-scaled representationof the aquifer permeability, because the GW flow is predominantly horizontal and alongthe unresolved fine-scale horizontal layering in the aquifer (Cherry & Freeze, 1979). Insection 5.2, we show that the inverse solution is robust irrespective of a wide range ofSD and ρ in the horizontal direction.We generate two unstructured tetrahedral meshes of differing resolution for our com-putational domain using Gmsh (Geuzaine & Remacle, 2009). The number of nodes in thesetwo meshes, and hence the number of DOFs approximating the log permeability m , are4 ,
081 and 16 ,
896 (Figure 2 a ). The overall number of DOFs at each time step of the dis-crete system state variables ( p , u and q ) are 72 ,
980 and 320 , , –12– able 1: 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.2axis as the well. The sink term is distributed over that volume such that it integratesto the pumping rate. We discretize the temporal evolution into 154 variable-length timesteps over a period of 175 days. At the start of the simulation, the time step is set to1.2 hours to accurately capture the rapid changes after the onset of pumping. The timestep length increases gradually to up to 5-day time interval toward the end of the sim-ulation.We implement the forward model, adjoint system, and gradient and Hessian eval-uations using the Python -based
FEniCS library (Logg et al., 2012) for FEM discretiza-tion in space. We use the hIPPYlib library for state-of-the-art Bayesian and determin-istic PDE-constrained inversion algorithms and prior construction (Villa et al., 2018).Using this implementation we infer the log-permeability MAP point, i.e. the point in theparameter space that maximizes the posterior distribution, from the InSAR data. Basedon this inferred permeability field, the forward model can simulate a three dimensionalevolution of the pore pressure in response to the pumping test (e.g. Figure 2 b ) and theexpected deformation in east, north, and vertical directions at each time interval. Wethen compare this model-based lateral deformation with the GPS data set for validationof the inverse solution. In section 4.1, we present results of applying the Bayesian inversion framework toinfer the permeability field of the Nevada test case described in section 3. Section 4.2demonstrates the superiority of Newton-type methods over steepest descent for these poroe-lastic inverse problems. –13– igure 2: (a) The computational mesh of 16 ,
896 nodes. The mesh is generated using
Gmsh and is refined around the well. (b) Simulated three-dimensional aquifer pore pres-sure profile at the time of the second InSAR data acquisition (5 . − ,
000 to − , − ,
000 to − , − ,
000 to 0 (Pa), respectively. These iso-volumes are cropped at z = . We solve for the log-permeability MAP point, m MAP , using the Envisat interfer-ogram that spans May 4, 2003 and October 26, 2003 (Figure 1 d ). The two-dimensionalhorizontal slice of the reconstructed permeability, κ = e m MAP , at mid-aquifer depth (Fig-ure 3 a ) reveals distinct features in the permeability field: A high permeability channelextending from south/southwest of the well (marked by red star) to the west; and a low-permeability flow barrier northwest of the well.Using this InSAR-inferred permeability field and the forward model, we simulatethe expected LOS deformation that occurred between May 2003 to October 2003. Thisreconstructed LOS displacement accurately captures the subsidence bowl due to the pump-ing test. In most of the region, the discrepancy between the InSAR data and the sim-ulated LOS deformation (Figure 3 c ) is roughly random noise of magnitude within theestimated noise level. We observe a relatively large discrepancy in the northwest, likelydue to the existence of a fault (Burbey, 2008) that is not accounted for in our aquifermodel. In section 5.3, we discuss the impact of the model error in our inversion solution.Note that the simulated LOS deformation is computed from the simulated east,north, and vertical deformations (Figure 3 d –3 f ) following equation (14). To validate ourinversion results, we further compare the simulated lateral deformation with the lateraldisplacements as recorded at 6 GPS stations (Figure 3 g ). We confirm that the recon-structed lateral displacements capture the GPS-observed nontrivial southeast-trendingsurface deformation during the first 22 days of the 2003 pumping test. It is worth not-ing that the northward component of these GPS-observed lateral displacements is well-approximated by our model despite the fact that the northward contribution to the LOSdeformation is negligible ( α = − . –14– .2 Convergence of the (Gauss–)Newton Method The (Gauss–)Newton conjugate gradient method, described in section 2.2.4, is knownto converge rapidly when applied to optimization problems for which the data misfit Hes-sian effective rank is relatively small and independent of the discretization dimension (Heinkenschloss,1993). As has been mentioned in section 2.2.4, this has been demonstrated either an-alytically or numerically for a wide spectrum of geophysical inverse problems. Althoughthis property of the Hessian has not been shown analytically for poroelastic inverse prob-lems, our numerical results confirm previous work by Hesse and Stadler (2014) that demon-strated the computational advantage of Newton-type methods, suggesting rapid decayof the eigenvalues of the data misfit Hessian for this problem.To verify the dimension-independent convergence of (Gauss–)Newton–CG itera-tion, in the form of equation (18), we solve for the permeability MAP point for both ourcoarse and fine meshes and compare their convergence in the blue and green curves inFigure 4 a ). In both cases, the (Gauss–)Newton–CG method converges in a small num-ber of iterations (less than 60). Even though the fine mesh has four times the numberof parameter DOFs as the coarse mesh, the number of iterations and PDE solves requiredfor both meshes are comparable (Table 2). This demonstrates that the algorithm is di-mension independent and hence scalable to larger problems. This is essential if the InSAR-based Bayesian poroelastic aquifer characterization is to be applied to GW managementat the basin scale (J. Chen et al., 2016, 2017). Additionally, the inferred permeabilityfields in both the coarse mesh (Figure 3 a ) and the fine mesh (Figure 4 b ) are consistent.Note that the method was set to use the Gauss–Newton Hessian in the first 50 iterations,then switch to using the full Newton Hessian in the remaining iterations.We also demonstrate the superiority of the (Gauss–)Newton–CG method’s conver-gence over the classical steepest descent method on the coarse mesh (red line in Figure 4 a ).As shown in Table 2, the steepest descent method terminated after 18 days of runtime(due to hitting the iteration limit of 5 , . × (butthe true speed up would be much larger had we allowed steepest descent to run to con-vergence). Comparing the number of PDE solves in each case, we see that steepest de-scent method required 30 times the number of PDE solves required by the (Gauss–)Newton–CG method and yet did not converge (Table 2).We point out that part of the computational savings obtained when using the (Gauss–)Newton–CG method are due to the cost of a PDE solve for the steepest descent methodbeing greater than that for the (Gauss–)Newton–CG method in our implementation. Onaverage, a PDE solve requires 20 seconds in the former case, while requiring only 7 sec-onds for the latter case on the coarse mesh. This is a result of using a direct solver forthe linear poroelasticity equations. When using the GN method, the direct solver’s LUfactors can be computed once at the l th GN iteration, stored, and reused in all the CGiterations required at that GN iteration. This is because the same incremental forwardand adjoint PDEs are solved at each CG iteration (equations B9 and B10). On the otherhand, for the steepest descent method, the parameter m l changes in each iteration andthus a new factorization of the poroelasticity operator is required at each steepest de-scent iteration. –15– igure 3: Results of characterizing Nevada aquifer permeability using InSAR data. In allthe panels, the white circles mark the GPS station locations and the red star marks thewell location. (a) Two-dimensional horizontal slice of the inferred permeability, e m MAP , atmid-aquifer depth (in decimal logarithm). (b) Reconstructed LOS displacement from theforward model using the permeability profile in Figure 3 a . (c) The discrepancy betweenInSAR data and the reconstructed LOS displacement: u disc = u InSAR − u LOS . (d)–(f)Displacement components ( u , u , and u ) of the LOS deformation in Figure 3 b . (g) Mapview of the reconstructed (blue) versus observed (black) cumulative lateral displacementsat the GPS station locations from day 1 to day 22 of the pumping. At each GPS station,each black dot shows the daily cumulative lateral displacement. In all the panels, x is theeast direction and y is the north direction relative to the well location (x,y)=(0,0) –16– igure 4: (a) Convergence of the (Gauss–)Newton method for finding the InSAR-data-based MAP point using the coarse mesh, 4 ,
081 nodes, (blue) and the fine mesh, 16 , norm of the gradient byfour orders of magnitude. (b) A two-dimensional horizontal slice of the InSAR-data-basedMAP point, e m MAP , at mid-aquifer depth (in decimal logarithm) using the 16 ,
896 nodemesh. –17– able 2: Convergence results for computing the MAP point for three cases. The first tworows report convergence results of the (Gauss–)Newton method using the coarse (4 , ,
896 nodes) mesh, respectively. The third row reports con-vergence results of steepest descent method (SDM) using the coarse mesh. The gradientreduction factor is defined as initial gradient norm / final gradient norm.Method, Gradient PDE solves a : Num. of global Num. of total Totalmesh reduction forward (fwd), iterations CG it. timefactor adjoint (adj)GN, ≥ fwd: 107 + 1 , ,
078 4.5hcoarse adj: 51 + 1 , ≥ fwd: 103 + 1 , ,
164 22.3hfine adj: 56 + 1 , ≈ . , ,
000 - 18dcoarse adj: 5 , a The number of forward and adjoint PDE solves are given as: number of forward(or adjoint) poroelasticity solves + number of incremental forward (or adjoint)poroelasticity solves.
In this section, we discuss how the essentially arbitrary pixel spacing of the InSARdata can be treated systematically in the Bayesian setting, without the introduction ofartificial weights. We then show that the inferred MAP point of the permeability fieldis not appreciably affected by the choice of prior. We further show that the discrepancybetween the model and the data beyond the noise identifies a region where our modelexhibits structural error. We conclude this section with a brief overview of Part II of thisstudy, which focuses on MCMC sampling of the Bayesian posterior, building on tools de-veloped in the present article.
Here we study the sensitivity of the inverse solution to the essentially arbitrary pixelspacing of the InSAR data. We downsample the original InSAR data multiple times (Fig-ure 5). Each level of downsampling is performed by averaging adjacent two-by-two pix-els into a single pixel. This is known as multi-looking, which reduces the InSAR phasedecorrelation noise by a factor of two. We incorporate this noise adjustment in (15). Fig-ure 6 shows the resulting inferred permeability MAP point at mid-aquifer depth usingeach of the six multi-looking InSAR data. When the number of InSAR pixels is suffi-cient to resolve the main deformation features (Level 0-3), the inferred permeability so-lutions are similar. Hence we conclude that the InSAR data redundancy can be addressedin a fully Bayesian setting, i.e., without introducing artificial weights, by treating the datanoise appropriately. However, if InSAR data are too coarse to resolve the main defor-mation features (Level 4-5), a deterioration in the validity of the inferred permeabilityfield is expected. –18– igure 5: Downsampled InSAR data sets labeled by the downsampling level. At eachdownsampling level, adjacent two-by-two data points are averaged into a single datapoint. Level zero is the original data set (31 ,
003 data points). The number of data pointsare 7 , , The limited data available for aquifer characterization has placed emphasis on meth-ods to inform the characterization with prior geostatistical knowledge (Goovaerts, 1997;Deutsch & Journel, 1998; Wackernagel, 2010). A large number of methods that allowimposition of additional knowledge from small-scale horizontal layering to the connec-tivity of large-scale fluvial channels (Remy et al., 2009; Mariethoz & Caers, 2014) havebeen developed. Unfortunately, most of these methods are not compatible with the scal-able framework used here that requires a sparse prior precision operator. However, weshow below that InSAR data are sufficiently informative that dense geostatistical pri-ors are not required.Here we limit ourselves to the class of Mat´ern covariances discussed in section 2.2.2.These priors have the typical parameters of a semi-variogram (Deutsch & Journel, 1998),the range, ρ ∝ (cid:112) δ/γ , and the sill (Figure 7 a ). The sill is the square of the pointwisemarginal standard deviation SD ∝ (cid:112) / ( δγ ). Due to the boundary effects, the exactsill cannot be determined a priori and is hence difficult to control exactly in Mat´ern pri-ors. Due to the continuous nature of the prior there is no nugget. In addition, the Mat´ern –19– igure 6: Two-dimensional horizontal slices of the permeability MAP point, e m MAP ,at mid-aquifer depth (in decimal logarithm) inferred from each of the six downsampledInSAR data sets in Figure 5.prior also allows us to specify a mean, which is taken from Burbey (2008). Two-dimensionalhorizontal slices of prior samples with SD ≈ . ρ = 2 to ρ =4, are shown in Figure 7 b and 7 c , respectively. We also note that in infinite dimensionalinverse problems the prior must act as a regularization. Table 3 notes priors that do notregularize the inverse problem sufficiently so that MAP point estimate does not convergeto the required tolerance.The high spatial resolution of the InSAR data provides good constraints on the lat-eral permeability variation in the aquifer. This can be seen in Figure 8, which shows theMAP points for each of the five prior choices in Table 3. Despite large ranges in both ρ and SD, the main features of the MAP point remain almost unchanged. This includesboth the high permeability channel south of the well and the low-permeability region north-west of the well. Of course, the amplitudes increase with SD and small features are lostas ρ increases, but main pattern is independent of the prior. This demonstrates that theMAP point is dominantly informed by the data and that the effect of the prior on theinferred permeability variation is small.This highlights the step change in the data availability for the characterization oflateral aquifer heterogeneity that is provided by InSAR. Given upcoming missions with –20– able 3: Parameters for five different prior choices. The SD values are reported as threevalues: minimum (third column), maximum (fourth column), and average over domain(fifth column). The values of γ and δ are 3 .
33 and 3 . −
6, respectively. The second tolast column reports whether inverting for the MAP point has converged (yes) or stoppedwhen reaching the maximum number of backtracking steps (no). The convergence crite-rion we impose is the reduction of the L norm of the gradient by 4 orders of magnitude. γ δ SD SD SD ρ (cid:112) / ( δγ ) Converged Gradientmin. max. avg. (km) reduction γ δ ≥ γ δ ≈ . γ δ ≥ γ δ ≈ . γ δ ≥ improved accuracy and increased frequency, InSAR-derived surface deformation mea-surements will provide a powerful and low-cost means of aquifer characterization. In section 4.1, we showed that the discrepancy between the model and the data isroughly a random noise signal of magnitude within the estimated noise level, except inthe northwest corner of the domain (Figure 3 c ). Burbey (2008) have postulated the ex-istence of a fault in this area, which is not accounted for in our aquifer model, presentedin section 3.3. Here we investigate whether the relatively large discrepancy in the north-west corner is due to this structural error in the model.Introducing the actual fault and the associated mechanics is beyond the scope ofthis article. Instead, we explore the possibility of model error by studying a syntheticinversion. For this model-error-free scenario we regard the permeability field inferred fromthe real InSAR data (Figure 3 a ) as the “true” permeability. Using this truth permeabil-ity field, a synthesized InSAR LOS displacement map is generated from the model-predictedLOS observations over the same period and area as for the actual InSAR observations.These synthesized observations are then polluted by a normally-distributed random noiseof the same magnitude assumed for the real InSAR data, 3 . a ).We then use the synthesized LOS data to solve the inverse problem of estimatingthe permeability MAP point. The overall pattern of the resulting permeability MAP point(Figure 9 b ) is very similar to those of the “true” medium (Figure 3 a ). There is, how-ever, a notable decrease in the magnitude of the permeability variation from the priormean, due to the regularizing action of the prior.Using the permeability field inferred from the synthetic data (Figure 9 b ), we solvethe forward problem and reconstruct the LOS observations, u LOS . The error in these re-constructed LOS observations (Figure 9 c ) is within the noise everywhere in the observedarea. This demonstrates that we should be able to reconstruct the LOS observations towithin data noise in the absence of model error. This suggests that the inability of theporoelatic model to reconstruct the LOS observations in the northwest (Figure 3c) is likelydue to a number of possible simplifications in the model, including the presence of a fault. –21– igure 7: (a) Averaged semivariance over 100 samples from the three prior choices with ρ = 1, 2, and 4 and an average-over-domain SD ≈ . ρ and similar SD.Again this demonstrates the high information content of the InSAR data. Not onlydo they provide detailed constraints on the lateral permeability variation, they also clearlyhighlight regions where the model likely has structural errors. This allows the targetedimprovement of the aquifer model. Our focus in part I of this two-part series of articles is on estimating the perme-ability maximum A posteriori point. In the preceding discussion, we have paid carefulattention to matters arising in the Bayesian formulation of this inverse problem, includ-ing treatment of multilooking data sets, the construction of the prior and its influenceon the MAP point, and analysis of model error. Our ultimate goal is to explore the pos-terior distribution and estimate the expected value along with the associated uncertain-ties of various quantities of interest. Estimating the MAP point is a critical first step be-cause of its intrinsic importance as the most likely permeability realization. Furthermore,the MAP point can be used to create proposals for MCMC sampling algorithms required –22– igure 8: Two-dimensional horizontal slices of the permeability MAP point, e m MAP , atmid-aquifer depth (in decimal logarithm) for the five chosen priors in Table 3. To makethe visual comparison easier, the MAP point for which (avg. SD, ρ ) = (1.4, 2 km) isshown twice: in the middle of the top and bottom rows. (a) The average-over-domainSD decreases from left to right: 2 .
8, 1 .
4, to 0 .
7, while ρ = 2 is constant. (b) ρ increasesfrom left to right: 1, 2, to 4 kilometers, while the average-over-domain SD ≈ . –23– igure 9: (a) Synthesized InSAR data generated from the forward model with the per-meability field in Figure 3 a . (b) Two-dimensional horizontal slice at mid-aquifer depthof the permeability, e m MAP , inferred using the data in Figure 9 a (in decimal logarithm).(c) The discrepancy between the synthesized InSAR data and the reconstructed LOS dis-placement: u disc = u InSAR − u LOS , where u LOS is reconstructed from the model with thepermeability field in Figure 9 b . We have shown that InSAR-based surface deformation measurements provide de-tailed information about lateral variations in aquifer permeability. In the Nevada testcase, the surface deformation due to a single well allowed us to identify both high per-meability pathways and flow barriers. Due to the high spatial resolution of the InSARdata, the main features of the inferred MAP permeability field are essentially indepen-dent of prior assumptions. The inferred permeability field also allows us to reconstructthe lateral deformations recorded at six GPS stations near the well. This provides val-idation, because the GPS data were not used in the inference of the permeability field.The above conclusions are derived from a single interferogram that captures thetotal subsidence due to the well test. Given the increasing availability of InSAR data fromSentinel-1 and upcoming NISAR missions, future aquifer characterization will be ableto utilize higher quality and more frequent surface deformation measurements. We there-fore believe that geodetic surface deformation measurements will dramatically increasethe information for aquifer characterization. Scalable and robust inversion frameworksthat integrate InSAR data into poromechanical aquifer models are likely to become crit-ical tools for regional groundwater management.The Bayesian inference framework was applied to a three-dimensional transient multi-physics problem with real data that requires the inference of up to 16 ,
896 parameters.This requires dimension invariant methods and was achieved by exploiting the compactnature of the parameter-to-observable map, adjoint-based derivatives, and Gauss–Newton–CG method, the combination of which capture the intrinsic low dimensionality of thisinverse problem and permit rapid, dimension-independent convergence. The power ofthese methods is available through the
Python -based Bayesian inverse problems library hIPPYlib . The scalability of this framework provides the basis for Bayesian uncertainty –24– uantification, based on Markov chain Monte Carlo sampling, which we will explore inPart II of this two-part series of articles.
Appendix A Variational Formulation and Discretization Details of the3-Field Formulation of the Biot System
For the variational form of the system (1)–(3) to be well-defined, we assume thestate variables p ( . ; t ), u ( . ; t ) and q ( . ; t ) belong to the following infinite dimensional func-tion spaces (note that these are regularity assumptions in space, we additionally requiresufficient regularity in time): P = (cid:8) p : Ω → IR | p ∈ L (Ω) (cid:9) (A1) U = (cid:8) u : Ω → IR d | u ∈ (cid:0) H (Ω) (cid:1) d and u ( x ; t ) = u d on ∂ Ω du (cid:9) (A2) Q = (cid:8) q : Ω → IR d | q ∈ H (div; Ω) and q ( x ; t ) · n = g on ∂ Ω np (cid:9) , (A3)respectively, where d = 1, 2, or 3 is the dimension of the physical domain. The space H (div; Ω) is defined as H (div; Ω) = (cid:110) f : f ∈ (cid:0) L (Ω) (cid:1) d and ∇ · f ∈ L (Ω) (cid:111) .To obtain the weak form of the system (1)–(3), we multiply the equations (1), (2)and (3) by test functions r ∈ P , v ∈ U and w ∈ Q , respectively, where the func-tion spaces P , U , and Q are defined as follows: P = (cid:8) p : Ω → IR | p ∈ L (Ω) (cid:9) (A4) U = (cid:8) u : Ω → IR d | u ∈ (cid:0) H (Ω) (cid:1) d and u ( x ; t ) = on ∂ Ω du (cid:9) (A5) Q = (cid:8) q : Ω → IR d | q ∈ H (div; Ω) and q ( x ; t ) · n = 0 on ∂ Ω np (cid:9) . (A6)We then integrate the three equations in space over the domain Ω. We discretize in timeusing backward Euler method. The variational problem is therefore the problem of find-ing p k ∈ P , u k ∈ U and q k ∈ Q that satisfy the following equations: (cid:90) Ω (cid:0) ( S (cid:15) p k + α ∇ · u k ) r + ∆ t k ∇ · q r (cid:1) d x = (cid:90) Ω (∆ t k f p + S (cid:15) p k − + α ∇ · u k − ) r d x (A7) − (cid:90) Ω ( σσσ ( u k ) − αp k I ) : ∇ v d x = − (cid:90) Ω f u · v d x − (cid:90) ∂ Ω n u g · v d s (A8) − (cid:90) Ω ∆ t k ( e m µ ) − q k · w d x + (cid:90) Ω ∆ t k p k ∇ · w d x = (cid:90) ∂ Ω dp ∆ t k p d w · n d s, (A9)for all the choices of test functions r ∈ P , v ∈ U and w ∈ Q and for time steps k = 1 , , , .., N t . The functions p and u are the initial values for the pressure and thedisplacement, respectively. Requiring that f p ( ., t ) ∈ L (Ω), p d ( ., t ) ∈ L ( ∂ Ω dp ), f u ( ., t ) ∈ ( H − (Ω)) d , and g ( ., t ) ∈ ( H − ( ∂ Ω n u )) d is necessary for the well-posedness of the weakproblem. We follow (Ferronato et al., 2010) in our choice of the finite element functionspaces in which the pressure p is approximated by piecewise constant function ˆ p ( x, t ) inˆ P ⊂ P and the velocity is approximated in the lowest order Raviart-Thomas space bythe function ˆ q ( x, t ) in ˆ Q ⊂ Q . The displacement is approximated by first-order Lagrangepolynomials, function ˆ u ( x , t ) in ˆ U ⊂ U . This choice of elements results in continuityof normal velocity across the elements facets and hence guarantees discrete mass con-servation.The discretized state functions can be written as follows (we adopt some of the no-tation from (Ferronato et al., 2010)): u ( x , t ) ≈ ˆ u ( x , t ) = ˆ u x ( x , t )ˆ u y ( x , t )ˆ u z ( x , t ) –25– (cid:80) n n i =1 φ i ( x ) u x,i ( t ) (cid:80) n n i =1 φ i ( x ) u y,i ( t ) (cid:80) n n i =1 φ i ( x ) u z,i ( t ) = B u ( x ) U ( t ) . where B u ( x ) is given by: B u ( x ) = φ ( x ) ... φ n n ( x ) 0 00 φ ( x ) ... φ n n ( x ) 00 0 φ ( x ) ... φ n n ( x ) , U ( t ) is a vector of the DOFs of the displacement in x , followed by the y -displacementDOFs and then by the z -displacement DOFs. The integer n n is the number of nodes.The basis functions φ i are given by the formula: φ i ( x ) = (cid:40) ( x − x l )( x i − x l )( x i − x l )( x i − x l ) , for x ∈ T ( l ) l ∈ { j i, , j i, , ..., j i,n φi } , for x ∈ Ω − ∪ l ∈{ j i, ,j i, ,...,j i,nφi } T ( l ) , where x i is the position vector of the node i associated with the basis function φ i and x l is the position vector of the normal projection of x i on the plane containing the faceopposite to the node x i in the tetrahedron T ( l ) . The integer n φ i is the number of thetetrahedra that share the node i and the set of indices { j i, , j i, , ..., j i,n φi } ⊂ { , , ..., n e } is the global indices of those tetrahedra, where n e is the number of tetrahedra in the mesh.The piecewise constant approximation of the pressure can be written as p ( x , t ) ≈ ˆ p ( x , t ) = (cid:80) n e j =1 χ j ( x ) p j ( t ) = B p ( x ) T P ( t ), B p ( x ) is the vector of the basis functions χ j ( x ), P ( t )is the DOFs of the pressure, and the basis functions χ j ( x ) are given by: χ j ( x ) = (cid:26) , for x ∈ T ( j ) , for x ∈ Ω − T ( j ) , where T ( j ) is the tetrahedron associated to the j th basis function, χ j ( x ). The lowest or-der RV approximation of the velocity is given by the following expression: q ( x , t ) = (cid:80) n f k =1 ψ k ( x ) q k ( t ) = B q ( x ) T Q ( t ), where n f is the number of faces, B q ( x ) is the vector of the basis functions,and Q ( t ) is the vector of the velocity DOFs. The vector-valued basis functions ψ k ( x )are given by: ψ k ( x ) = (cid:40) ± ( x − x k )3 | V ( T ( m ) ) | , for x ∈ T ( m ) m ∈ { j k, , j k, } , for x ∈ Ω − ∪ m ∈{ j k, ,j k, } T ( m ) , where { j k, , j k, } ⊂ { , , ..., n e } are the indices of the tetrahedrons that share the k th face, and x k is the position of the node that is not shared by k th face in the tetrahedron T ( m ) . The conventional choice of the sign is such that the vector ψ k ( x ) points outwardthe element T ( m ) with smallest index m .We substitute the approximated functions in the weak form (A7)–(A9) to obtainthe discretized system. The resulting linear system of equations that we need to solveeach time step is as follows: MP k + α DU k + ∆ t k D q Q k =∆ t k F pk + MP k − + α DU k − (A10) − α GP k − EU k = − F u k − F u ,bcsk (A11)∆ t k G p P k − ∆ t k KQ k =∆ t k F p,bcsk , (A12)where the subscript k denote that the vector is evaluated (or to be evaluated) at time t k . We combine the equations (A10)–(A12) in one system as follows: L k X k = F k + NX k − , (A13)where L k = M α D ∆ t k D q − α G − E 0 − ∆ t k G p ∆ t k K , –26– = M α D 00 0 00 0 0 , and F k = ∆ t k F pk − F u k − F u ,bcsk − ∆ t k F p,bcsk . X k = P k U k Q k . For derivation purposes, we combine solving for all the N t time steps, system (A13), inone “global” system (following Hesse and Stadler (2014)) which we write as: ¯S ¯X = ¯F , (A14)where ¯S = L ... − N L ... : : : ... : : ... L N t −
00 0 0 ... − N L N t , ¯X = X X : X N t − X N t , and ¯F = F F : F N t − F N t . The bar is used in the notation to distinguish the space-time discretized operators andvectors from the space-only discretized operators.
Appendix B Derivation of Newton Iteration for Estimating the MAPPoint
We form the Lagrangian for the constrained optimization problem (17) to be: L ( m , ¯X , ¯Y ) = 12 ( ¯B ¯X − d obs ) T Γ − ( ¯B ¯X − d obs ) (B1)+ 12 ( m − ¯ m ) T R ( m − ¯ m ) + ¯Y T ( ¯S ( m ) ¯X − ¯F ) . (B2)The adjoint equation, the state equation and the gradient are given by: ∂ L ∂ ¯X ( ¯X , ¯Y , m ) = ¯B T Γ − ( ¯B ¯X − d obs ) + ¯S T ( m ) ¯Y = 0 (adjoint) , (B3) ∂ L ∂ ¯Y ( ¯X , ¯Y , m ) = ¯S ( m ) ¯X − ¯F = 0 (state) , (B4) ∂ L ∂ m ( ¯X , ¯Y , m ) = R ( m − ¯ m ) + ¯C T ( m ) ¯Y (gradient) . (B5)Consider the Lagrangian L H formed as follows: L H := δ ¯X T (cid:0) ¯B T Γ − ( ¯B ¯X − d obs ) + ¯S T ( m ) ¯Y (cid:1) (B6)+ δ ¯Y T (cid:0) ¯S ( m ) ¯X − ¯F (cid:1) (B7)+ δ m T (cid:0) R ( m − ¯ m ) + ¯C T ( m ) ¯Y (cid:1) , (B8)where δ ¯X , δ ¯Y and δ m are incremental directions. The incremental forward problem isdefined as follows: ∂ L H ∂ ¯Y = ¯S ( m ) δ ¯X + ¯C ( m ) δ m = 0 . (B9)And similarly the incremental adjoint problem is: ∂ 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 . (B10) –27– he k th time step in the incremental forward problem will require solving the system: L ( m ) δ X k = N δ X k − − C k ( m ) δ m , (B11)and the k th time step in the incremental adjoint problem requires solving the followingsystem: L T ( m ) δ Y k = N T δ Y k +1 − W XX δ X k − W k Xm ( m ) δ m . (B12)The Hessian action on δ m is given by: ∂ L H ∂ m = ∂∂ m ( δ ¯X T ¯S T ( m ) ¯Y ) + ∂∂ m ( δ ¯Y T ¯S ( m ) ¯X ) + R δ m + ∂∂ m ( δ m T ¯C T ( m ) ¯Y )= W m ¯X ( m ) δ ¯X + ¯C T ( m ) δ ¯Y + R δ m + R mm ( m ) δ m . (B13)Substituting the values of δ X and δ X from the systems (B9) and (B10) into the expres-sion (B13) gives the Hessian expression: H = R + 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 ) . (B14) Acknowledgments
This work was supported by National Science Foundation (NSF) Grants CBET1508713and ACI-1550593, and DOE grant DE-SC0019303. A.A. would like to acknowledge fund-ing from the Ministry of Education in Saudi Arabia. The authors would like to thankDr. Thomas Burbey for providing the Nevada experiment GPS and well head data andfor his helpful comments. The authors would also like to thank Dr. Umberto Villa, a maindeveloper of hIPPYlib , for the many helpful discussions during the course of this research,and Dr. Jeonghun (John) Lee for his help with the convergence verification of the 3D3-field poroelasticity forward solver to a manufactured solution. The Envisat SAR im-agery used in this study can be downloaded through the UNAVCO Data Center SARarchive.
References
Alexanderian, A., Petra, N., Stadler, G., & Ghattas, O. (2014). A-optimal designof experiments for infinite-dimensional Bayesian linear inverse problems withregularized (cid:96) -sparsification. SIAM Journal on Scientific Computing , (5),A2122–A2148. doi: 10.1137/130933381Alexanderian, A., Petra, N., Stadler, G., & Ghattas, O. (2016). A fast and scalablemethod for A-optimal design of experiments for infinite-dimensional Bayesiannonlinear inverse problems. SIAM Journal on Scientific Computing , (1),A243–A272. doi: 10.1137/140992564Alexanderian, A., Petra, N., Stadler, G., & Ghattas, O. (2017). Mean-variance risk-averse optimal control of systems governed by PDEs with random parameterfields using quadratic approximations. SIAM/ASA Journal on UncertaintyQuantification , (1), 1166–1192. (arXiv preprint arXiv:1602.07592) doi:10.1137/16M106306XAmelung, F., Galloway, D. L., Bell, J. W., Zebker, H. A., & Laczniak, R. J. (1999).Sensing the ups and downs of Las Vegas: InSAR reveals structural control ofland subsidence and aquifer-system deformation. Geology , (6), 483-486.Bashir, O., Willcox, K., Ghattas, O., van Bloemen Waanders, B., & Hill, J. (2008).Hessian-based model reduction for large-scale systems with initial conditioninputs. International Journal for Numerical Methods in Engineering , ,844-868. –28– ell, R. E. (2008). The role of subglacial water in ice-sheet mass balance. NatureGeoscience , (5), 297–304. doi: 10.1038/ngeo186Biot, M. A. (1941). General Theory of Three- Dimensional Consolidation. Journal ofApplied Physics . doi: 10.1063/1.1712886Bohling, 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., 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. (Gordon Bell Prize finalist)Bui-Thanh, T., & Ghattas, O. (2012a). Analysis of the Hessian for inverse scatteringproblems. Part II: Inverse medium scattering of acoustic waves.
Inverse Prob-lems , (5), 055002. doi: 10.1088/0266-5611/28/5/055002Bui-Thanh, T., & Ghattas, O. (2012b). Analysis of the Hessian for inverse scatteringproblems. Part I: Inverse shape scattering of acoustic waves. Inverse Problems , (5), 055001. doi: 10.1088/0266-5611/28/5/055001Bui-Thanh, T., & Ghattas, O. (2013). Analysis of the Hessian for inverse scatter-ing problems. Part III: Inverse medium scattering of electromagnetic waves. Inverse Problems and Imaging , (4), 1139–1155.Bui-Thanh, T., & Ghattas, O. (2015). A scalable MAP solver for Bayesian inverseproblems with Besov priors. Inverse Problems and Imaging , (1), 27-54.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/WR022i002p00211Chaussard, E., B¨urgmann, R., Shirzaei, M., Fielding, E. J., & Baker, B. (2014).Predictability of hydraulic head changes and characterization of aquifer-system and fault properties from InSAR-derived ground deformation. Journal of Geophysical Research: Solid Earth , (8), 6572–6590. doi:10.1002/2014JB011266Chen, J., Knight, R., & Zebker, H. A. (2017). The temporal and spatial variabil-ity of the confined aquifer head and storage properties in the san luis valley,colorado inferred from multiple insar missions. Water Resources Research , (11), 9708-9720. doi: 10.1002/2017WR020881Chen, J., Knight, R., Zebker, H. A., & Schreder, W. A. (2016). Confined aquiferhead measurements and storage properties in the san luis valley, colorado, fromspaceborne insar observations. Water Resources Research , (5), 3623–3636.doi: 10.1002/2015WR018466 –29– hen, P., Villa, U., & Ghattas, O. (2017). Hessian-based adaptive sparse quadraturefor infinite-dimensional Bayesian inverse problems. Computer Methods in Ap-plied Mechanics and Engineering , , 147-172. Retrieved from https://doi.org/10.1016/j.cma.2017.08.016 Cherry, J., & Freeze, A. (1979).
Groundwater (1st ed.). Prentice Hall.Crestel, B., Alexanderian, A., Stadler, G., & Ghattas, O. (2017). A-optimalencoding weights for nonlinear inverse problems, with application to theHelmholtz inverse problem.
Inverse Problems , (7), 074008. Retrievedfrom http://iopscience.iop.org/10.1088/1361-6420/aa6d8e Deutsch, C. V., & Journel, A. G. (1998).
GSLIB: geostatistical software library anduser’s guide. Second edition .Du, J., & Olson, J. E. (2001, 9). A poroelastic reservoir model for pre-dicting subsidence and mapping subsurface pressure fronts.
Journal ofPetroleum Science and Engineering , (3-4), 181–197. Retrieved from http://linkinghub.elsevier.com/retrieve/pii/S0920410501001310 doi: 10.1016/S0920-4105(01)00131-0Eaton, 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.002Eisenstat, S. C., & Walker, H. F. (1996). Choosing the forcing terms in an inexactNewton method. SIAM Journal on Scientific Computing , (1), 16–32. Re-trieved from http://dx.doi.org/10.1137/0917003 doi: 10.1137/0917003Ferronato, 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 quantificationin large-scale linear inverse problems based on low-rank partial Hessian ap-proximations. SIAM Journal on Scientific Computing , (1), 407-432. doi:10.1137/090780717Galloway, D., & Hoffmann, J. (2007). The application of satellite differential SARinterferometry-derived ground displacements in hydrogeology. HydrogeologyJournal , (1), 133-154. doi: 10.1007/s10040-006-0121-5Galloway, D. L., & Burbey, T. J. (2011). Review: regional land subsidence accompa-nying groundwater extraction. Hydrogeology Journal , (8), 1459–1486.Geuzaine, C., & Remacle, J.-F. (2009). Gmsh: A 3-D finite element mesh generatorwith built-in pre- and post-processing facilities. International Journal for Nu-merical Methods in Engineering , (11), 1309-1331. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2579 doi: 10.1002/nme.2579Goovaerts, P. (1997). Geostatistics for natural reources evaluation .Gunzburger, M. D. (2003).
Perspectives in flow control and optimization . Philadel-phia: SIAM.Haga, 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.Hanasaki, N., Kanae, S., Oki, T., Masuda, K., Motoya, K., Shirakawa, N., . . .Tanaka, K. (2008). An integrated model for the assessment of global waterresources - Part 2: Applications and assessments. Hydrology and Earth SystemSciences , (4), 1027–1037. doi: 10.5194/hess-12-1027-2008Hanssen, R. (2001). Radar interferometry: Data interpretation and error analysis .Springer. –30– einkenschloss, M. (1993). Mesh independence for nonlinear least squares problemswith norm constraints.
SIAM Journal on Optimization , , 81–117.Hesse, M., & Stadler, G. (2014). Joint inversion in coupled quasistatic poroelasticity. Journal of Geophysical Research: Solid Earth , (2), 1425–1445.Hoekstra, A. Y., & Mekonnen, M. M. (2012). The water footprint of humanity. Proceedings of the National Academy of Sciences , (9), 3232–3237. Re-trieved from doi:10.1073/pnas.1109936109Hoffmann, J., Zebker, H. A., Galloway, D. L., & Amelung, F. (2001). Seasonalsubsidence and rebound in Las Vegas Valley, Nevada, observed by syntheticaperture radar interferometry. Water Resources Research , (6), 1551–1566.doi: 10.1029/2000WR900404Iglesias, M. A., & McLaughlin, D. (2012). Data inversion in coupled subsurface flowand geomechanics models. Inverse Problems , (11), 115009.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.047Kaipio, J., & Somersalo, E. (2005). Statistical and computational inverse prob-lems (Vol. 160). Springer-Verlag New York. Retrieved from https://link.springer.com/978-0-387-27132-3 doi: 10.1007/b138659Lindgren, 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. (2012). Automated solution of differentialequations by the finite element method: The FEniCS book (Vol. 84). SpringerScience & Business Media.Lu, Z., & Danskin, W. R. (2001). InSAR analysis of natural recharge to define struc-ture of a ground-water basin, San Bernardino, California.
Geophysical ResearchLetters , (13), 2661–2664.Mariethoz, G., & Caers, J. (2014). Multiple-point Geostatistics: Stochastic Modelingwith Training Images . doi: 10.1002/9781118662953Martin, 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/96WR00160Miller, M. M., Shirzaei, M., & Argus, D. (2017). Aquifer Mechanical Properties andDecelerated Compaction in Tucson, Arizona. Journal of Geophysical Research:Solid Earth , (10), 8402–8416. doi: 10.1002/2017JB014531Nocedal, J., & Wright, S. J. (2006). Numerical optimization (second ed.). Berlin,Heidelberg, New York: Springer Verlag.Oliver, 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. –31– hillips, P. J. (2005). Finite element methods in linear poroelasticity: theoretical andcomputational results (Unpublished doctoral dissertation).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.Remy, N., Boucher, A., & Wu, J. (2009). Applied Geostatistics with SGeMS .Cambridge: Cambridge University Press. Retrieved from http://ebooks.cambridge.org/ref/id/CBO9781139150019 doi: 10.1017/CBO9781139150019Rosen, P., Hensley, S., Joughin, I., K.Li, F., Madsen, S., Rodriguez, E., & Goldstein,R. M. (2000). Synthetic aperture radar interferometry.
Proceedings of theIEEE , (3), 333-382. doi: 10.1109/5.838084Rost, S., Gerten, D., Bondeau, A., Lucht, W., Rohwer, J., & Schaphoff, S.(2008). Agricultural green and blue water consumption and its influenceon the global water system. Water Resources Research , (9), 1–17. doi:10.1029/2007WR006331Scanlon, B. R., Faunt, C. C., Longuevergne, L., Reedy, R. C., Alley, W. M.,McGuire, V. L., & McMahon, P. B. (2012). Groundwater depletion andsustainability of irrigation in the us high plains and central valley. Proceedingsof the national academy of sciences , (24), 9320–9325.Schmidt, D. A., & Br¨ugmann, R. (2003). Time-dependent land uplift and sub-sidence in the Santa Clara Valley, California, from a large interferometricsynthetic aperture radar data set. Journal of Geophysical Research: SolidEarth , (B9).Segall, P., Rubin, A. M., Bradley, A. M., & Rice, J. R. (2010). Dilatant strengthen-ing as a mechanism for slow slip events. Journal of Geophysical Research: SolidEarth , (12), 1–37. doi: 10.1029/2010JB007449Showalter, R. E. (2000). Diffusion in poro-elastic media. Journal of MathematicalAnalysis and Applications , (1), 310–340.Stuart, A. M. (2010). Inverse problems: A Bayesian perspective. Acta Numerica , ,451-559. doi: 10.1017/S0962492910000061Tarantola, A. (2005). Inverse problem theory and methods for model parameter esti-mation . Philadelphia, PA: SIAM.UNESCO. (2012).
World Water Development Report Volume 4: ManagingWater under Uncertainty and Risk (Vol. 1). Retrieved from doi: 10.1608/FRJ-3.1.2Vasco, D., Karasaki, K., & Kishida, K. (2001). A coupled inversion of pressure andsurface displacement.
Water Resources Research , (12), 3071–3089. Retrievedfrom http://onlinelibrary.wiley.com/doi/10.1029/2001WR000391/full Vasco, D. W., Ferretti, A., & Novali, F. (2008, 11). Estimating permeability fromquasi-static deformation: Temporal variations and arrival-time inversion.
Geo-physics , (6), O37-O52. Retrieved from http://library.seg.org/doi/abs/10.1190/1.2978164 doi: 10.1190/1.2978164Villa, U., Petra, N., & Ghattas, O. (2018). hIPPYlib: an Extensible Software Frame-work for Large-scale Deterministic and Bayesian Inverse Problems. Journal ofOpen Source Software , (30). doi: 10.21105/joss.00940Villa, U., Petra, N., & Ghattas, O. (2019). hIPPYlib: An extensible softwareframework for large-scale inverse problems governed by PDEs; Part I: Deter-ministic inversion and linearized Bayesian inference. Submitted . Retrieved from http://arxiv.org/abs/1909.03948
Vorosmarty, C. J. (2000, 7). Global Water Resources: Vulnerability from ClimateChange and Population Growth.
Science , (5477), 284–288. Retrieved from doi: 10.1126/science.289.5477.284 –32– ackernagel, H. (2010). Multivariate Geostatistics: An Introduction with Appli-cations . Springer. Retrieved from
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), 1–5. doi: 10.1029/2010GL044571Wang, H. (2000). Theory of Linear Poroelasticity with Applications to Geomechanicsand Hydrogeology . Princeton, NJ: Princeton University Press.Yeh, W. W.-G. (1986). Review of parameter identification procedures in ground-water hydrology: The inverse problem.
Water Resources Research , (2), 95–108.Zebker, H. A., & Villasenor, J. (1992, Sep.). Decorrelation in interferometric radarechoes. IEEE Transactions on Geoscience and Remote Sensing , (5), 950-959. doi: 10.1109/36.175330(5), 950-959. doi: 10.1109/36.175330