In an earlier work, we introduced the two-grid method for unstructured tetrahedra applied to staggered solution of coupled ﬂow and poromechanics and demonstrated its accuracy for the benchmark Mandel’s problem. In this work, we employ the two-grid method to study surface deformation and fault stability in a ﬁeld scale carbon storage project. Keywords: Carbon storage, Multiphase geomechanics, Fault stability, Surface deformation, Unstructured tetrahedral grids 1. Motivation Figure 1: Locations of RCSP large-scale projects. Source: [1]. The Southwest Regional Partnership on Carbon Sequestration (SWP) is one of the seven RCSPs. Geologic extent of the Farnsworth ﬁeld, Source: [2] Carbon capture and geologic storage is the separation and capture of CO 2from emissions of industrial processes, and the transport and permanent storage of the CO 2in deep underground formations. The U.S. Department of Energy’s (DOE) Regional Carbon Sequestration Partnership (RCSP) Initiative was launched in 2003 to advance development and testing of e ective CO 2storage technologies across the United States. This national initiative consists of seven partnerships as shown in Fig. 1. As of September 2016, almost 10 million metric tons of CO 2have been stored in geologic formations via RCSP large-scale ﬁeld projects. With the further huge amount of planned CO 2injection in these large scale projects, the risk of induced seismicity due to fault reactivation is too important to ignore. These seismic risks can be assessed in silico by developing a computationally inexpensive framework with the capability to model multiphase geomechanics and fault stability in a huge domain. These in silico simulations lend to o ering guiding protocols for the design of these operations [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1, 14]. The major bugbear to developing such a framework is the computationally intractable size of the geomechanical grid no matter which numerical method is used to resolve the coupled system of equations. Although coupled ﬂow and geomechanics commercial simulators such as Eclipse (Schlumberger) [15] and STARS (CMG) [16] have a long history in practical reservoir management, the calculations associated with storing large Preprint submitted to arXiv amounts of CO 2for climate mitigation are usually not connected in a uniform framework. The ﬂow simulations require detailed models honouring the ﬂow connectivity of the geological formation, while geomechanical models must account for the mechanical connectivity, including discrete discontinuities such as faults or fractures. The physical processes involved in the creation of the reservoir formations include sedimentation, which generates layers, erosion, which leads to pinch-out or removal of layers, and faulting, which leads to coupling o di erent types of layers. The result is a complex geometrical model where the physical properties in the model and the grid can not be deﬁned independently. The two critical features in problems involving faults are (1) Fault curvature and intersections must be honored because the fault stresses depend on the fault geometry by deﬁnition, and (2) Faults must be extended above and below the production or injection interval i.e. in the overburden and basement to allow modeling of induced seismicity in basement and mechanical integrity failure in caprock. For a compact representation of the physical parameters, reservoir models generally require the use of unstructured grids [17, 18, 19, 20] In lieu of all of the coupled problem free surface g Solve geomechanics on geomechanics mesh Solve ow on ow mesh project project (a) New time step New coupling iteration Solve ow on ow grid if not converged if converged keeping rate of volumetric stress xed Project pressures onto geomechanics grid Solve geomechanics on geomechanics grid Project volumetric strains onto ow grid Check for convergence(b) Figure 2: (a) The two-grid method enables the study of induced seismicity in a huge geomechanical domain. Projection operators are used to transfer information from one grid to the other. Source: [21]. (b) Two grid staggered solution algorithm New time step New coupling iteration Solve ow on ow grid if not converged if converged keeping rate of volumetric stress xed Project pressures onto geomechanics grid Solve geomechanics on geomechanics grid Project volumetric strains onto ow grid Check for convergence Construct jacobian in absence of faults Add zero thickness elements on fault Employ fault rupture model to modify jacobian in presence of faults Solve for displacements on all nodes and lagrange multipliers on fault nodes Figure 3: Procedural framework to incorporate faults in the geomechanical system of equations. A brief discussion of the concept of zero thickness elements is provided in Appendix B. above, we employ a technique which allows for spatial decoupling of the ﬂow and geomechanics domains with the subproblems being resolved on separate unstructured tetrahedral grids (Fig. 2a). As shown in Fig. 2b, the two-grid approach is built on top of a staggered solution algorithm in which the ﬂow and geomechanics subproblems are solved sequentially and iteratively using a ﬁxed stress split iterative scheme [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 2 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 21, 46, 47, 48, 49, 50, 51, 52]. The accuracy of the two-grid method has been demonstrated for the classical Mandel’s problem [53, 54] in our earlier work [21]. In the presence of faults though, there is an added set of operations that need to be performed for the solution of the geomechanical problem as shown in Fig. 3. 1.1. Computational framework Seismic and geological information Reservoir ow properties, multiphase uid properties Reservoir ow model Geomechanical model Displacements and stress Fluid pressures and saturations owsim geomechanics grid nite volume grid vtk .out ow grid hdf5 MRST Figure 4: Workﬂow. Properties like porosity, permeability, ﬂuid density, viscosity, and compressibility are obtained from an older uncoupled model and transferred to the geomechanical model using the Matlab Reservoir Simulation Toolbox [55]. Source: [21] The governing equations of multiphase poromechanics are laid out in Appendix A. The workﬂow of the in- house simulator is elucidated in Fig. 4. The ﬂow simulator ﬂowsim is a general purpose, object-oriented, ﬁeld-scale simulator for multiphase (oil, water, and gas) and multicomponent (e.g. methane, ethane, decane, CO 2, N2, and water) ﬂow through subsurface. It is capable of handling unstructured grids and complex production and injection scenarios in the ﬁeld, such as wells perforated at multiple depths and ﬂowing under variable rate and pressure controls. The geomechanics simulator PyLith is a open-source ﬁnite element code for the simulation of static and dynamic large- scale deformation problems [56, 57]. It uses an implicit formulation to solve quasi-static problems and an explicit formulation to solve dynamic rupture problems. It allows localized deformation along discrete features, such as faults, and is well integrated with meshing codes, such as CUBIT for both tetrahedral and hexahedral meshes [58]. The details of the space and time discretization and construction of two-grid projection operators are provided in Appendix C. The multiphase ﬂow model is solved using ﬁnite volume method whereas the geomechanical equations are solved using ﬁnite element method. 2. Numerical simulations The study area for this work is the CO 2storage site chosen by the SWP called the Farnsworth Field Unit (FWU). It is a deep structural basin approximately 70,000 square miles in size and primarily located in the northern Texas panhandle and western Oklahoma, and extending into eastern Colorado and western Kansas (see Fig. 1). The ﬁeld was discovered in 1955 and was unitized in 1963. Secondary recovery commenced in 1964 and CO 2injection began in December 2010. The ﬁeld had sequestered approximately 901 ;556 metric tonnes of CO 2by June, 2016. The CO 2 for this project is anthropogenically sourced from a fertilizer and an ethanol plant as shown in Fig. 5a. The SWP has a goal of permanently sequestering more than 1 million metric tonnes of CO 2in this site. A conservative estimate of 3 (a) (b) Figure 5: (a) (Left) The border of FWU and some of the well locations. (Right) CO 2is anthropogenically sourced from the Arkalon Ethanol Plant in Liberal, Kansas and the Agrium Fertilizer Plant in Borger, Texas. Source: [59, 60]. (b) Faults within or near the Farnsworth Unit (outline shown in blue). Source: [61] fault 1 fault 3 fault 6 x y 10000 m 10000 m (a) Top view of ﬂow domain, no ﬂow boundary conditions on all boundaries fault 1 fault 3 fault 6 x y(b) Top view of geomechanics domain y z z=0, traction free morrow shale morrow B top morrow B bottom fault 3 2000 m (c) Side view of geomechanics domain y z morrow shale morrow B top morrow B bottom fault 3 700 m z=1300 m(d) Side view of ﬂow domain Figure 6: Boundary conditions on ﬂow and geomechanics models CO 2storage capacity of the injection zone exceeds 10 million tonnes [59, 62, 63, 64, 65, 66]. As shown in Fig. 5b, nine faults which intersected the Morrow B sandstone reservoir and the limestone caprock were interpreted within the FWU [62, 63, 64, 65, 66]. Because of their locations, Faults 1-3 have the greatest potential for impacting production and carbon sequestration at the FWU [61]. As shown in Fig. 6, we consider faults 1, 3 and 6 in our calculations. Faults 1 and 3 are cutting through the model while fault 6 ends below the reservoir layer. Fault 3 is located near the center of the model and dips toward south with a 6oangle. No ﬂow boundary condition is imposed on all boundaries of the ﬂow domain. We employ the prescribed fault rupture model (see Appendix B) and prescribe zero slip which means all faults are locked. Our focus in this study is on quantifying the likelihood of failure or the change in mechanical stability of the faults. As shown in Fig. 7a, we employ di erent grids for ﬂowsim and PyLith to go with our two-grid method. The grids are generated using CUBIT mesh generation software [58]. Since the numerical method employed for the ﬂow model is ﬁnite volume instead of ﬁnite element, we employ some further processing of the ﬂow grid ﬁle to render it for the ﬁnite volume calculations. Reservoir permeability values are obtained by mapping the permeability ﬁeld of an older model constructed in 4 (a) (Above) Flow mesh, 31882 elements (Below) Pylith mesh, 42108 elements (b) Fault 3 is broken into fault 3L (western part) and fault 3R (eastern part). Injectors are shown with blue dots, and producers are shown with red dots. Source: [60] Figure 7: Areal extent of both the domains is 10000 m10000 m. The geomechanics domain goes all the way to the free surface z=0 while the ﬂow domain is truncated at z=1300 m. Both the domains extend upto z=2000 m Eclipse on to the coupled geomechanical model [60]. The older model is uncoupled to mechanics, and has been used to track CO 2migration and oil recovery in FWU [2, 67]. We use a three phase-three component black-oil model [68] to represent the multiphase ﬂow system during ﬂow simulation. The simulation spans approximately 60 years of production and injection from January 1956 to July 2016. As shown in Fig. 7b for the single grid case, there are 76 wells in the ﬁeld: 45 producers and 31 injectors. 2.1. Observations The entire simulation runs over roughly 24000 days, but we plot the results at multiple time steps along the way to observe the fault behavior and the e ect of the injection /production activity on uplift /subsidence at the free surface. We have plotted the results at the end of 2400 days, which represent the time stamp at around 10 percent of the total simulation time. The surface deformation plot shown in Fig. 8a is counter-intuitive. As the number of producers are more than the number of injectors, we would typically expect subsidence to be observed at the free surface. Instead, we observe uplift at the free surface, which just demonstrates the fact that back-of-the-envelope guesses about ground deformation would be erroneous in most other cases as well. Furthermore, we see three patches of more uplift across the free surface as shown by the red color in Fig. 8a. These observations would be hard to predict without the presence of the simulator that models the overburden and goes all the way to the free surface. Moreover, this capability is critical in cases where multiple forward runs are performed for a subsequent inverse analysis based estimation of the subsurface porosity /permeability ﬁelds from ground deformation data obtained from InSAR and /or GPS. We also note that the faults being ‘locked’ does not mean we impose zero normal displacements at the fault nodes. As shown in Fig. 8b, at the bottom of the domain, the faults 1 and 3 have thin streaks of negative z displacement compared to the rest of the botoom of the domain which is held back by the zero displacement boundary condition. The faults being ‘locked’ means that the slip vector at the fault ﬁnite elements nodes is assigned a zero magnitude. This means that the fault is not allowed to slide, as the purpose of this study is fault stability and not fault slip. We then plot the z-traction (z being the depth direction) for faults 1 and the western part of fault 3 also refered to as fault 3L. We see that the 5 (a) ground deformation (b) fault kinematics Figure 8: z-displacement plots after roughly 2400 days (a) fault 1 (b) fault 3L Figure 9: z-traction plots after roughly 2400 days stresses in fault 1 are predominantly compressive which could be a consequence of the pressure of the overlying rock, or it could be a consequence that there are many producers lines up along fault 1 (see Fig. 7b). We do not see as many producers or injectors near fault 3L, which implies the stresses observed on fault 3L could primarily be a consequence of how the overlying rock behaves to the stimulus provided by the wells across the domain. It is important to note though, that the fault 3L does intersect fault 1, and the inﬂuence of the fault 1 stress ﬁeld on fault 3L cannot be ignored as well. It can be hypothesized that the stress ﬁeld of the predominant fault 1 has an inﬂuence on the stress ﬁeld of fault 3L and vice versa. 3. Continuing work As we have mentioned, the plots represent the stress ﬁeld and deformation ﬁeld at roughly 10 percent of the total simulation time. As and when the further time steps generate more output, we shall plot more results and summarize the observations in a future piece. We shall also compare and contrast the results with the single grid results, albeit the single grid results are plotted on the grid which was truncated at a depth from the free surface. We also intend to have plots for the single grid case where the grid goes all the way to the free surface. 6 Appendix A. Governing equations of single phase and multiphase poromechanics Appendix A.1. Balance laws We use a classical continuum representation in which the ﬂuids and the solid skeleton are viewed as overlapping continua [69, 70]. The governing equations for coupled ﬂow and geomechanics are obtained from conservation of mass and balance of linear momentum. We assume that the deformations are small, that the geomaterial is isotropic, and that the conditions are isothermal. Let pbe our domain of interest for the geomechanics problem and let @ pbe its closed boundary. Under the quasistatic assumption for earth displacements, the governing equation for linear momentum balance of the solid /ﬂuid system can be expressed as r+bg=0; (A.1) whereis the Cauchy total stress tensor, gis the gravity vector, and b= Pnphase  S +(1 )s;is the bulk density,  andS are the density and saturation of ﬂuid phase ,sis the density of the solid phase, is the true porosity, andnphase is the number of ﬂuid phases. The true porosity is deﬁned as the ratio of the pore volume Vpto the bulk volume Vbin the current (deformed) conﬁguration. Let fbe our domain of interest for the ﬂow problem and let @ fbe its closed boundary. Assuming that the ﬂuids are immiscible, the mass-conservation equation for each phase is dm dt+rw = f ; (A.2) where the accumulation term dm =dtdescribes the time variation of ﬂuid mass relative to the motion of the solid skeleton, w is the mass-ﬂux of ﬂuid phase relative to the solid skeleton, and f is the volumetric source term for phase . Balance equations (A.1) and (A.2) are coupled by virtue of poromechanics. On one hand, changes in the pore ﬂuid pressure lead to changes in e ective stress, and induce deformation of the porous material—such as ground subsidence caused by groundwater withdrawal. On the other hand, deformation of the porous medium a ects ﬂuid mass content and ﬂuid pressure. Appendix A.2. Single-phase poromechanics For isothermal single-phase ﬂow of a slightly compressible ﬂuid in a poroelastic medium with no stress depen- dence of permeability, the single-phase ﬂuid mass conservation equation reduces to dm dt+rw=ff; (A.3) where mis the ﬂuid mass content (ﬂuid mass per unit bulk volume of porous medium), fis the ﬂuid density, w=fv is the ﬂuid mass ﬂux (ﬂuid mass ﬂow rate per unit area and time), and vis the seepage velocity relative to the deforming skeleton, given by Darcy’s law: v=k  rpfg ; (A.4) where kis the intrinsic permeability tensor, is the ﬂuid dynamic viscosity and pis the pore-ﬂuid pressure [69]. It is useful to deﬁne the ﬂuid content variation , :=m f;0; (A.5) wherem=mm0is the increment in ﬂuid mass content with respect to the initial reference state, and f;0is the reference ﬂuid density. The self-consistent theory of poroelastic behavior proposed by [71] links the changes in total stress and ﬂuid pressure with changes in strain and ﬂuid content. Following [72], the poroelasticity equations can be written in incremental form as =Cdr:"bp1 =b"v+1 Mp;(A.6) 7 where Cdris the rank-4 drained elasticity tensor, 1is the rank-2 identity tensor, "is the linearized strain tensor, deﬁned as the symmetric gradient of the displacement vector u, ":=1 2 ru+rTu ; (A.7) and"v=tr(") is the volumetric strain. Note that we use the convention that tensile stress is positive. It is useful to express the strain tensor as the sum of its volumetric and deviatoric components: "=1 3"v1+e; (A.8) from which it follows that the volumetric stress v=tr()=3 satisﬁes: v=Kdr"vbp: (A.9) Equation (A.6) implies that the e ective stress in single-phase poroelasticity, responsible for skeleton deformation, is deﬁned in incremental form as 0:=+bp1: (A.10) Biot’s theory of poroelasticity has two coupling coe cients: the Biot modulus Mand the Biot coe cient b. They are related to rock and ﬂuid properties as [72] 1 M= 0cf+b 0 Ks;b=1Kdr Ks; (A.11) where cf=1=Kfis the ﬂuid compressibility, Kfis the bulk modulus of the ﬂuid, Ksis the bulk modulus of the solid grain, and Kdris the drained bulk modulus of the porous medium. To set the stage for the numerical solution strategy of the coupled problem, it is useful to write the ﬂuid mass balance equation (A.3) (the pressure equation) in a way that explicitly recognizes the coupling with mechanical deformation. Equations (A.6) state that the increment in ﬂuid mass content has two components: increment due to expansion of the pore space and increment due to increase in the ﬂuid pressure. Assuming small elastic deformations and applying linearization from the reference state to the current state, we can write Eqs. (A.6) as 0=Cdr:"b(pp0)1; (A.12) 1 f;0(mm0)=b"v+1 M(pp0): (A.13) Substituting Eq. (A.13) into Eq. (A.3), we obtain the ﬂuid mass balance equation in terms of the pressure and the volumetric strain: 1 M@p @t+b@"v @t+rv=f: (A.14) Linearizing the relation between volumetric total stress and volumetric strain with respect to the reference state, vv;0=Kdr"vb(pp0); (A.15) allows us to express the change in porosity as the sum of a volumetric stress component and a ﬂuid pressure compo- nent. From m=f and Eq. (A.13), f f;0  0=b Kdrvv;0+ b2 Kdr+1 M! (pp0): (A.16) Using the e ective stress equation, Eq. (A.15), we can rewrite Eq. (A.14) in terms of pressure and volumetric total stress: b2 Kdr+1 M!@p @t+b Kdr@v @t+rv=f: (A.17) 8 Appendix A.3. Multiphase poromechanics In the multiphase or partially saturated ﬂuid system, it is not possible to linearize Eqs. (A.6) around a reference state because [72]: 1. Gases are very compressible, 2. Capillary pressure e ects are intrinsically nonlinear, and 3. Phase saturations vary between 0 and 1 and, therefore, a typical problem samples the entire range of nonlinearity. Therefore, following [72], we use the incremental formulation of poromechanics for multiphase systems, which does not assume physical linearization of total stress from the initial state to the current (deformed) state. Under this formulation, we split the total stress on the porous material into two parts: one that is responsible for deformation of the solid skeleton (the e ective stress), and another component that is responsible for changes in the ﬂuid pressures, =Cdr:"X b p 1; (A.18) where b are the Biot coe cients for individual phases such thatP b =b, where bis the Biot coe cient of the saturated porous material. It is common to further assume that b are proportional to the respective saturations S [73, 74, 75]. The e ective stress concept allows us to treat a multiphase porous medium as a mechanically equivalent single-phase continuum [76, 77]. The appropriate form of the e ective stress equation in a multiphase system is still an active area of research [78, 79, 77, 80, 81, 82]. Here we use the concept of equivalent pressure [79] in the e ective stress equation (Eq. (A.18)), pE=X S p U; (A.19) where U=P R p dS is the interfacial energy computed from the capillary pressure relations [82]. The equivalent pressure accounts for the interface energy in the free energy of the system, and leads to a thermodynamically consistent and mathematically well-posed description of the multiphase ﬂuid response to the solid deformation [82]. For a system with two phases, the wetting phase wand the non-wetting phase g, the capillary pressure is Pc(Sw)Pwg(Sw)=pgpw; (A.20) and the interfacial energy is U=R1 SwPwgdS. Assuming b =bS [73, 74, 75], and using Eq. (A.19) in Eq. (A.18), we obtain the stress-strain relationship for multiphase linear poroelasticity: =0bpE1; 0=Cdr:": (A.21) In the deformed conﬁguration, the mass of phase per unit volume of porous medium is m = S (1+"v): (A.22) Note that, by deﬁnition, the sum of all ﬂuid phase saturations satisﬁesPnphase S 1. Extending Eq. (A.13) for multiphase systems [72, 83], we have dm ! =b d"v+X N dp ; (A.23) where N=M1is the inverse Biot modulus. In a multiphase system, the Biot modulus is a symmetric positive deﬁnite tensor M=[M ], and the Biot coe cient is a vector. To determine the coupling coe cients N as a function of the primary variables (pressure, saturations, and displacement), and rock and ﬂuid properties, we develop an alternate expression for the di erential increment in ﬂuid mass. Using Eq. (A.22), dm =d( S (1+"v)); (A.24) which can be expanded as dm ! = @S @P dP + S c dp + S d"v+S d ; (A.25) 9 where c =1  d dp is the compressibility of the ﬂuid phase , and@S =@P is the inverse capillary pressure derivative. Above, repeated indices do not imply summation and we have assumed inﬁnitesimal deformations. We can express the increment in porosity, d , as a function of the volumetric e ective stress, d0 v, to obtain a closed-form expression of Eq. (A.25). Let Vs=VbVpbe the volume of the solid matrix, and d"sv=dVs=Vs=dsv=Ksbe the volumetric dilation of the solid matrix, where svis the volumetric matrix stress. From d =d(Vp=Vb)= (dVp=VpdVb=Vb), we can write the incremental form of strain partition as (1 )d"v=(1 )d"sv+d : (A.26) Similarly, the volumetric Cauchy total stress can be partitioned into the volumetric matrix stress and the ﬂuid pressure as dv=(1 )dsv dpE: (A.27) Substituting dsvfrom Eq. (A.27) into Eq. (A.26), we obtain d =b Kdrd0 v+(1b)dpE: (A.28) Equation (A.28) implies that an increment in porosity is related to increments in volumetric e ective stress and ﬂuid pressures, similar to Eq. (A.16) in the single-phase case. Substituting d"vfrom Eq. (A.21) and d from Eq. (A.28) into Eq. (A.25) allows us to express the increment in the phase mass as a function of the increments in the total volumetric stress and phase pressures. Equating this to Eq. (A.23) yields the desired expressions for the coupling coe cients N , which for a water-gas system are Ngg= @Sw @Pwg+ Sgcg+S2 gN; Ngw=Nwg= @Sw @Pwg+SgSwN; Nww= @Sw @Pwg+ Swcw+S2 wN;(A.29) where N=(b )(1b)=Kdr, and the subscripts wandgdenote water and gas phases, respectively. Finally, we obtain the multiphase ﬂow equation for phase in a poroelastic medium by substituting the two constitutive relations, the e ective stress equation, Eq. (A.21), and the ﬂuid mass increment equation, Eq. (A.23), in the mass balance equation, Eq. (A.2): @ @t0BBBBBB@ X N +b b Kdr! p 1CCCCCCA+1 Kdr@ @t( b v)+rw = f ; 8 =1;:::; nphase: (A.30) Because we assume that the ﬂuids are immiscible, the mass-ﬂux of phase isw = v , where we adopt the traditional multiphase-ﬂow extension of Darcy’s law [84, 69]: v =kkr  (rp  g); (A.31) where andkr are the dynamic viscosity and the relative permeability of phase in presence of other ﬂuid phases. Appendix B. Fault mechanics We treat faults as surfaces of discontinuity embedded in the continuum across which displacement is allowed to be discontinuous to recognize the possibility of fault slip [22, 57]. As shown in Fig. B.10, each node on the fault is triplicated to create a “ +” node, a “–” node, and a Lagrange node in the middle. The side nodes store the positive side 10 fault nodes   +  +Figure B.10: The two sides of the fault surface, which need not be planar, are designated as the ‘ +’ side and the ‘’ side. Each fault node is triplicated to create a “ +” node, a “–” node, and a Lagrange node in the middle displacement u+and the negative side displacement u, respectively. All three nodes are physically collocated on the initial grid, so elements representing the fault are zero thickness elements. Slip on the fault is the displacement of the positive side relative to the negative side, (u+u)d=0on; (B.1) where dis the fault slip vector. We impose the e ective traction on the fault by introducing a Lagrange multiplier, l, which is a force per unit area required to satisfy equation B.1. The Lagrange node stores the Lagrange multiplier land the fault slip d. The reader is refered to [85, 86, 87] for an understanding of Lagrange multipliers with similar intent in a di erent context. The Jacobian of the geomechanical system after discretization using ﬁnite elements takes the form (see [57, 22] for a more detailed discussion) 266666666666664Knn Knn+Knn 0 Kn+nKn+n+ 0 LT p Knn 0 K nnLT p 0 L pLp 0377777777777775 where ndenotes degrees of freedom (DOF) not associated with the fault, ndenotes DOF associated with the negative side of the fault, n+denotes DOF associated with the positive side of the fault, and pdenotes DOF associated with the Lagrange multipliers. The vector being solved for is 8>>>>>><>>>>>>:un un+ un lp9>>>>>>=>>>>>>; Appendix B.1. Prescribed fault rupture In a prescribed (kinematic) fault rupture, we specify the slip-time history d(x;y;z;t) at every location on the fault surfaces. The vector being solved for is 8>>><>>>:un un lp9>>>=>>>; as the displacements un+on the “ +” nodes can be expressed in terms of displacements unon the “-” nodes using the fault slip. As a result, the number of DOF is considerably reduced. The slip-time history enters into the calculation of the residual as do the Lagrange multipliers, which are available from the current trial solution. In prescribing the slip-time history, we do not specify the initial tractions on the fault surface so the Lagrange multipliers are the change in the tractions on the fault surfaces corresponding to the slip. 11 Appendix B.2. Spontaneous fault rupture In contrast to prescribed fault rupture, in spontaneous fault rupture a constitutive model controls the tractions on the fault surface. The shear traction on the fault is computed as =jl(ln)nj where nis the fault normal vector which points from the negative side to the positive side of the fault. The frictional stress fon the fault as f=8>><>>: cfln;ln<0; c; ln0;(B.2) where cis the cohesive strength of the fault and fis the coe cient of friction. The coe cient of friction fis a function of the slip magnitude jdj, and drops linearly from its static value, s, to its dynamic value, d, over a critical slip distance dc, f=8>><>>:s(sd)jdj dc;jdjdc; d;jdj>dc:(B.3) We use the Mohr-Coulomb theory to deﬁne the stability criterion for the fault [88]. When the shear traction on the fault is below the friction stress,  f, the fault does not slip. When the shear traction is larger than the friction stress, > f, the contact problem is solved to determine the Lagrange multipliers and slip on the fault, such that the Lagrange multipliers are compatible with the frictional stress [22]. Appendix C. Fully discrete set of equations Appendix C.1. Fully discrete multiphase ﬂow problem The ﬂow problem is discretized using the ﬁnite volume method on a ﬁner grid than the geomechanics problem, f=Snf i=1 iwhere nfis the number of ﬂow elements. For simplicity, let us consider two ﬂuid phases, water and gas. Note that in each element, Sw+Sg1. We integrate the ﬂuid phase mass conservation equation, Eq. (A.30), for each phase over each element i. For the water phase, this yields @ @tZ iw Nww+Nwg+bbw Kdr! pg Nww+b2 w Kdr! Pwg! d +1 Kdr@ @tZ iwbwvd Z @ iwwnid =Z iwfwd ; (C.1) where we used the capillary pressure relation, Pwg=pgpw, to eliminate the water phase pressure, the Biot coe cient relation, bg+bw=b, and integration by parts for the mass-ﬂux term to express it as a surface integral. niis the outward normal to the boundary of element i. Similarly, we have a mass balance equation for the gas phase: @ @tZ ig Nwg+Ngg+bbg Kdr! pg Nwg+bgbw Kdr! Pwg! d +1 Kdr@ @tZ igbgvd Z @ iwgnid =Z igfgd : (C.2) We approximate both the pressure and the saturation ﬁelds with a piecewise constant interpolation function, ', such that'itakes a constant value of 1 over element iand 0 at all other elements. Phase pressures and saturations are approximated as p p ;h=nfX i=1'ip ;i;S S ;h=nfX i=1'iS ;i (C.3) 12 where the discrete pressures, p ;i, and phase saturations, S ;i, are located at the center of element i. We can further express the mass ﬂux term as a sum of integral ﬂuxes between element iand its adjacent elements j: Z @ iwwnid =nface;iX j=1Z i jwwni jd =nface;iX j=1Ww;i j; (C.4) where nface;iis the number of faces of element i, and ni jis the outward normal at the face i j. The inter-element ﬂux of water, Ww;i j, can be evaluated from Darcy’s law (Eq. (A.31)) as a function of the rock and ﬂuid properties, pressures, and saturations of element iand its adjacent elements j, using either a two-point or a multipoint ﬂux approximation [89, 90, 91]. The two-point ﬂux approximation method uses the ﬂow potentials, w;i=pw;iw;igzi andw;j=pw;jw;jgzj, to approximate the ﬂux through interface i jbetween two neighboring elements as follows: Ww;i j=Ti jw;i j w;iw;j ; (C.5) where ziandzjare the centroid depths of the elements (with z-axis pointing downward), w;i jis the water phase mobility and Ti jis the geometric transmissibility of the interface i j. The phase mobility is calculated from upstream cells, w;i j=8>><>>:w;ikr w;i=w;i;ifw;i>w;j; w;jkr w;j=w;j;otherwise:(C.6) The geometric transmissibility is estimated from the harmonic average of the respective element permeabilities and element sizes, Ti j=Ai j=(`i=ki+`j=kj), where Ai jis the area of the interface, `iand`jare the distances between the centroid of the interface and the centroid of the respective adjacent elements, and kiandkjare the element permeabilities, assumed isotropic. After substitution, the semi-discrete water phase mass balance equation is 0=@ @t" Vb;iw;i Nww+Nwg+bbw Kdr! ipg;i Nww+b2 w Kdr! iPwg;i!# +@ @t w;i bwv Kdr! iVb;i! nface;iX j=1Ww;i jw;ifw;iVb;i;8i=1;:::; nf; (C.7) where subscript irefers to the value at element i. Similarly, we could write the semi-discrete equation for the gas phase: 0=@ @t" Vb;ig;i Ngw+Ngg+bbg Kdr! ipg;i Ngw+bgbw Kdr! iPwg;i!# +@ @t g;i bgv Kdr! iVb;i! nface;iX j=1Wg;i jg;ifg;iVb;i;8i=1;:::; nf; (C.8) The fully discrete system of equations is Rn+1 g;i=t" Vbg Ngw+Ngg+bbg Kdr! pg Ngw+bgbw Kdr! Pwg!# i +t" gbg KdrvVb# inface;iX j=1Wn+1 g;i jh n+1 gfn+1 gVbi ii=1;:::; nf; (C.9) Rn+1 w;i=t" Vbw Nww+Nwg+bbw Kdr! pg Nww+b2 w Kdr! Pwg!# i +t" wbw KdrvVb# inface;iX j=1Wn+1 w;i jh n+1 wfn+1 wVbi ii=1;:::; nf: (C.10) 13 The global residual vectors are obtained by assembly of the element residual vectors: Rg=[Rg;1;:::; Rg;nf]T,Rw= [Rw;1;:::; Rw;nf]T. In vector form, the system of nphasenfalgebraic equations of the ﬂow problem reads: "Rg Rw#n+1 ="0 0# ; (C.11) which must be solved at every time step for the vector of element gas-phase pressures and water-phase saturations Pn+1 g=[pn+1 g;1;:::; pn+1 g;nf]T Sn+1 w=[Sn+1 w;1;:::; Sn+1 w;nf]T The water phase pressure is determined with the help of the capillary pressure relation (Eq. (A.20)), and the gas saturation is determined from the constraint Sg+Sw1. We solve the linear system of equations of the multiphase ﬂow problem using Newton’s method. The correction vector is the solution of the linear system of equations: 26666664@Rg @Pg@Rg @Sw @Rw @Pg@Rw @Sw37777775(k)"Pg Sw#(k) ="Rg Rw#(k) ; (C.12) where the partial derivatives of the residuals are evaluated using the constitutive equations (Eq. (A.29) and Eq. (A.31)) that relate rock and ﬂuid properties to the ﬂuid pressures and saturations. Flow problem is solved using the SAMG multigrid pre- conditioner [92, 93]. Appendix C.2. Fully discrete geomechanics problem Find ( u;l) belonging to appropriate functional spaces satisfying the essential boundary conditions ( u=uonu) such that Z prs: (0bpE1)d p+Z f+(lbpfn)dZ f(lbpfn)dZ pbgd pZ td =0; (C.13)Z f+u+dZ fudZ fdd =0; (C.14) for all test functions (ndim1 vector) belonging to the appropriate functional space satisfying =0onu. Here, we used the multiphase e ective stress equation (Eq. (A.21)), pEis the equivalent pressure (Eq. (A.19)), pfis the fault pressure, and bis the bulk density, all of which depend on the phase pressures and saturations, and, therefore, on the solution of the ﬂow problem. Let the domain be partitioned into non-overlapping elements (grid blocks), p=Snelem j=1 jwhere nelemis the number of geomechanics elements. For the remainder of this section, we shall use the variable instead of pfor the sake of notational convenience. A fault is treated as an interior boundary with its domain, f, partitioned into non-overlapping fault elements, f=Snf;elem j=1f;j, where the subscript findicates variables associated with a fault. The displacement, the Lagrange multiplier, and the slip ﬁelds are approximated as follows: uuh=nnodeX b=1bUb;llh=nf;nodeX ¯b=1¯bL¯b;ddh=nf;nodeX ¯b=1¯bD¯b; (C.15) where subscript hindicates the ﬁnite element approximation, nnodeis the number of displacement nodes (including the duplicated nodes at the fault elements, and nf;nodeis the number of Lagrange nodes. Ubis the ndim1 vector of nodal displacements for node bin the domain coordinate system, Ub=(Ub;x;Ub;y;Ub;z).L¯bandD¯bare the ndim1 Lagrange multiplier and slip vectors at Lagrange node ¯bof the fault, and they are in the fault coordinate system (left-lateral, reverse, and opening, respectively: L¯b=(L¯b;1;L¯b;2;L¯b;3) and D¯b=(D¯b;1;D¯b;2;D¯b;3). The interpolation functions, band¯b, are the usual C0-continuous isoparametric functions, such that they take a value of 1 at the respective nodes, and 0 at all other nodes [94]. In a Galerkin method, the test functions are chosen to be equal to the interpolation functions. 14 After substitution of the ﬁnite element approximations into the weak form of the problem (Eqs. (C.13)–(C.14)), we obtain the discrete equations in residual form at all displacement nodes aand all Lagrange nodes ¯ a: 0=Z BT a(0 hbpE;h1)d +Z f+T a(lhbpf;hn)dZ fT a(lhbpf;hn)d Z T ab;hgd Z T atd8a=1;:::; nnode; (C.16) 0=Z f+T ¯auh+dZ fT ¯auhdZ fT ¯adhd8¯a=1;:::; nf;node; (C.17) where we have deﬁned the nodal matrices for the shape function a=2666666664a0 0 0a0 0 0a3777777775: (C.18) and its symmetric gradient Ba=2666666666666666666666664@xa 0 0 0@ya 0 0 0 @za @ya@xa 0 0@za@ya @za 0@xa3777777777777777777777775; (C.19) corresponding to the compact engineering notation for stress and strain inside an element, 0 h=[0 h;xx;0 h;yy;0 h;zz;0 h;xy;0 h;yz;0 h;xz]T; "h=["h;xx;"h;yy;"h;zz;2"h;xy;2"h;yz;2"h;xz]T respectively [94]. The identity tensor in compact notation is 1=[1;1;1;0;0;0]T. Equation (C.16) results in ndimnnode equations, and Eq. (C.17) results in ndimnf;nodeequations. Since we use the prescribed fault rupture model, D¯ais prescribed at all fault nodes. The fully discrete system of equations is Rn+1 u;a=Z BT a(0n+1 hbpn+1 E;h1)d +Z f+T a(ln+1 hbpn+1 f;hn)d Z fT a(ln+1 hbpn+1 f;hn)dZ T an+1 b;hgd Z T atd =0 8a=1;:::; nnode; (C.20) Rn+1 l;¯a=Z f+T ¯aun+1 h+dZ fT ¯aun+1 hdZ fT ¯adn+1 hd =08¯a=1;:::; nf;node; (C.21) The global residual vectors are obtained by assembly of the nodal residual vectors: Ru=[Ru;1;:::; Ru;nnode]T,Rl= [Rl;1;:::; Rl;nf,node]T. In vector form, the system of ( ndimnnode+ndimnf;node) algebraic equations of the mechanical problem reads:"Ru Rl#n+1 ="0 0# ; (C.22) which needs to be solved for the displacements Un+1=[Un+1 1;:::; Un+1 nnode]Tat the displacement nodes, and the Lagrange multipliers Ln+1=[Ln+1 1;:::; Ln+1 nf,node]Tat the Lagrange nodes, subject to the constraint (Eq. (B.1)) at every time step. We solve the linear systems of equations of the mechanical problem using Newton’s method. Given an approximation [Un+1;Ln+1](k)to the solution at tn+1, an improved solution is obtained as [Un+1;Ln+1](k+1)=[Un+1;Ln+1](k)+[Un+1;Ln+1](k) 15 where the correction vector is the solution to the system of linear equations (removing the explicit reference to the time level n+1):" K CT C0#(k)"U L#(k) ="Ru Rl#(k) ; (C.23) where the block matrices are obtained via element-by-element assembly of the individual nodal contributions to the Jacobian: Kab=@Ru;a @Ub=Z BT aDBbd ; (C.24) CT a¯b=@Ru;a @L¯b=Z f+T a¯bdZ fT a¯bd: (C.25) Matrix Cis the part of the Jacobian associated with the slip constraint, Eq. (C.17), and consists of direction cosine matrices to convert from the global coordinate system to the fault coordinate system (vice versa for CT). Note that for a linear elastic material with time-independent material properties and boundary conditions, the Jacobian matrix does not change with time, although the residuals may change due to coupling with the ﬂow and due to fault slip. This is a saddle-point problem. Custom preconditioners are required to solve the linear system e ciently using the Portable, Extensible Toolkit for Scientiﬁc Computation (PETSc) [95] multigrid preconditioner for the elasticity submatrix in conjunction with a custom fault preconditioner for the Lagrange multiplier submatrix [56, 57]. Appendix C.3. Two-grid operators Ep Ef Ef Ep Figure C.11: Depiction of couple of intersecting ﬂow-geomechanics element pairs. We refer to the ﬂow element as Efand geomechanics element asEp. Depiction is to reiterate the point that there is no restriction of whether the ﬂow or geomechanics element needs to be smaller than the other The basic idea of the two-grid approach is to transfer pore pressures from the ﬂow domain to the geomechanics do- main and conversely transfer volumetric strain from the geomechanical domain to the ﬂow domain. The transfered pore pressure at every geomechanics element is a volume average of the pore pressures at the ﬂow elements inter- acting with the geomechanics element. Likewise, the transfered volumetric strain at each ﬂow element is a volume average of the volumetric strains at the geomechanics elements interacting with the ﬂow 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 ﬁnite elements depicted in Fig. C.11. The relationship between the ﬂow-geomechanics element pair is obtained by determining the barycentric coordinates of the vertices of the ﬂow element with respect to the geomechanics element and vice versa as explained in Algorithm 1. 16 Algorithm 1 Determining the relationship between ﬂow-geomechanics element pairs forLoop over all geomechanics elements do forLoop over all ﬂow elements do count1 0 forLoop over four vertices of ﬂow element do ifVertex inside geomechanics element then .Invoke barycentric coordinates count1 count1 +1 ifcount1! =0then Flow to geomechanics projection operator Meas (Ef) Geomechanics to ﬂow projection operator Meas (Ep) count2 0 forLoop over four vertices of geomechanics element do ifVertex inside ﬂow element then .Invoke barycentric coordinates count2 count2 +1 ifcount2! =0then Flow to geomechanics projection operator Meas (Ef) Geomechanics to ﬂow projection operator Meas (Ep)