We present a computational framework for fast monitoring of fault stabil- ity and ground deformation in multiphase geomechanics and demonstrate its efficacy for a carbon sequestration–enhanced oil recovery case study. The staggered solution algorithm for the coupled problem is augmented with a feature that allows for the flow and geomechanics sub-problems to be solved on different unstructured tetrahedral grids. For the field scale problem, the geomechanics grid goes all the way to the free surface while the flow grid is truncated at a depth above which the layers are impermeable. This frame- work avoids the unnecessary computational burden associated with equili- brating the initial pressure solution in the overburden, allows for a study of the critical interaction between overburden and faults, allows for fast rendi- tions of ground deformation, and allows a choice of resolution for the flow and geomechanics grids independently to capture disparate length scales of the underlying physics. Keywords: Carbon sequestration, Multiphase geomechanics, Fault stability, Ground deformation, Unstructured tetrahedral grids, Computational geometry 1. Introduction The US Department of Energy Regional Carbon Sequestration Partner- ship (RCSP) Initiative (see Fig. 1) was launched in year 2003 to advance development and testing of effective carbon capture, utilization and storage (CCUS) technologies across the United States. With the amount of planned CO2injection in these projects in the order of millions of metric tonnes, the quantification of geomechanical impact becomes critical in order to create guiding protocols for the design and monitoring of CCUS operations [3–14]. Figure 1: Locations of RCSP large-scale projects (Source: [1]) and geologic extent of Farnsworth (Source: [2]) Given the issues of poor resolution, interpretation uncertainty and partial coverage that are associated with geodetic (GPS and InSAR) measurements and recent push towards continuous monitoring of injection sites as part of the traffic light system [15], the need for predicting ground deformation to monitor subsidence or uplift has grown. Moreover, this capability is critical in cases where multiple forward simulation runs are performed as part of an inverse analysis to estimate the subsurface porosity/permeability fields from ground deformation data obtained from InSAR and/or GPS. Typically, the numerical grid in the region where the bulk of the flow-geomechanical coupling activity takes places carries a resolution that is required for the simulation accuracy. That region is typically the thin sliver that represents the reservoir. To simulate ground deformation, and any local geomechanical activity in the overburden would require a mesh that is extruded from the reservoir grid and extended all the way to the ground surface. This makes the geomechanical grid computationally expensive. As a result, it is optimal to truncate the flow grid to a certain depth above which the flow activity is not the subject of interest. But the problem is that in case of the single grid approach, since both flow and geomechanics are resolved on the same grid, even the geomechanical grid would need to be truncated at the same depth, thereby obviating a study of ground de- formation, and indeed any local geomechanical activity in the overburden. As we shall see in the field case, the overburden carries a layered structured with different rock lithologies, and any study that truncates the geomechan- ical grid at the reservoir depth would not allow encompassing geomechanical processes, e.g. caprock fracturing, above the reservoir for the concomitant study from a geomechanical standpoint. 2 coupled problem free surface g Solve geomechanics on geomechanics mesh Solve ow on ow mesh project project(a) Enabling simultaneous study of fault stability and ground deformation 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) Flowchart Figure 2: The two grid method In lieu of all of the above, we present a technique which allows for spa- tial decoupling of the flow and geomechanics domains with the subproblems being resolved on separate grids (Fig. 2a). The geomechanical equations are solved using the finite element method whereas the multiphase flow equa- tions are solved using the finite volume method. Typically, in such problems, the geomechanics mesh is expected to be coarser than, or the same, as the flow mesh everywhere the two meshes co-exist, but we generalize the method to the cases where the geomechanics elements can be smaller than the flow elements in localized regions where capturing the mechanics is more perti- nent, e.g. a basement fault below the reservoir, a bounding fault outside the reservoir, or fractures in the caprock, all of which may slip due to stresses transferred from the reservoir. Furthermore, we would like to generalize the notion of the two-grid from structured hexahedral meshes [16] to unstruc- tured tetrahedral meshes. As shown in Fig. 2b, the two-grid approach is built on top of the popular fixed stress split iterative scheme [16–52]. 2. Governing equations Under the quasi-static assumption for earth displacements, the govern- ing equation for linear momentum balance of the solid/fluid system can be expressed as ∇ ·σ+ρbg=0, (1) 3 where σis the Cauchy total stress tensor, gis the gravity vector, and ρb= ϕPnphase β ρβSβ+ (1−ϕ)ρsis the bulk density, ρβandSβare the density and saturation of fluid phase β,ρsis the density of the solid phase, ϕis the true (Eulerian) porosity, and nphaseis the number of fluid phases. Note that, by definition, the sum of all fluid phase saturations satisfiesPnphase β Sβ≡1. Assuming that the fluids are immiscible, the mass-conservation equation for each phase αis dmα dt+∇ ·wα=ραfα, (2) where the accumulation term dmα/dtdescribes the time variation of fluid mass relative 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 αthat can model the effect of wells. Eqs. (1) and (2) are coupled by virtue of poromechanics. Following [53], the poroelasticity equations can be written in incremental form as δσ=Cdr:δε−bδ˜p1≡Cdr:δε−X βbβpβ1, (3) where ˜ p=P βSβpβ, and bβ=bSβare the Biot coefficients for individual phases such thatP βbβ=b,Cdris the rank-4 drained elasticity tensor, 1is the rank-2 identity tensor, bis the Biot coefficient of the saturated porous material, εis the linearized strain tensor, defined as the symmetric gradient of the displacement vector u, ε:=1 2 ∇u+∇Tu ≡1 3εv1+e, (4) where εv= tr(ε) is the volumetric strain and eis the deviatoric strain tensor. In the deformed configuration, the mass of phase αper unit volume of porous medium is mα=ραSαϕJ (5) where J= 1 + εvis the Jacobian of deformation under infinitesimal defor- mation theory. In lieu of Eqs. (3)-(5), Eq. (2) can be rewritten as ∂ ∂t ραX β Nαβ+bαbβ Kdr pβ! +1 Kdr∂ ∂t(ραbασv) +∇ ·wα=ραfα(6) where Nαβis given in terms of saturation, porosity, fluid compressibility, Biot coefficients, drained modulus and the slope of the capillary pressure curve 4 The derivation of Eq. 6 is given in Appendix A. Because we assume that the fluids are immiscible, the mass-flux of phase αiswα=ραvα, where we adopt the traditional multiphase-flow extension of Darcy’s law [54, 55]: vα=−kkr α µα(∇pα−ραg), (7) where µαandkr αare the dynamic viscosity and the relative permeability of phase αin presence of other fluid phases. For single-phase water flow, α=β=w,Pαβ= 0, Sα= 1, cα=cw, and kr α= 1, as a result of which Eq. (6) reduces to ∂ ∂t ρw1 M+b2 Kdr pw +1 Kdr∂ ∂t(ρwbσv) +∇ ·(ρwvw) =ρwfw.(8) where1 M=ϕcw+(b−ϕ)(1−b) Kdris the inverse Biot modulus. 3. Two-grid algorithm and computational framework Algorithm 1: Two grid fixed-stress split algorithm for multiphase geomechanics 1Construct two-grid operators; 2Initialize pore pressures and equilibrate stresses at t=0; 3while t<Tdo 4 while Not converged do /* Loop over coupling iterations */ 5 while Not converged do /* Loop over Newton iterations */ 6 Solve multiphase flow keeping volumetric stress rate fixed; 7 Project pressures onto geomechanics mesh; 8 Solve geomechanics; 9 Project volumetric strains onto flow mesh; 10 t←t+dt; The fully discrete system of equations is presented in Appendix B. As shown in Algorithm 1, the basic idea of the two-grid approach is to transfer pore pressures from the flow mesh to the geomechanics mesh after the flow 5 Algorithm 2: Construction of flow-to-mechanics operators 1forLoop over all geomechanics elements do /* index igeo */ 2 forLoop over all flow elements do /* index iflow */ 3 count1 ←0;// Cases where flow element vertices are inside geomechanics element but not vice-versa 4 forLoop over four vertices of flow element do 5 ifVertex inside geomechanics element then 6 count1 ←count1 + 1; 7 ifcount1!=0 then 8 v1.push back( e−d); v2.push back(iflow); // d is distance between cell centers of element pair 9 count2 ←0;// Cases where geomechanics element vertices are inside flow element but not vice-versa 10 forLoop over four vertices of geomechanics element do 11 ifVertex inside flow element then 12 count2 ←count2 + 1; 13 ifcount2!=0 then 14 v1.push back( e−d); v2.push back(iflow); 15 flowtogeo.push back(v1/accumulate(v1.begin(), v1.end(), 0.0)); projection index.push back(v2); geo toflow.push back(v1); solve and conversely transfer volumetric strain from the geomechanics mesh to the flow mesh after the geomechanics solve in every coupling iteration in every time step. The transferred pore pressure at every geomechanics element is a weighted average of the pore pressures at the collection of flow elements interacting with the geomechanics element. Likewise, the transferred volu- metric strain at each flow element is a weighted average of the volumetric strains at the collection of geomechanics elements interacting with the flow element. The weights are the inverse exponential of the distances between the cell centers of the element pair. The reason for choosing the inverse expo- nential is to factor in the possibility of extremely low values for this distance for element pairs very close to each other. The collection of elements is ob- 6 Algorithm 3: Construction of mechanics to flow operators 1forLoop over all flow elements do 2 sum←0; 3 forLoop over all geomechanics elements do 4 forj=0→projection index[igeo].size()-1 do 5 ififlow==projection index[igeo][j] then 6 sum = sum + geo toflowvol[igeo][j]; 7 forLoop over all geomechanics elements do 8 forj=0→projection index[igeo].size()-1 do 9 ififlow==projection index[igeo][j] then 10 geotoflowvol[igeo][j] ←geotoflowvol[igeo][j]/sum; Algorithm 4: Usage of barycentric coordinates to determine if a point is inside or outside a tetrahedron 1forLoop over all four vertices of tetrahedron do 2 call baryth(x(1,:), x(2,:), x(3,:), x(4,:), p, s); // x(i,:) are coordinates of the ith vertex of tetrahedron, p is the point, and s(1:4) are barycentric coordinates of p 3 ifs(1)≥0And s(2) ≥0And s(3) ≥0And s(4) ≥0then 4 inside=.true.; // Point inside tetrahedron 5 else 6 outside=.true.; // Point outside tetrahedron tained by checking all flow elements which have atleast one vertex inside the geomechanics element as well as flow elements which have atleast one vertex of the geomechanics element inside it, as explained in Algorithms 2 and 3. We use barycentric coordinates to determine if a point is inside or outside a tetrahedron, as explained in Algorithm 4. The computation of pressures at a geomechanics element and volumetric strains at a flow elements are eluci- dated in Figs.3a and 3b respectively. We use STL vectors [56] to deploy the two grid projection operators, as shown in Algorithm 2. 7 ow element geo element ow element p1 p2 pgeo=p1exp(d1)+p2exp(d2) exp(d1)+exp( d2)+ ow element p3 d1 d2 d3(a) Flow to geomechanics operator ow element geo element 1 flow=1exp(d1)+2exp(d2) exp(d1)+exp( d2)+ d1 2 geo element d2(b) Geomechanics to flow operator Figure 3: dis a generic variable used to denote the distance between the cell centers of an element pair 4. Farnsworth Unit field simulation Figure 4: Faults within and near the Farnsworth Unit boundary (blue outline) along the Morrow layer. Source: [57] The study area for this work is the CO 2storage site called the Farnsworth Unit (FWU) oilfield, a deep structural basin approximately 70,000 square miles in size covering the northern Texas panhandle, western Oklahoma, eastern Colorado and western Kansas (see Fig. 1). As shown in Fig. 4, nine faults which intersected the Morrow-B sandstone reservoir and the Thirteen Finger limestone caprock were identified within FWU during 3-D seismic in- terpretation [58–62]. Because of their locations, Faults 1 through 3 have the greatest potential for impacting hydrocarbon production and carbon seques- 8 (a) Flow mesh, 31882 elements (b) Geomechanics mesh, 42108 elements Figure 5: Separate flow and geomechanics grids for two-grid method fault 1 fault 3 fault 6 x y 10000 m 10000 m (a) Top view of flow domain, no flow boundary con- ditions 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 flow domain Figure 6: Boundary conditions for the two grid simulation Figure 7: Injectors are shown with blue dots, and producers are shown with red dots tration at FWU [57]. A conservative estimate of CO 2storage capacity of the injection zone exceeds 10 million tonnes [58–63]. As shown in Fig. 5, we employ two different grids for Flowsim and PyLith 9 to use our two-grid method. 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 6◦angle. Similarly, Fault 1 is bending to west with an angle of 6◦. Both faults 1 and 3 cross through the three horizons. However, Fault 6 is younger, that is, more recent than Morrow-B Sandstone, and terminates just below the reservoir. We assume fault 6 to be vertical based on regional geology. Since Faults 1 and 3 are intersecting with each other, they will share the same grid elements around the intersection region. Thus, for these elements, it is difficult to determine which fault movement to follow during simulations of fault slip. We address this situation for both meshes by splitting Fault 3 into two smaller faults, Fault 3L and 3R, resulting in a gap between Faults 3L and 1 and between Faults 3R and 1 . It is important to note that while the flow and geomechanics meshes are different, the faults are at the same physical location in both meshes. Transmissibilities of the three faults are set to zero. Compared to the scenario with conductive faults, zero-transmissibility faults cause the reservoir pressure changes resulting from fluid injection and production to be larger. This induces larger changes in the fault stresses, which allows us to simulate a mechanically conservative scenario in the field with regard to the maximum allowable CO 2injection pressure and storage capacity The simulation spans approximately 57 years of production and injection from January 1956 to June 2013. It is important to keep in mind the field was discovered in 1955 and was unitized in 1963. Secondary recovery commenced in 1964 and CO 2injection began in December 2010. So, although the field is now primarily a carbon storage site, the production/injection history reflects a net fluid withdrawal since it was producing oil since early 1956 to late 2010. As shown in Fig. 7, there are 76 wells in the field: 45 producers and 31 injectors. We use a three phase-three component black-oil model [64] to represent the multiphase flow system during flow simulation. The black-oil flow model must be initialized with phase-dependent fluid properties: com- pressibility, viscosity, and the dissolved phase volume fraction as functions of the fluid phase pressure. We provide these properties for oil, water, and gas phases that are representative of FWU fluids at reservoir conditions. No-flow boundary condition is imposed on all boundaries of the flow domain. Two of the lateral sides of the geomechanics domain is subject to zero normal dis- placement boundary condition while the other two lateral sides are subject to traction boundary conditions. The top surface of geomechanics domain, 10 which is the ground surface, is traction free while the bottom surface is sub- ject to zero normal displacement boundary condition. The flow domain is truncated at z=−1300 m. The geomechanics domain goes all the way to the free surface z= 0. Areal extent of both the domains is 10000 m ×10000 m. Both the domains extend up to z=−2000 m. The simulations were per- formed using the Amazon Web Services instance c5a.8xlarge using 16 MPI processes. The model is populated with poromechanical properties obtained using two representative well logging data sets. Different reservoir properties are assigned along and across the layers to model an anisotropic response. The average bulk density is ρb= 2480 kg/m3, the average sonic compressional velocity is Vp= 3600 m/s, and the average sonic shear velocity is Vs= 1920 m/s. Dynamic values of the bulk modulus and Young’s modulus can be calculated as follows: K=ρb V2 p−4 3V2 s , E =ρbV2 s 3V2 p−4V2 s V2 p−V2 s The values are K= 19 .95 GPa and E= 23 .79 GPa. The Biot coefficient is assumed to be b= 0.6, which is representative of consolidated clastic formations in FWU. Reservoir permeability values are obtained by mapping the permeability field of an older simulation model [58–62] on to the coupled geomechanical model. The older model is constructed in Eclipse, is uncoupled to geomechanics, and has been used to track CO 2migration and oil recovery in FWU. Also, the uncoupled Eclipse model covers a smaller area and has smaller mesh elements compared to the coupled geomechanical model. As a result, the coupled model overlaps only partially with the uncoupled model. We average the property value of all the elements of the uncoupled model that are intersected by a coupled model element to obtain the property value of the coupled element. In the non-overlapping region, the geomechanical elements are assigned the field-averaged value of the respective property. We assume a uniform initial porosity of 0.3 in the flow model. This honors the historical production and injection volumes and avoids complexity due to heterogeneity in the porosity field. 4.1. DCFF time series In this study, the FWU faults do not fail, which agrees with observations in the field. Our focus is on quantifying the likelihood of failure or the 11 Fault 1 Fault 3R Upward DCFF trajectory Upward DCFF trajectoryFigure 8: DCFF time series for nodes along the intersection of faults 1 and 3, and for nodes on fault 3R along the line closest to fault 1 change in mechanical stability induced by production and injection. The likelihood of fault failure is assessed by evaluating the change in Coulomb Failure Function, DCFF, which is a function of the changes in pressure and stress around the fault DCFF = CFF( t)−CFF(0) ≡(|τ|+µfσ′ n)t−(|τ|+µfσ′ n)0= ∆|τ|+µf∆σ′ n ≡∆|τ|+µf∆σ′ nz}| { (∆σn+b∆pf) (9) 12 where ∆ |τ|is the magnitude of change in the shear traction vector, µfis the fault friction coefficient, ∆ σnis the change in total normal traction and ∆pfis the change in fault pressure. Here, we consider tension to be positive and compression to be negative, so a positive DCFF indicates a destabilizing effect on the fault. We assume a uniform static friction coefficient value of µf= 0.6 on all the faults. It is evident from Fig. 8 that the DCFF time series for nodes inside the reservoir closely follow the trends of the time series for average reservoir pres- sure. When the average reservoir pressure drops, the DCFF reduces whereas when the average reservoir pressure increases, the DCFF increases. This im- plies that the fault is stable along the reservoir when there is net fluid with- drawal and tends to be unstable when there is net fluid injection. This also follows the trend that is projected by Eq. (9) which states that the DCFF in- creases when the fault pressure increases and reduces when the fault pressure reduces. We also note that as the CO 2storage is proceeding, the DCFF at the reservoir nodes is continuing to follow a monotonically increasing trend. We also observe that the fault nodes at the overburden and underburden, the trend is that the DCFF increases when the average reservoir pressure decreases and decreases when the average reservoir pressure increases. This suggests that fault is stable along the overburden and underburden when there is net fluid injection and tends to be unstable along the overburden and underburden when there is net fluid withdrawal. The established hy- pothesis in [3] that fluid injection destabilizes a fault located in the region where the pore pressure has increased, whereas fluid withdrawal destabilizes a fault at the region generally outside the reservoir, lends credence to this observation. 4.2. Ground deformation As shown in Fig. 9, the free surface experiences subsidence, with more subsidence observed at regions where the producer wells are present. From the elasticity theory, the ground deformation pattern is expected to show a smoother version of the subsurface deformation pattern, which results from superposition of deformation fields of individual well activities. This demonstrates the power of our two-grid framework; it allows building mod- els that extend up to the ground surface and predict horizontal and vertical ground deformation, which can be compared against geodetic (GPS and In- SAR) measurements. Such comparisons are useful for calibration of coupled models, assimilation of multiphysics data, and estimation of uncertain rock 13 top surface slice Figure 9: Time series of z-displacement for two nodes on the ground surface. Point 1 (solid line) is located in a fault compartment occupied by many production wells; whereas Point 2 (dash line) is in a larger, hydraulically-connected region with fewer wells. This explains the difference in their displacement magnitudes and trends. properties or geologic structures. This is difficult to do with a single-grid framework where the grid does not extend to the ground for computational reasons. 5. Conclusions and outlook We presented a computational framework for fast estimation of fault sta- bility and ground deformation in multiphase geomechanics pertaining to a field-scale CO 2storage project in an oilfield. The staggered solution algo- rithm is augmented with a computational geometry-based feature (Section 3) that allows for the flow and geomechanics sub-problems to be solved on dif- ferent unstructured tetrahedral grids. Fault stability is estimated by tracking the time series of the change in Coulomb failure function DCFF (Eq. (9)) at critical points on the fault. We observed that the fault stability in the reservoir region is dictated by the pore pressure perturbations while the fault stability in the overburden region is dictated more by the local fault stress field. We also noted from the trends in Fig. 8, that as the CO 2storage is proceeding, the faults will continue to experience a destabilizing effect in the reservoir region. This points to the fact that at a point, the destabiliz- ing effect may increase enough to trigger seismic/aseismic fault slip activity depending on the background stress level and fault’s frictional properties. 14 Appendix A. Equation of mass conservation for multiphase flow in deformable porous media The increment in mean total stress is given by dσv=1 3tr(δσ) =1 3tr (Cdr:δε−bδ˜p1) =Kdrdεv−bδ˜p (A.1) where Kdris drained bulk modulus of porous medium. Eq. (A.1) can also be written as dεv=dσv Kdr+b d˜p Kdr(A.2) From the solid mass conservation equation, the incremental porosity update is given as [65, 66] dϕ=b−ϕ Kdr (dσv+d˜p) (A.3) In the deformed configuration, the mass of phase αper unit volume of porous medium is mα=ραSαϕJ where J= 1 + εvis the Jacobian of deformation under infinitesimal defor- mation theory. The incremental fluid mass can be written as dmα ρα=ϕJ∂Sα ∂PαβdPαβ+ϕJS αcαdpα+ϕSαdεv+JSαdϕ, (A.4) where cα=1 ραdρα dpαis the compressibility of the fluid phase α,Pαβ=pβ−pαis the capillary pressure and ∂Sα/∂P αβis the inverse capillary pressure deriva- tive which can be obtained from one of the Pαβvs.Sαmodels commonly used in flow simulations, e.g. the van Genuchten model. Substituting Eqs. (A.3) and (A.2) in Eq. (A.4) and assuming εv≪1 and noting that bα=bSα, we get dmα ρα=ϕ∂Sα ∂PαβdPαβ+ϕSαcαdpα+ ϕb+b−ϕSα Kdrd˜p+bαdσv Kdr 15 which can be written as dmα ρα=ϕ∂Sα ∂PαβdPαβ+ϕSαcαdpα+ ϕb+b−ϕSα KdrX βSβδpβ+bαdσv Kdr =ϕ∂Sα ∂PαβdPαβ+ϕSαcαdpα+NSαX βSβδpβ+b2Sα KdrX βSβδpβ+bαdσv Kdr where N= (b−ϕ)(1−b)/Kdr. Defining Nαβfrom the first three terms on the right hand side and substituting the above in Eq. (2), we obtain the multiphase flow equation for phase αin a poroelastic medium as: ∂ ∂t ραX β Nαβ+bαbβ Kdr pβ! +1 Kdr∂ ∂t(ραbασv) +∇ ·wα=ραfα Appendix B. Fully discrete system of equations Let Ωpand Ωfbe our domains of interest for the geomechanics and flow problems with closed boundaries ∂Ωpand∂Ωfrespectively. The geome- chanical equations are solved using the finite element method whereas the multiphase flow equations are solved using the finite volume method. Appendix B.1. Fully discrete multiphase flow problem For the two-phase water-gas flow problem, we get Ngg=−ϕ∂Sw ∂Pwg+ϕSgcg+S2 gN, Ngw=Nwg=ϕ∂Sw ∂Pwg+SgSwN, Nww=−ϕ∂Sw ∂Pwg+ϕSwcw+S2 wN,(B.1) Let Ωf=Snf i=1Ωiwhere nfis the number of flow elements. Note that in each element, Sw+Sg≡1. We integrate the fluid phase mass conservation equation, Eq. (6), for each phase over each element i. For the water and gas 16 phases, this yields ∂ ∂tZ Ωiρw Nww+Nwg+bbw Kdr pg− Nww+b2 w Kdr Pwg dΩ +1 Kdr∂ ∂tZ ΩiρwbwσvdΩ−Z ∂Ωiww·nidΓ =Z ΩiρwfwdΩ, (B.2) ∂ ∂tZ Ωiρg Nwg+Ngg+bbg Kdr pg− Nwg+bgbw Kdr Pwg dΩ +1 Kdr∂ ∂tZ ΩiρgbgσvdΩ−Z ∂Ωiwg·nidΓ =Z ΩiρgfgdΩ (B.3) where we used the capillary pressure relation, Pwg=pg−pw, to eliminate the water phase pressure, the Biot coefficient relation, bg+bw=b, and integration by parts for the mass-flux term to express it as a surface integral. niis the outward normal to the boundary of element i. We approximate both the pressure and the saturation fields 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 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 flux term as a sum of integral fluxes between element iand its adjacent elements j: Z ∂Ωiww·nidΓ =nface,iX j=1Z Γijww·nijdΓ =nface,iX j=1Ww,ij, where nface,iis the number of faces of element i(nface,i= 4 for a tetrahedral mesh such as the ones we used for the Mandel and Farnsworth models above andnface,i= 6 for a hexahedral mesh), and nijis the outward normal at the face Γ ij. The inter-element flux Ww,ij, is evaluated using a two-point flux approximation method which uses the flow potentials, Φ w,i=pw,i− ρw,igziand Φ w,j=pw,j−ρw,jgzj, to approximate the flux through interface ij between two neighboring elements calculated from upstream cells as follows: Ww,ij=( Aij/(ℓi/ki+ℓj/kj)ρw,ikr w,i/µw,i(Φw,i−Φw,j),if Φ w,i>Φw,j, Aij/(ℓi/ki+ℓj/kj)ρw,jkr w,j/µw,j(Φw,i−Φw,j),otherwise 17 where ziandzjare the centroid depths of the elements (with z-axis pointing downward), Aijis the area of the interface, ℓiandℓjare the distances be- tween the centroid of the interface and the centroid of the respective adjacent elements, and kiandkjare the element permeabilities. Using the notation δt(·)≡(·)n+1−(·)n ∆t, 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 i−nface,iX j=1Wn+1 g,ij− ρn+1 gfn+1 gVb i, Rn+1 w,i=δt Vbρw Nww+Nwg+bbw Kdr pg− Nww+b2 w Kdr Pwg i +δt ρwbw KdrσvVb i−nface,iX j=1Wn+1 w,ij− ρn+1 wfn+1 wVb i where subscript i∈[1, nf] refers to the value at element i. It is important to note that δtσvis the same across all Newton iterations within a fixed stress iteration in lieu of the fixed stress constraint. The global residual vectors are obtained by assembly of the element residual vectors: Rg= [Rg,1, . . . , R g,nf]T,Rw= [Rw,1, . . . , R w,nf]T. In vector form, the system of nphase×nfalgebraic equations of the flow problem reads: Rg Rwn+1 =0 0 , (B.4) 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 We solve the linear system of equations of the multiphase flow problem using Newton’s method. The correction vector is the solution of the linear system of equations:"∂Rg ∂Pg∂Rg ∂Sw ∂Rw ∂Pg∂Rw ∂Sw#(k)δPg δSw(k) =−Rg Rw(k) , (B.5) where the partial derivatives of the residuals are evaluated using Eq. (B.1). 18 Γf Γf−Γf+ τlf,dfu+ u−+ −nf,σnMesh processing fault nodes in the meshFigure B.10: A fault surface Γ f(shown as a line in a 2D domain) is processed during mesh processing to create three surfaces: the positive side surface Γ f+containing u+, the negative side surface Γ f−containing u−, and the middle slip surface containing the fault effective traction vector lfand the slip vector df. The fault surface need not be planar i.e. it can be curved with the normal vector nfvarying in space on the fault surface. Appendix B.2. Fully discrete geomechanics problem 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 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 displacement u+and the negative side displacement u−, respectively, of the adjacent rock hosting the fault. All three nodes are physically collocated on the initial grid, so the cohesive elements representing the fault are zero thickness elements at t= 0 before any fault slip or opening. Slip on the fault is the displacement of the positive side relative to the negative side, (u+−u−)−df=0on Γ f, (B.6) where dfis the fault slip vector. We impose the effective traction on the fault by introducing a Lagrange multiplier, lf, which is a force per unit area required to satisfy Eq. B.7. The Lagrange node stores the Lagrange multiplier lfand the fault slip df. The problem statement is to find ( u,lf) belonging to appropriate functional spaces satisfying the essential boundary conditions 19 (u=uon Γ u) such that Z Ωp∇sη: (σ′−b˜p1)dΩp+Z Γf+η·(lf−bpfn)dΓ −Z Γf−η·(lf−bpfn)dΓ−Z Ωpη·ρbgdΩp−Z Γση·tdΓ = 0 , (B.7) Z Γf+η·u+dΓ−Z Γf−η·u−dΓ−Z Γfη·dfdΓ = 0 , (B.8) for all test functions η(ndim×1 vector) belonging to the appropriate func- tional space satisfying η=0on Γ u.tis the prescribed traction on the traction boundary Γ σof the domain, e.g. t=0on the traction-free ground surface. Here, we used the multiphase effective stress equation (Eq. (3)), pf is 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 flow 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 bound- ary with its domain, Γ f, partitioned into non-overlapping fault elements, Γf=Snf,elem j=1Γf,j, where the subscript findicates variables associated with a fault. The displacement, the Lagrange multiplier, and the slip fields are approximated as follows: u≈uh=nnodeX b=1ηbUb,lf≈lh=nf,nodeX ¯b=1η¯bL¯b,df≈dh=nf,nodeX ¯b=1η¯bD¯b, where subscript hindicates the finite element approximation, nnodeis the number of displacement nodes (including the positive and negative side nodes at the fault elements, and nf,nodeis the number of Lagrange nodes. Ubis the ndim×1 vector of nodal displacements for node bin the domain coordinate system. L¯bandD¯bare the ndim×1 Lagrange multiplier and slip vectors at Lagrange node ¯bof the fault, and they are in the fault coordinate system de- fined by left-lateral, reverse, and opening directions. After substitution of the finite element approximations into the weak form of the problem (Eqs. (B.7)– (B.8)), we obtain the discrete equations in residual form at all displacement 20 nodes aand all Lagrange nodes ¯ a: 0=Z ΩBT a(σ′ h−b˜ph1)dΩ +Z Γf+ηT a(lh−bpf,hn)dΓ −Z Γf−ηT a(lh−bpf,hn)dΓ−Z ΩηT aρb,hgdΩ−Z ΓσηT atdΓ, (B.9) 0=Z Γf+ηT ¯auh+dΓ−Z Γf−ηT ¯auh−dΓ−Z ΓfηT ¯adhdΓ (B.10) where ηaandBaare nodal matrices for the shape function and its gradient respectively. Equation (B.9) results in ndim×nnodeequations, and Eq. (B.10) 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(σ′n+1 h−b˜pn+1 h1)dΩ +Z Γf+ηT a(ln+1 h−bpn+1 f,hn)dΓ −Z Γf−ηT a(ln+1 h−bpn+1 f,hn)dΓ−Z ΩηT aρn+1 b,hgdΩ−Z ΓσηT atdΓ =0 Rn+1 l,¯a=Z Γf+ηT ¯aun+1 h+dΓ−Z Γf−ηT ¯aun+1 h−dΓ−Z ΓfηT ¯adn+1 hdΓ =0 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 me- chanical problem reads:Ru Rln+1 =0 0 , (B.11) which needs to be solved for the displacements Un+1= [Un+1 1, . . . ,Un+1 nnode]T and the Lagrange multipliers Ln+1= [Ln+1 1, . . . ,Ln+1 nf,node]T, subject to the con- straint Eq. (B.6) 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) 21 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) , (B.12) 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Ω, CT a¯b=∂Ru,a ∂L¯b=Z Γf+ηT aη¯bdΓ−Z Γf−ηT aη¯bdΓ Matrix Cis the part of the Jacobian associated with the slip constraint, Eq. (B.10), 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 prop- erties and boundary conditions, the Jacobian matrix does not change with time, although the residuals may change due to coupling with the flow and due to fault slip. This is a saddle-point problem.