Two-grid method on unstructured tetrahedra: Applying computational geometry to staggered solution of coupled flow and mechanics problems
NNoname manuscript No. (will be inserted by the editor)
Two-grid method on unstructured tetrahedra:Applying computational geometry to staggeredsolution of coupled flow and mechanics problems
Saumik Dana · Xiaoxi Zhao · Birendra Jha
Received: date / Accepted: date
Abstract
We develop a computational framework that leverages the features ofsophisticated software tools and numerics to tackle some of the pressing issues inthe realm of earth sciences. The algorithms to handle the physics of multiphaseflow, concomitant geomechanics all the way to the earth’s surface and the complexgeometries of field cases with surfaces of discontinuity are stacked on top of eachother in a modular fashion which allows for easy use to the end user. The currentfocus of the framework is to provide the user with tools for assessing seismicrisks associated with energy technologies as well as for use in generating forwardsimulations in inversion analysis from data obtained using GPS and InSAR. In thiswork, we focus on one critical aspect in the development of the framework: theuse of computational geometry in a two-grid method for unstructured tetrahedralmeshes
Keywords
Energy technologies · Coupled multiphase flow and geomechanics · Computational geometry · Staggered solution algorithm
The development of a computationally inexpensive framework with the capabilityto model deformation of the earth’s surface and fault slip due to deep subsurface
Saumik DanaUniversity of Southern California, Los Angeles, CA 90007E-mail: [email protected] ZhaoUniversity of Southern California, Los Angeles, CA 90007E-mail: [email protected] JhaUniversity of Southern California, Los Angeles, CA 90007E-mail: [email protected] a r X i v : . [ c s . C E ] F e b Saumik Dana et al. overburdensideburden non-pay free surface d H l R l R l G l G reservoir l G > l R d >> H overburden pressure??uplift/subsidence?? Fig. 1: Direct imposition of heuristic overburden pressures on the flow domaincompletely disregards the mechanical behavior of the surrounding rock and obvi-ates a study of fault slip away from the reservoir as well as deformation of earth’ssurfacepressure perturbations associated with multiphase flow is critical from two differentstandpoints – Assessing the seismic risks associated with carbon sequestration, enhancedgeothermal systems, waste water disposal, enhanced oil recovery thereby of-fering guiding protocols for the design of these operations [19,43,52,47,39,54,28,61] – Estimating subsurface properties using inversion analysis of ground deforma-tion data obtained from Global Positioning System (GPS) and interferometricsynthetic aperture radar (InSAR). This typically requires multiple forward sim-ulations with the geomechanical domain extending all the way to the earth’ssurface [27,26,57,55,56,18,60,33,31,34]The major bugbear to developing such a framework is the computationally in-tractable size of the geomechanical grid no matter which numerical method isused to resolve the coupled system of equations. The typical approaches to avoidsuch issues are to impose an overburden pressure directly on the reservoir thustreating it as a coupled problem domain (see Fig. 1) or to model flow on hugedomain with zero permeability cells mimicking the no flow boundary conditionbetween the flow and non-flow region. The former approach precludes a study ofsurface deformation, does not mimic the true effect of the overburden on the stresssensitive reservoir, and is incapable of capturing induced seismicity inside faultsaway from reservoir whereas the latter approach is computationally intractablefor large field scale problems due to memory requirements. In order to addressthis, we develop a two-grid coupled multiphase flow and geomechanics framework which allows for spatial decoupling of the flow and geomechanics domains withthe geomechanics subproblem being resolved on a separate grid with a larger spa-tial extent going all the way to the free surface (see Fig. 2). This computationalframework is built on top of a staggered solution algorithm that solves the flowand mechanics subproblems sequentially and iteratively. itle Suppressed Due to Excessive Length 3 coupled problem free surface g Solve geomechanics on geomechanics mesh Solve flow on flow mesh project project
Fig. 2: The two-grid method enables the study of induced seismicity in faults awayfrom reservoir, takes into account the mechanical behavior of surrounding rockand allows for determination of surface deformation due to subsurface pressureperturbationTypically, in such problems, the geomechanics mesh is expected to be coarserthan the flow mesh everywhere the two meshes exist, but we generalize the methodto the cases where the geomechanics elements can be smaller than the flow ele-ments in small localized region where capturing the mechanics is more pertinent.Furthermore, we would like to generalize the notion of the two-grid from struc-tured hexahedral meshes (as was done in [22]) to unstructured tetrahedral meshes.A depiction of the intersection of two tetrahedral elements in given in Fig. 3. Inthis document, we focus on demonstrating the convergence of the two-grid methodfor unstructured tetrahedral grids using the classical Mandel’s problem [44,5] an-alytical solution. We work two cases, one in which we have a finer mesh for flow,and the other in which we have a finer mesh for mechanics.This paper is structured as follows: the governing equations are provided insection 2, the solution strategy is explained in section 3, the salient features ofour software are elucidated in section 4, the numerical simulations for Mandel’sproblem are provided in section 5, and the conclusions and outlook are providedin section 6.
We use a classical continuum representation in which the fluids and the solidskeleton are viewed as overlapping continua [10,21]. The governing equations forcoupled flow and geomechanics are obtained from conservation of mass and bal-ance of linear momentum. We assume that the deformations are small, that thegeomaterial is isotropic, and that the conditions are isothermal.
Saumik Dana et al. E p E f E f E p Fig. 3: Depiction of couple of intersecting flow-geomechanics element pairs. Werefer to the flow element as E f and geomechanics element as E p . Depiction is toreiterate the point that there is no restriction of whether the flow or geomechanicselement needs to be smaller than the otherUnder the quasistatic assumption for earth displacements, the governing equa-tion for linear momentum balance of the solid/fluid system can be expressed as ∇ · σ + ρ b g = , (2.1)where σ is the Cauchy total stress tensor, g is the gravity vector, and ρ b = φ (cid:80) n phase β ρ β S β + (1 − φ ) ρ s , is the bulk density, ρ β and S β are the density andsaturation of fluid phase β , ρ s is the density of the solid phase, φ is the trueporosity, and n phase is the number of fluid phases. The true porosity is defined asthe ratio of the pore volume V p to the bulk volume V b in the current (deformed)configuration.Assuming that the fluids are immiscible, the mass-conservation equation foreach phase α is dm α dt + ∇ · w α = ρ α f α , (2.2)where the accumulation term dm α /dt describes the time variation of fluid massrelative to the motion of the solid skeleton, w α is the mass-flux of fluid phase α relative to the solid skeleton, and f α is the volumetric source term for phase α .Balance equations (2.1) and (2.2) are coupled by virtue of poromechanics. On onehand, changes in the pore fluid pressure lead to changes in effective stress, and induce deformation of the porous material—such as ground subsidence caused bygroundwater withdrawal. On the other hand, deformation of the porous mediumaffects fluid mass content and fluid pressure. The simplest model of this two-waycoupling is Biot’s macroscopic theory of poroelasticity [11,29,20]. The specificequations of single phase poromechanics are provided in Appendix A itle Suppressed Due to Excessive Length 5 New time stepNew coupling iterationSolve flow on flow grid if not convergedif convergedkeeping rate of volumetric stress fixedProject pressures onto geomechanics gridSolve geomechanics on geomechanics gridProject volumetric strains onto flow gridCheck for convergence
Fig. 4: Two grid staggered solution algorithmThe two-grid approach is built on top of a staggered solution algorithm (Fig.4) in which the flow and geomechanics subproblems are solved sequentially anditeratively using a fixed stress split iterative scheme [36,37,46,45,17,59,16,24,25,23,7,8,14,58,13,53,51,6,42].3.1 Space and time discretizationWe use the finite volume method for the discretization of the flow problem [9],and the nodal-based finite element method for the discretization of the mechanicsproblem [32,62]. This space discretization is locally mass conservative at the ele-ment level, and enjoys excellent stability properties [35,49,50,38]. The pressuresand saturations degrees of freedom are located at the element centers of the flowgrid, and the displacement vector degrees of freedom are located at the elementnodes of the mechanics grid. In quasi-static poromechanics, the time derivativeappears only in the accumulation term of the fluid mass balance equation, and wetreat it using a fully implicit Backward Euler time integration scheme.
Saumik Dana et al.
Algorithm 1
Determining the relationship between flow-geomechanics elementpairs for
Loop over all geomechanics elements dofor
Loop over all flow elements do count1 ← for Loop over four vertices of flow element doif
Vertex inside geomechanics element then count1 ← count1 + 1 if count1!=0 then Flow to geomechanics projection operator ← Meas ( E f )Geomechanics to flow projection operator ← Meas ( E p )count2 ← for Loop over four vertices of geomechanics element doif
Vertex inside flow element then count2 ← count2 + 1 if count2!=0 then Flow to geomechanics projection operator ← Meas ( E f )Geomechanics to flow projection operator ← Meas ( E p ) from the geomechanical domain to the flow domain. The transfered pore pressureat every geomechanics element is a volume average of the pore pressures at theflow elements interacting with the geomechanics element. Likewise, the transferedvolumetric strain at each flow element is a volume average of the volumetric strainsat the geomechanics elements interacting with the flow element. The global two-grid operators are constructed in the pre-processing step by applying a local two-grid procedure to each pair of intersecting finite elements of the two grids as shownin Algorithm 1. The relationship between the flow-geomechanics element pair isobtained by determining the barycentric coordinates of the vertices of the flowelement with respect to the geomechanics element and vice versa as explained inAlgorithm 1. We developed a coupled multiphase flow and geomechanical simulator by integrat-ing our version of Stanford’s General Purpose Research Simulator (GPRS) [15,48]called flowsim which operates as the flow simulator into an open source frameworkcalled PyLith [1,3] which operates as the mechanics simulator. The workflow ofour software are elucidated in Fig. 5. Available seismic and geological informationare used to construct a geomechanical model and populate it with poromechanicalproperties. Reservoir flow properties (porosity and permeability) and fluid prop-erties (fluid density, viscosity, and compressibility) are obtained from an olderuncoupled model and transferred to the geomechanical model using the MatlabReservoir Simulation Toolbox [41]. , itle Suppressed Due to Excessive Length 7 Seismic andgeological informationReservoir flow properties,multiphase fluid properties Reservoir flowmodelGeomechanicalmodel Displacementsand stressFluid pressuresand saturationsflowsimgeomechanics grid finite volume grid vtk .outflow grid hdf5MRST
Fig. 5: WorkflowN , and water) flow through subsurface. It treats element connections through ageneral connection list, which allows for both structured and unstructured grids.It is capable of handling complex production and injection scenarios in the field,such as wells perforated at multiple depths and flowing under variable rate andpressure controls.4.2 Mechanics SimulatorPyLith is a open-source finite element code for the simulation of static and dy-namic large-scale deformation problems [2,4]. It uses an implicit formulation tosolve quasi-static problems and an explicit formulation to solve dynamic ruptureproblems. The salient features of Pylith are elucidated in Fig. 6. Some of the ad- vantages of PyLith are: (a) It is written using C++ and Python languages and isextendable, (b) it is suitable for parallel computing, (c) it allows localized deforma-tion along discrete features, such as faults, (d) it is well integrated with meshingcodes, such as LaGriT for tetrahedral meshes [40] and CUBIT for both tetrahedraland hexahedral meshes [12] Saumik Dana et al.
Fig. 6: PyLith requires a finite-element mesh, simulation parameters, and spatial databases (defining the spatial variation of various parameters). It writes the so-lution output to either VTK or HDF5/Xdmf files, which can be visualized withParaView or Visit. Post-processing is generally done using the HDF5 files withPython or Matlab scripts [2]4.3 GriddingWe employ different grids for flowsim and PyLith to go with our two-grid method.The finite element grids are generated using CUBIT [12] mesh generation softwareand exported in the Exodus-II format. Since the numerical method employed forthe flow model is finite volume instead of finite element, we employ some furtherprocessing of the flow grid file to render it for the finite volume calculations. Weprocess the Exodus-II grid file using a MATLAB script to generate the equivalentfinite volume grid in the domain with element centroid coordinates, element bulkvolumes, and face transmissibilities in the Corner Point Geometry format [30]4.4 Coupling
We design a C++ class, iflowsim, to allow communication between the flow and themechanics simulators. The two-grid operators built on top of iflowsim to projectinformation between the two grids. PyLith supports distributed memory paral-lelization (Message Passing Interface or MPI) whereas flowsim’s parallelization isbased on the shared memory architecture (Multiprocessing or OpenMP) itle Suppressed Due to Excessive Length 9
No-flow No-flowNo-flow p=0 σ = 2 MP a (0,0) (50,0)(50,10)(0,10)
Fig. 7: Schematic for Mandel’s problem (a) Flow mesh, 5706 elements (b) Pylith mesh, 3697 elements
Fig. 8: Finer mesh for flow compared to Pylith (a) Pylith mesh, 5706 elements (b) Flow mesh, 3697 elements
Fig. 9: Finer mesh for Pylith compared to flow
Numerical: flow 5706, pylith 3697Analytical (a) Flow mesh: 5706 elements and pylith mesh: 3697 elements
Numerical: flow 3697, pylith 5706Analytical (b) Flow mesh: 3697 elements and pylith mesh: 5706 elements
Fig. 10: Comparison of non-monotonic pore pressure responseThe Mandel’s problem [44,5] has been used as a benchmark problem for testingthe validity of numerical codes of coupled poroelasticity. Its main feature, the itle Suppressed Due to Excessive Length 11
Mandel-Cryer effect, is that the pore pressure at the center of a loaded specimenrises above its initial value because of the two-way coupling between fluid flowand solid deformation. It involves a long specimen of rectangular cross sectionpressed on one side with an impermeable plate that applies a constant compressivestress and fixed on two sides using impermeable roller boundaries (Fig. 7). Thefourth side of the cross section is free from normal and shear stresses (traction-free boundary) and is open to the atmosphere (constant pressure boundary). Theporous medium is saturated with a slightly compressible fluid, water, with initialpressure set at the reference value. Since the specimen is long, we assume planestrain conditions, namely, that the displacement and fluid flux vanish in the zdirection (perpendicular to the 2-D domain).With these boundary conditions, the three-dimensional equations of poroelas-ticity reduce to one-dimensional equations for σ xx ( y, t ) and p ( y, t ), which can besolved analytically [44,5]. At t = 0 + , a uniform undrained pressure is generatedby the Skempton effect, along with uniform stress σ xx = − σ . The specimen ex-pands toward the top boundary due to the Poisson effect. As time progresses, thepressure near the top boundary decreases because of fluid drainage, which makesthe specimen more compliant there. If the hydraulic diffusivity is small, the effectof drainage is not observed immediately near the no-flux bottom boundary. Thisresults into load transfer of compressive total stress toward the bottom bound-ary, in response to which the pressure there continues to rise above its undrainedvalue. At long times, all excess pressure vanishes and a uniform horizontal stress, σ xx = − σ , returns. Hence, the pressure evolution at points away from the drainedboundary is nonmonotonic, a phenomenon not observed in a purely diffusive pro-cess such as that modeled by the Terzaghi theory, where the pressure is uncoupledfrom the solid deformation.In one simulation, we employ a finer mesh for the flow problem as shown inFig. 8 whereas in the other simulation, we employ a finer mesh for the mechanicsproblem as shown in Fig. 9. Fig. 10 compares the pressure response at the rightbottom corner ( x = 50 , y = 0) at different times. We developed a two-grid method which leverages computational geometry on un-structured tetrahedral meshes that allows us to solve the flow and geomechanicssubproblems on separate grids. This allows use for flexibility in the choice of thetwo separate meshes depending on the needs of the specific problem. There isno restriction on whether the flow element needs to be finer or coarser than thegeomechanics element in any region of the subdomains. This method also allowsus to have separate domains for flow and geomechanics, which is critical for largescale field problems typically encountered in energy technologies. We shall now employ our computational framework to solve field scale problems with faults infuture articles.
Acknowledgements
The first author would like to thank Sreekanth Arikatla (Staff R&DEngineer, Kitware Inc.) for discussions on computational geometry for unstructured tetrahedra2 Saumik Dana et al.
Conflict of interest
The authors declare that they have no conflict of interest.
A Single-phase poromechanics
For isothermal single-phase flow of a slightly compressible fluid in a poroelastic medium withno stress dependence of permeability, the single-phase fluid mass conservation equation reducesto dmdt + ∇ · w = ρ f f, (A.1)where m is the fluid mass content (fluid mass per unit bulk volume of porous medium), ρ f isthe fluid density, w = ρ f v is the fluid mass flux (fluid mass flow rate per unit area and time),and v is the seepage velocity relative to the deforming skeleton, given by Darcy’s law: v = − k µ (cid:0) ∇ p − ρ f g (cid:1) , (A.2)where k is the intrinsic permeability tensor, µ is the fluid dynamic viscosity and p is thepore-fluid pressure [10]. It is useful to define the fluid content variation ζ , ζ := δmρ f, , (A.3)where δm = m − m is the increment in fluid mass content with respect to the initial referencestate, and ρ f, is the reference fluid density. The self-consistent theory of poroelastic behaviorproposed by [11] links the changes in total stress and fluid pressure with changes in strain andfluid content. Following [20], the poroelasticity equations can be written in incremental formas δ σ = C dr : ε − bδp ζ = bε v + 1 M δp, (A.4)where C dr is the rank-4 drained elasticity tensor, is the rank-2 identity tensor, ε is thelinearized strain tensor, defined as the symmetric gradient of the displacement vector u , ε := 12 (cid:16) ∇ u + ∇ T u (cid:17) , (A.5)and ε v = tr( ε ) is the volumetric strain. Note that we use the convention that tensile stress ispositive. It is useful to express the strain tensor as the sum of its volumetric and deviatoriccomponents: ε = 13 ε v + e , (A.6)from which it follows that the volumetric stress σ v = tr( σ ) / δσ v = K dr ε v − bδp. (A.7)Equation (A.4) implies that the effective stress in single-phase poroelasticity, responsible forskeleton deformation, is defined in incremental form as δ σ (cid:48) := δ σ + bδp . (A.8)Biot’s theory of poroelasticity has two coupling coefficients: the Biot modulus M and the Biotcoefficient b . They are related to rock and fluid properties as [20]1 M = φ c f + b − φ K s , b = 1 − K dr K s , (A.9)itle Suppressed Due to Excessive Length 13where c f = 1 /K f is the fluid compressibility, K f is the bulk modulus of the fluid, K s isthe bulk modulus of the solid grain, and K dr is the drained bulk modulus of the porousmedium. To set the stage for the numerical solution strategy of the coupled problem, it isuseful to write the fluid mass balance equation (A.1) (the pressure equation) in a way thatexplicitly recognizes the coupling with mechanical deformation. Equations (A.4) state that theincrement in fluid mass content has two components: increment due to expansion of the porespace and increment due to increase in the fluid pressure. Assuming small elastic deformationsand applying linearization from the reference state to the current state, we can write Eqs. (A.4)as σ − σ = C dr : ε − b ( p − p ) , (A.10)1 ρ f, ( m − m ) = bε v + 1 M ( p − p ) . (A.11)Substituting Eq. (A.11) into Eq. (A.1), we obtain the fluid mass balance equation in terms ofthe pressure and the volumetric strain:1 M ∂p∂t + b ∂ε v ∂t + ∇ · v = f. (A.12)Linearizing the relation between volumetric total stress and volumetric strain with respect tothe reference state, σ v − σ v, = K dr ε v − b ( p − p ) , (A.13)allows us to express the change in porosity as the sum of a volumetric stress component anda fluid pressure component. From m = ρ f φ and Eq. (A.11), ρ f ρ f, φ − φ = bK dr ( σ v − σ v, ) + (cid:18) b K dr + 1 M (cid:19) ( p − p ) . (A.14)Using the effective stress equation, Eq. (A.13), we can rewrite Eq. (A.12) in terms of pressureand volumetric total stress: (cid:18) b K dr + 1 M (cid:19) ∂p∂t + bK dr ∂σ v ∂t + ∇ · v = f. (A.15) References
1. Aagaard, B., Kientz, S., Knepley, M., Strand, L., Williams, C.: PyLith User Manual,Version 1.8.0. Computational Infrastructure for Geodynamics, University of California,Davis (2012)2. Aagaard, B., Kientz, S., Knepley, M., Strand, L., Williams, C.: Pylith user manual: version2.1. 0. Davis, CA: Computational Infrastructure of Geodynamics (2013)3. Aagaard, B.T., Knepley, M.G., Williams, C.A.: A domain decomposition approach toimplementing fault slip in finite-element models of quasi-static and dynamic crustal de-formation. J. Geophys. Res. , 3059–3079 (2013)4. Aagaard, B.T., Knepley, M.G., Williams, C.A.: A domain decomposition approach toimplementing fault slip in finite-element models of quasi-static and dynamic crustal de-formation. Journal of Geophysical Research: Solid Earth (6), 3059–3079 (2013)5. Abousleiman, Y., Cheng, A.H.D., Cui, L., Detournay, E., Roegiers, J.C.: Mandel’s problemrevisited. G´eotechnique (2), 187–195 (1996)6. Ahmed, E., Nordbotten, J.M., Radu, F.A.: Adaptive asynchronous time-stepping, stop-ping criteria, and a posteriori error estimates for fixed-stress iterative schemes for cou-pled poromechanics problems. Journal of Computational and Applied Mathematics ,112312 (2020)7. Almani, T., Kumar, K., Dogru, A., Singh, G., Wheeler, M.F.: Convergence analysis of mul-tirate fixed-stress split iterative schemes for coupling flow with geomechanics. ComputerMethods in Applied Mechanics and Engineering , 180–207 (2016)4 Saumik Dana et al.8. Almani, T., Kumar, K., Wheeler, M.F.: Convergence and error analysis of fully discrete it-erative coupling schemes for coupling flow with geomechanics. Computational Geosciences(2017)9. Aziz, K., Settari, A.: Petroleum Reservoir Simulation. Elsevier, London (1979)10. Bear, J.: Dynamics of Fluids in Porous Media. Wiley, New York (1972)11. Biot, M.A.: General theory of three-dimensional consolidation. J. Appl. Phys. , 155–164(1941)12. Blacker, T.D., Bohnhoff, W.J., Edwards, T.L.: Cubit mesh generation environment. vol-ume 1: Users manual. Tech. rep., Sandia National Labs., Albuquerque, NM (United States)(1994)13. Borregales, M., Kumar, K., Radu, F.A., Rodrigo, C., Gaspar, F.J.: A partially parallel-in-time fixed-stress splitting method for biot’s consolidation model. Computers and Math-ematics with Applications (6), 1466 – 1478 (2019). 7th International Conference onAdvanced Computational Methods in Engineering (ACOMEN 2017)14. Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A.: Anderson accelerated fixed-stresssplitting schemes for consolidation of unsaturated porous media. Computers and Math-ematics with Applications (6), 1479 – 1502 (2019). 7th International Conference onAdvanced Computational Methods in Engineering (ACOMEN 2017)15. Cao, H.: Development of techniques for general purpose simulators. Ph.D. thesis, StanfordUniversity (2002)16. Castelletto, N., White, J.A., Ferronato, M.: Scalable algorithms for three-field mixed finiteelement coupled poromechanics. Journal of Computational Physics , 894–918 (2016)17. Castelletto, N., White, J.A., Tchelepi, H.A.: Accuracy and convergence properties of thefixed-stress iterative solution of two-way coupled poromechanics. International Journal forNumerical and Analytical Methods in Geomechanics (14), 1593–1618 (2015)18. Chang, H., Chen, Y., Zhang, D., et al.: Data assimilation of coupled fluid flow and geome-chanics using the ensemble kalman filter. SPE Journal (02), 382–394 (2010)19. Council, N.R., et al.: Induced seismicity potential in energy technologies. NationalAcademies Press (2013)20. Coussy, O.: Mechanics of Porous Continua. John Wiley and Sons, Chichester, England(1995)21. Coussy, O.: Poromechanics of freezing materials. J. Mech. Phys. Solids , 1689–1718(2005)22. Dana, S., Ganis, B., Wheeler, M.F.: A multiscale fixed stress split iterative scheme forcoupled flow and poromechanics in deep subsurface reservoirs. Journal of ComputationalPhysics , 1–22 (2018)23. Dana, S., Ita, J., Wheeler, M.F.: The correspondence between voigt and reuss bounds andthe decoupling constraint in a two-grid staggered algorithm for consolidation in heteroge-neous porous media. Multiscale Modeling & Simulation (1), 221–239 (2020)24. Dana, S., Wheeler, M.F.: Convergence analysis of fixed stress split iterative scheme foranisotropic poroelasticity with tensor biot parameter. Computational Geosciences (5),1219–1230 (2018)25. Dana, S., Wheeler, M.F.: Convergence analysis of two-grid fixed stress split iterativescheme for coupled flow and deformation in heterogeneous poroelastic media. ComputerMethods in Applied Mechanics and Engineering , 788–806 (2018)26. Galloway, D.L., Hoffmann, J.: The application of satellite differential sar interferometry-derived ground displacements in hydrogeology. Hydrogeology Journal (1), 133–154(2007)27. Galloway, D.L., Hudnut, K.W., Ingebritsen, S., Phillips, S.P., Peltzer, G., Rogez, F., Rosen,P.: Detection of aquifer system compaction and land subsidence using interferometric syn-thetic aperture radar, antelope valley, mojave desert, california. Water Resources Research (10), 2573–2585 (1998)28. Gaucher, E., Schoenball, M., Heidbach, O., Zang, A., Fokker, P.A., van Wees, J.D., Kohl,T.: Induced seismicity in geothermal reservoirs: A review of forecasting approaches. Re-newable and Sustainable Energy Reviews , 1473 – 1490 (2015)29. Geertsma, J.: The effect of fluid pressure decline on volumetric change of porous rocks.Trans. AIME , 331–340 (1957)30. GeoQuest, S.: Eclipse reference manual. Schlumberger, Houston, Texas (2014)31. Hesse, M.A., Stadler, G.: Joint inversion in coupled quasi-static poroelasticity. Journal ofGeophysical Research: Solid Earth (2), 1425–1445 (2014)itle Suppressed Due to Excessive Length 1532. Hughes, T.J.R.: The Finite Element Method: Linear Static and Dynamic Finite ElementAnalysis. Prentice-Hall, Englewood Cliffs, NJ (1987)33. Iglesias, M.A., McLaughlin, D.: Data inversion in coupled subsurface flow and geomechan-ics models. Inverse Problems (11), 115009 (2012)34. Jha, B., Bottazzi, F., Wojcik, R., Coccia, M., Bechor, N., McLaughlin, D., Herring, T.,Hager, B.H., Mantica, S., Juanes, R.: Reservoir characterization in an underground gasstorage field using joint inversion of flow and geodetic data. International Journal forNumerical and Analytical Methods in Geomechanics (14), 1619–1638 (2015)35. Jha, B., Juanes, R.: A locally conservative finite element framework for the simulation ofcoupled flow and reservoir geomechanics. Acta Geotechnica , 139–153 (2007)36. Jha, B., Juanes, R.: Coupled multiphase flow and poromechanics: A computational modelof pore pressure effects on fault slip and earthquake triggering. Water Resources Research (5), 3776–3808 (2014)37. Kim, J., Tchelepi, H., Juanes, R.: Stability, accuracy and efficiency of sequential methodsfor coupled flow and geomechanics. SPE Journal (2), 249–262 (2011)38. Kim, J., Tchelepi, H.A., Juanes, R.: Stability, accuracy and efficiency of sequential methodsfor coupled flow and geomechanics. Soc. Pet. Eng. J. (2), 249–262 (2011)39. Lackner, K.S.: A guide to co2 sequestration. Science (5626), 1677–1678 (2003)40. LaGriT: Los alamos grid toolbox,(lagrit) los alamos national laboratory (2013)41. Lie, K.A.: An introduction to reservoir simulation using MATLAB/GNU Octave: Userguide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge UniversityPress (2019)42. Lu, X., Wheeler, M.F.: Three-way coupling of multiphase flow and poromechanics inporous media. Journal of Computational Physics , 109053 (2020)43. Majer, E., Nelson, J., Robertson-Tait, A., Savy, J., Wong, I.: Protocol for addressing in-duced seismicity associated with enhanced geothermal systems. US Department of Energyp. 52 (2012)44. Mandel, J.: Consolidation des sols (´etude math´ematique)*. G´eotechnique (7), 287–299(1953)45. Mikeli´c, A., Wang, B., Wheeler, M.F.: Numerical convergence study of iterative couplingfor coupled flow and geomechanics. Computational Geosciences (3), 325–341 (2014)46. Mikeli´c, A., Wheeler, M.: Convergence of iterative coupling for coupled flow and geome-chanics. Computational Geosciences (3), 455–461 (2013)47. Orr, F.M.: Onshore geologic storage of co2. Science (5948), 1656–1658 (2009)48. Pan, H., Cao, H.: User Manual for General Purpose Research Simulator. Stanford Uni-versity Petroleum Engineering Institute, Stanford, CA (2010)49. Phillips, P.J., Wheeler, M.F.: A coupling of mixed and continuous Galerkin finite elementmethods for poroelasticity I: the continuous in time case. Comput. Geosci. , 131–144(2007)50. Phillips, P.J., Wheeler, M.F.: A coupling of mixed and continuous Galerkin finite elementmethods for poroelasticity II: the discrete-in-time case. Comput. Geosci. , 145–158(2007)51. Radu, F.A., Nordbotten, J.M., Pop, I.S., Kumar, K.: A robust linearization scheme for fi-nite volume based discretizations for simulation of two-phase flow in porous media. Journalof Computational and Applied Mathematics , 134 – 141 (2015)52. Rutqvist, J., Rinaldi, A.P., Cappa, F., Moridis, G.J.: Modeling of fault reactivation andinduced seismicity during hydraulic fracturing of shale-gas reservoirs. Journal of PetroleumScience and Engineering , 31–44 (2013)53. Storvik, E., Both, J.W., Kumar, K., Nordbotten, J.M., Radu, F.A.: On the optimization ofthe fixed-stress splitting for biot’s equations. International Journal for Numerical Methodsin Engineering p. nme.6130 (2019)54. Szulczewski, M.L., MacMinn, C.W., Herzog, H.J., Juanes, R.: Lifetime of carbon cap-ture and storage as a climate-change mitigation technology. Proceedings of the NationalAcademy of Sciences (14), 5185–5189 (2012)55. Vasco, D.: Estimation of flow properties using surface deformation and head data: Atrajectory-based approach. Water Resources Research (10) (2004)56. Vasco, D., Karasaki, K., Kishida, K.: A coupled inversion of pressure and surface displace-ment. Water resources research (12), 3071–3089 (2001)57. Vasco, D., Rucci, A., Ferretti, A., Novali, F., Bissell, R., Ringrose, P., Mathieson, A.,Wright, I.: Satellite-based measurements of surface deformation reveal fluid flow associatedwith the geological storage of carbon dioxide. Geophysical Research Letters (3) (2010)6 Saumik Dana et al.58. White, J.A., Castelletto, N., Klevtsov, S., Bui, Q.M., Osei-Kuffuor, D., Tchelepi, H.A.: Atwo-stage preconditioner for multiphase poromechanics in reservoir simulation. ComputerMethods in Applied Mechanics and Engineering , 112575 (2019)59. White, J.A., Castelletto, N., Tchelepi, H.A.: Block-partitioned solvers for coupled porome-chanics: A unified framework. Computer Methods in Applied Mechanics and Engineering , 55–74 (2016)60. Wilschut, F., Peters, E., Visser, K., Fokker, P.A., van Hooff, P., et al.: Joint history match-ing of well data and surface subsidence observations using the ensemble kalman filter: afield study. In: SPE Reservoir Simulation Symposium. Society of Petroleum Engineers(2011)61. Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E., Bruhn,D.: Analysis of induced seismicity in geothermal reservoirs – an overview. Geothermics52