University of Southern California Dissertations and Theses

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

(USC Thesis Other)

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content A COUPLED MULTIPHASE FLOW-GEOMECHANICS FRAMEWORK FOR MODELING MATERIAL, KINEMATIC AND CONTACT NONLINEARITIES IN FAULTED RESERVOIRS by Xiaoxi Zhao A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) December 2020 Copyright 2020 Xiaoxi Zhao Acknowledgments I would like to offer my deepest thanks to my advisor Prof. Birendra Jha, who introduced me to the field of Computational Geoscience. It has been an honor to workwithhimduringmyPh.D.study. Heisalwayssupportiveandalwaysprovides mewithvaluableandcriticalsuggestions. Iappreciateallhiscontributionsoftime, ideas, and supports to make my Ph.D. experience fruitful and insightful. I would like to express my sincere appreciation to Prof. Iraj Ershaghi who motivated me for starting this journey. The enthusiasm he has for study and research always inspires me in all the up and down times in the Ph.D. pursuit. He is always there to listened to my challenges and gives me continuous support, advice and encouragement. I thank all my committee members: Prof. Birendra Jha, Prof. Iraj Ershaghi and Prof. Muhammad Sahimi for their time and helpful comments. I would also like to thank the other three members of my qualifying exam committee, Prof. Kristian Jessen, Prof. Jincai Chang and Prof. Felipe De Barros, for their time and insightful questions. I would like to express my sincere gratitude to Dr. Fred Aminzadeh for his guidance, patience, and support during the first few years of my studies. I thank the members of the Jha Research Group for their invaluable friendship. ii Last but not least, I would like to thank my family for all their love and encouragement. For my parents, Baosheng Zhao and Changli He, who provide me endless love and support and always believe in me. For my husband Tongxin Zhang who are always there with me on this journey. This dissertation could not be possible without their unconditional support and love. iii Abstract Mitigating the hazard of fault reactivation during water, gas injection and oil production is important for ensuring environmental sustainability of geologic car- bon sequestration and hydrocarbon recovery in the long term. In this Dissertation, we present a novel multiphase flow-geomechanics framework that can capture the nonlinearities of material (poroplasticity), kinematic (finite strain) and contact (dynamic fault slip). Using our high-resolution framework, we compare the fault destabilization effects of gas injection with those of oil production. Focusing on the confined mul- tiphase reservoir located in the hanging wall block of a normal fault, we analyze the fault tractions, slip, and reservoir pressure to understand how a fault is desta- bilized in successive steps, how hypocenter and aftershock locations are selected, and what role the well type parameter (injection or production) can play in con- trolling the fault stability. We extract and summarize the characteristic features of induced fault activation in terms of the dependence of Coulomb Failure Stress (CFS) and the Fault Stress Ratio (FSR) on the dimensionless reservoir pressure perturbation caused by injection or production. By combining the finite strain poroplasticity with Mohr-Coulomb theory of fault failure, we conduct a fully-coupled multiphase flow and geomechanical study of production-induced fault slip in a poroplastic reservoir. We evaluate the impact iv of plasticity on reservoir pressure, deformation, induced stress, and fault slip by comparing infinitesimal strain elastic and finite strain poro-elastoplastic models. For real-world applications, we consider different scenarios where the reservoir is either mechanically weaker or stronger than the caprock. We also apply our workflow to incorporate complex geological structure, petro- physical heterogeneity, and multi-well schedule of a real oilfield. Our multiphase geomechanics formulation captures development of a field across different recov- ery mechanisms, from primary depletion to waterflooding to CO2 flooding, and predicts the geomechanical effects of fluid transport and storage on fault stability. We show that our workflow is able to characterize geomechanical signatures of flow compartmentalization and bypassed oil due to faults, inadequate well pattern configurations, and well rate variations. v Contents Acknowledgments ii Abstract iv List of Tables viii List of Figures ix 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Scope of Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Mathematical Model 5 2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Balance laws . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Poroplastic deformation . . . . . . . . . . . . . . . . . . . . 10 2.1.3 Fault slip problem . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Constitutive relations . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Flow rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2 Drucker-Prager model . . . . . . . . . . . . . . . . . . . . . 23 2.3 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.1 Weak Forms . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.2 Mechanics model . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.3 Flow model . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.4 Sequential iterative solution strategy . . . . . . . . . . . . . 37 2.3.5 Simulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4 Benchmark Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.1 The Terzaghi problem . . . . . . . . . . . . . . . . . . . . . 39 2.4.2 The Mandel problem . . . . . . . . . . . . . . . . . . . . . . 40 2.4.3 The rigid footing problem . . . . . . . . . . . . . . . . . . . 42 vi 3 Coupled Multiphase Flow and Geomechanics with Elastic mate- rials 47 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Physical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.1 Mesh initialization and boundary conditions . . . . . . . . . 51 3.2.2 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.3.1 Vertical and lateral positions of the hypocenter . . . . . . . 57 3.3.2 Stress path analysis . . . . . . . . . . . . . . . . . . . . . . . 60 3.3.3 Production-induced fault slip . . . . . . . . . . . . . . . . . 62 3.3.4 Injection-induced fault slip . . . . . . . . . . . . . . . . . . . 65 3.3.5 Characteristics of fault activation . . . . . . . . . . . . . . . 66 3.3.6 The slip onset time . . . . . . . . . . . . . . . . . . . . . . . 75 3.4 A real-world application . . . . . . . . . . . . . . . . . . . . . . . . 76 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4 Coupled Multiphase Flow and Geomechanics with Elastoplastic materials 84 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.2 Effect of finite strain plasticity on the flow problem . . . . . . . . . 86 4.3 Analogy between fault slip and plastic yielding . . . . . . . . . . . . 87 4.4 Production-induced plastic failure and fault slip . . . . . . . . . . . 88 4.4.1 Spatial profiles and time evolution . . . . . . . . . . . . . . . 90 4.4.2 Effect of plastic failure on fault slip . . . . . . . . . . . . . . 96 4.4.3 Effect of rock type . . . . . . . . . . . . . . . . . . . . . . . 98 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5 Role of Well Operations and Multiphase Geomechanics on Fault Stability 103 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2 Farnsworth Unit Geology . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3 Numerical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.3.1 Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.3.2 Fault model . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.3.3 Petrophysical properties . . . . . . . . . . . . . . . . . . . . 112 5.3.4 Initial and boundary conditions . . . . . . . . . . . . . . . . 115 5.3.5 Production and injection history . . . . . . . . . . . . . . . . 116 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.4.1 Field results . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.4.2 Fault results . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . 127 vii Reference List 132 viii List of Tables 3.1 Input parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.1 Four cases of reservoir-caprock system defined in terms of their elas- tic moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 ix List of Figures 2.1 Model problem of well-induced deformation in a faulted reservoir- caprock system. Well operation causes a change in the pore pressure (|Δp|> 0) and a change in the state of stress (|Δσ|> 0) inside the permeable reservoir and a change in the state of stress in caprock and basement. If the stress change is large enough, it can induce plastic/permanent deformation in the rock. If the fault is critically- stressed (depends on fault’s dip, cohesion and friction, and tectonic compression), thestresschangecaninducefaultslipatthehypocenter. 5 2.2 Kinematics of a porous body saturated with two fluids. ϕ is the solid deformation mapping. Mappings of the two fluid phases are not shown for brevity. Multiplicative elastoplasticity is implied by the decomposition of the deformation gradientF into elasticF e and plasticF p parts and the stress-free intermediate configuration. . . . 6 x 2.3 Flow chart showing the workflow and connections between different modules. Available seismic and geological information are used to construct a geomechanical model and populate it with poromechan- ical properties. Reservoir flow properties (porosity and permeabil- ity) and fluid properties (fluid density, viscosity and compressibil- ity) are obtained from an older uncoupled model and transferred to the geomechanical model using the Matlab Reservoir Simulation Toolbox (MRST (Lie, 2016)). The flow-geomechanical simulator, PyLith-GPRS, outputs 3D fields of the effective stress tensor, the displacement vector, fluid phase pressures and saturations. . . . . . 38 2.4 Terzaghi problem. (a) Model geometry with boundary conditions. The observation point is shown in a black circle. (b) Comparison of dimensionless pressure vs. time from the numerical simulation and analytical solution. (c) Comparison of displacement vs. time from the numerical simulation and analytical solution. . . . . . . . . . . 39 2.5 Comparison of finite strain simulation and infinitesimal strain simu- lation under the cases of small and large deformation. (a) Displace- ment vs. time (b) Dimensionless pressure vs. time. . . . . . . . . . 40 2.6 The Mandel problem. (a) Model geometry with boundary condi- tions. The observation point A is shown in a black circle. The observation line A’A is shown in dash line. (b) Comparison of dimensionless pressure vs. time from the numerical simulation and analytical solution. (c) Comparison of displacement vs. time from the numerical simulation and analytical solution. . . . . . . . . . . . 41 xi 2.7 Time evolution of the (a) vertical stress, (b) vertical strain and (c) evolution of the stress state at the observation point A in Fig. 2.6a. Verticalplasticstrain(bluedotline)startstoincreaseandtheelastic simulation (EL, dash line) shows higher vertical stress and lower vertical strain than Drucker-Prager plastic simulation (DP, solid line) after stresses reach the yield line (blue line). . . . . . . . . . . 43 2.8 Time evolution of (a) vertical displacement and (b) dimensionless pressure at the observation point A. Elastic simulation (EL, dash line) shows lower vertical displacement and higher pressure than Drucker-Prager plastic simulation (DP, solid line) after plastic fail- ure initiates. (c) Profile of vertical strain along AA’ at t=0.02 day. Vertical strain decreases from point A to A’. Total vertical strain from DP simulation (solid line) is higher than that from EL sim- ulation (dash line) in the cells that have plastic failure (dot line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 2.9 The rigid footing problem. (a) Model geometry with boundary con- ditions. The observation point, shown as a red dot, is below the footing. (b) Comparison of finite strain and infinitesimal strain simulations for two types of rocks: weak (E = 127 MPa) and strong (E = 254 MPa). (c) Number of sequential iterations for conver- gence between flow and mechanics problems in finite and infinitesi- mal strain simulations for weak rock (E = 127 MPa) . . . . . . . . 44 2.10 Results from finite strain (a to c) and infinitesimal strain (d to f) simulations. (a) Overpressure, (b) vertical displacement (exagger- ated by a factor of 10), and (c) plastic strain magnitude at t = 0.2 day. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 xii 3.1 A sketch of the faults and reservoirs in the St. Johns Dome nat- ural CO 2 reservoir located in a faulted anticline. CO 2 migrates from magma across the Coyote Wash Fault (CWF) to the reservoir and the surface, where travertines are formed. This is Figure 1(b) from (Miocic et al., 2019). Used with permission. Licensing infor- mation is available at (CC, 2019). . . . . . . . . . . . . . . . . . . . 50 3.2 (a) Physical model with a well injecting or producing from a reser- voir (gray color) confined between caprock and basement (white color) and intersected by a fault. (b) A vertical cross section of the model showing the reservoir, well and fault. (c) The fault surface viewed along the x-axis direction. Projection of the well on the fault is shown. Position of a line where detailed analysis of flow-induced stress and slip will be carried out in the Results section below. . . 52 3.3 Left: Relativepermeabilitiesofoil(k ro )andwater(k rw )asfunctions of the water saturation S w . Right: The formation volume factor of oil (B o ) and the oil viscosity (μ o ) as functions of the oil phase pres- sure p o . The oil compressibility is defined as c o = (1/B o )dB o /dp o . These properties affect the evolution of pressure and, therefore, that of induced stress. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 xiii 3.4 (a) Relative permeability of waterk rw increases and relative perme- ability of sCO 2 k rg decreases with increasing water saturation S w . The capillary pressure P gw = p g −p w decreases with S w because the rock is assumed to be water-wet. (b) The sCO 2 formation vol- ume factor B g , which is the ratio of the sCO 2 volume at reser- voir conditions to the volume at surface conditions, decreases with increasing sCO 2 pressure p g , due to sCO 2 compressibility: c g = (1/B g )dB g /dp g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.5 Twocasesofthemodelproblemstudiedinthispaper. Anellipsoidal reservoir is confined within a low permeability caprock and inter- sected by a sealing fault. Injection-induced expansion of the reser- voir applies compression on the fault, i.e. Δσ n < 0, and production- induced contraction applies tension, i.e. Δσ n > 0. . . . . . . . . . . 56 3.6 Oilproduction-inducedchangesintheeffectivenormaltraction Δσ 0 n , updip shear traction Δτ, fault pressure Δp f , and updip slip plotted over the fault surface. The plots are at t = 800 day, which is 200 days after the first event time. Slip nucleates in the top reservoir layer (asterisk denotes hypocenter) and propagates laterally along the layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.7 sCO 2 injection-inducedchangesintheeffectivenormaltraction Δσ 0 n , updip shear traction Δτ, fault pressure Δp f , and downdip slip plot- ted over the fault surface. The plots are at t=264 day, which is 169 days after the first event time. Slip nucleates in the bottom reser- voir layer (asterisk denotes hypocenter), propagates laterally in the reservoir, then updip in the reservoir and downdip in the basement. 58 xiv 3.8 Stress paths (a,b) and shear stress drop (c,d) at six points on the fault as CO 2 injection or oil production continues. Negative σ 0 n indicates compression. Negative τ indicates downdip shear. μ s and μ d lines are the Mohr-Coulomb failure lines corresponding to the static and dynamic coefficients of friction as per the slip-weakening model. The fault surface inset in (b) shows the locations of the points and the projection of the well. (c) and (d) show the drop in shear stress with increasing slip, which is related to the amount of seismic/aseismic energy released. . . . . . . . . . . . . . . . . . . . 63 3.9 3D fields of (a) overpressure Δp, (b) oil saturationS o , (c) horizontal displacement U x , and (d) vertical displacement U z resulting from 1000 days of oil production. Well location is shown with a black line in (c) and (d). Δp is almost zero (yellow color) in the footwall block due to the low permeability of the fault. S o decreases due to oil production. Cut-away views ofU x andU z fields are shown using a slice along the fault surface and a y = 1000 m slice. Production- induced contraction of the reservoir causes negative (blue color) U x displacements on the fault in the reservoir interval, negative U z in the overburden above the reservoir, and positive (red color) U z in the basement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.10 Depth profiles of changes in tractions and slip on the fault due to oil production. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 xv 3.11 3D fields of (a) overpressure Δp, (b) sCO 2 saturation S g , (c) hori- zontal displacement U x , and (d) vertical displacement U z resulting from 264 days of sCO 2 injection. Well location is shown with a black line in (c) and (d). Δp is almost zero (blue color) in the footwall block due to the low permeability of the fault. The sCO 2 saturation S g builds up near the crest of the dome due to the buoyancy and mobility of the injected sCO 2 . Cut-away views of U x and U z fields are shown using a slice along the fault surface and a y = 1000 m slice. Injection-induced expansion of the reservoir causes positive (red color) U x displacements on the fault in the reservoir interval, positiveU z intheoverburdenabovethereservoir, andnegative(blue color) U z in the basement. . . . . . . . . . . . . . . . . . . . . . . . 79 3.12 Depth profiles of changes in tractions and slip on the fault due to sCO 2 injection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3.13 Evolutions of the effective normal traction, σ 0 n (dash lines), and the updip shear traction, τ (solid lines), as functions of the change in average reservoir pressure, Δp res , for production and injection cases. Evolutions are shown at the respective hypocenter locations shown in Fig. 3.8 (circle for production and diamond for injection). t s1 marks the slip onset times for the two cases. The linear behavior in the tractions before t s1 is predicted by Eq. (3.14). . . . . . . . . . . 80 xvi 3.14 Characteristiccurvesofinducedfaultslip. (a)ThechangeinCoulomb failure stress ΔCFS and (b) the change in fault stress ratio ΔFSR as functions of the dimensionless reservoir pressure perturbation Δp rd . The evolutions are plotted at the hypocenters of production-induced (circle) and injection-induced (diamond) slip. Gray region indi- cates the region of induced fault destabilization. (a) In the gray region, before slip, ΔCFS at the event hypocenters evolve approx- imately linearly. (b) In the gray region, before slip, ΔFSR shows nonlinear behavior for injection and approximately linear behav- ior for production. For both injection- and production-induced slip, ΔFSR(t s1 ) =−(μ s +FSR 0 ), which is a constant value because FSR 0 is independent of the hypocenter depth. . . . . . . . . . . . . . . . . 81 3.15 ΔFSR vs. ΔCFS characteristic curves of induced slip for (a) super- critical CO 2 injection and (b) oil production. The symbols–square, triangles, circle, cross and plus–denote different locations on the fault (see the inset figure). Gray region is the region of induced fault destabilization. The curves start at (0, 0) and evolve with time along a straight line of negative slope, with deviations from the line indicating how fast the pressure is changing in the vicinity of that point. ΔCFS after failure is equal to the negative of initial CFS, which is different for points at different depths. . . . . . . . . 82 xvii 3.16 (a), (b), and (c) show the co-evolution of FSR and CFS on three faults in the Farnsworth oilfield. The grayscale of the dots, with the color bar on the right, and the arrows indicate the direction of time. The evolutions are plotted on fault points at the top of the oil reservoir, which are identified as blue dots in (d). The ellipse in (a) highlights the cloud of points during 21300-22350 days when ΔCFS is relatively constant. . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.1 (a) Physical model with a well producing from a reservoir (gray color) confined between caprock and basement (white color) and intersected by a fault. (b) A vertical cross section of the reservoir showing the layered topography, well and fault. Position of the Top point and Hanging wall cell, Footwall cell, where detailed analysis of stress and slip will be carried out in the Results section below. . 89 4.2 (a) Relative permeability of water k rw (solid) and relative perme- ability of oil k ro (dash) as functions of the water saturation s w . (b) The oil formation volume factor Bo and the oil viscosity μ o . . . . . 90 4.3 3D fields of (a) pressure change and (b) water saturation change at day850. Pressurechangeisalmostzero(yellowcolor)inthefootwall block due to the low permeability of the fault. Water saturation increase due to oil production. . . . . . . . . . . . . . . . . . . . . 91 4.4 3D fields of (a) vertical displacement, (b) plastic strain magnitude in the footwall block, and (c) slip magnitude on the fault surface at day 850. (a) shows a cut-away view of the u z field using a slice along the fault surface and a y = 1000 m slice. . . . . . . . . . . . 92 xviii 4.5 Case1stresscomponentsevolutionof(a)effectivenormalstress, (b) and (c) deviatoric stress versus pressure for selected footwall block cell (black line) and corresponding hanging wall block cell (blue line). (d) Stress paths of the (a) footwall block cell (black points) and hanging wall block cell (blue points) from Case 1. . . . . . . . 93 4.6 Time evolution of (a) the average reservoir pressure. Evolution of (b)theverticaldisplacementatthetopofthereservoirand(c)atthe top of the model, (d) plastic strain magnitude at the footwall block cell, and(c)slipareaonthefaultsurfaceversustheaveragereservoir pressure for the four cases with finite strain and infinitesimal strain simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.7 Stress paths of the (a) footwall block cell and (b) top point on the fault surface from Case 1 finite strain and infinitesimal strain simulations. The solid line in each plot is the M-C failure line. The arrows indicate the direction of time. The red star marks the onset of plastic failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.8 Effect of rock type. Evolution of (a) effective stress invariantI 0 1 and (b) deviatoric stress invariant √ J 2 in the footwall cell, and evolution of (c) effective normal fault traction σ 0 n and (d) shear traction τ at the top point as functions of the average reservoir pressurep res from thefourcasesinTable4.1. Thearrowsindicatethedirectionoftime. We highlight the initiation of plastic failure with red star and fault slip with red pentagon in Case 1. . . . . . . . . . . . . . . . . . . . 99 xix 4.9 Effect of rock type on the rate of failure. Evolution of (a) yield functionf p and(b)CoulombFailureFunctionCFF, asafunctionof average reservoir pressure p res from four cases. The arrows indicate the direction of time. . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.10 Stress paths of the (a) footwall block cell and (b) top point on the fault surface from the four cases. The solid straight line in (a) is the D-P failure line and in (b) is the M-C failure line. The arrows indicate the direction of time. . . . . . . . . . . . . . . . . . . . . . 101 5.1 Map of the Farnsworth Unit (FWU) field. Left: The border of FWU and some of the well locations are shown. Right: FWU is located in the northern Texas within the Anadarko Basin (shaded region), and it is shown within a green square. The main source of anthropogenic CO 2 used during injection is highlighted by two star- shaped symbols: the northern star is the Arkalon Ethanol Plant in Liberal, Kansas, and the southern star is the Agrium Fertilizer Plant in Borger, Texas. . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.2 Map view showing the faults within FWU at the Morrow Sandstone depth level. The FWU boundary is shown with thick black lines. Township and range boundaries are shown in a light gray color. Nine faults are interpreted using FWU’s 3D seismic survey. Except Fault-6, all other faults intersect both the Morrow-B Sandstone and Thirteen Finger Limestone. . . . . . . . . . . . . . . . . . . . . . . 106 xx 5.3 (a) Top view and (b) side view along cross-section AA’ show the model geometry and the mechanical boundary conditions. Fault-1 and Fault-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 towards south with a 6 ◦ angle as shown in (b). Three layers are highlighted in (b): Morrow Shale, Morrow-B Sandstone Top and Morrow-B Sandstone Base. (a) The minimum principal stress is applied on the y-positive surface. The maximum principal stress is applied on the x-positive surface. (b) A normal compression equal to the overburden is applied on the top boundary. . . . . . . 108 5.4 Computational domain mesh showing the locations of the fault sur- faces and the wells. We use tetrahedral elements with smaller-sized elements within the reservoir layers and next to the fault surfaces and larger-sized elements near the domain boundaries. 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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.5 3D mapping of rock properties from an existing flow-only uncoupled model into the geomechanical model. Histograms of log-normalized permeability in (a) x-direction, (b) y-direction, and (c) z-direction from the uncoupled Eclipse model (blue color) and the coupled geomechanical model (orange color). . . . . . . . . . . . . . . . . . 113 xxi 5.6 Formation volume factors for (a) oil, (b) water, (c) gas phases. The bubble point pressure of oil is approximately 1900 psi (1 psi = 6894.76 Pa) below which dissolved gas is released from the oil phase. We assume a constant water compressibility, which corresponds to a constant dB w /dp. We use CO 2 properties for the gas phase instead of the reservoir gas properties. This allows us to honor the CO 2 injection and transport behavior in the field more accurately. . . . 114 5.7 The three-phase (oil, water and gas) relative permeability model in oursimulationisbasedontwotwo-phaserelativepermeabilities. (a) Relative permeability of water k rw (dash) and relative permeability of oil in presence of waterk row (solid) as functions of the water sat- urationS w . (b) Relative permeability of gask rg (dash) and relative permeability of oil in presence of gas k rog (solid) as functions of the free gas phase saturation S g . The relative permeability of oil k ro is calculated dynamically during the simulation as a function of k row , k rog , and the saturations. . . . . . . . . . . . . . . . . . . . . . . . . 115 5.8 TypicalproductionandinjectionrateprofilesofFWUwells. (a)The liquid production rate (Qo+Qw, oil+water) of a producer, (b) the water injection rate Qw of a water injector, and (c) the gas injection rate Qg of a gas injector. Production starts at the beginning of the simulation. Water injection starts in year 1995. Gas injection starts in year 2010. Note the variability and discontinuity in the rate profiles, which will be related to fault stability patterns in the analysis below. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 xxii 5.9 (a) Time evolution of the average reservoir pressure (psi) due to production,injectionandflow-inducedstresses. Themaximumpres- sure drop is approximately 200 psi. Cumulative liquid production is shown by the red dot curve. (b) Time evolution of the vertical displacement at two selected points on Morrow-B Top (blue line) and Morrow-B Base (green dash line). Vertical displacement on the two surfaces have opposite signs because of their opposite motion in response to volumetric contraction/expansion of the reservoir with pressure decrease/increase. Displacement at the reservoir bottom has smaller magnitude due to the zero displacement boundary con- dition at the bottom boundary of the model. . . . . . . . . . . . . 118 5.10 Changes in the reservoir pressure compared to the initial pressure at the selected time steps. Pressure decreases in regions with net pro- duction and increases in regions with net injection. The maximum pressure drop is 200 psi approximately and occurs around t=13250 day. The final pressure drop att=22250 day is 150 psi approximately.119 5.11 Changesintheoilsaturationwithrespecttotheinitialoilsaturation at the selected time steps. The oil saturation decreases around both producers and injectors, as expected. The maximum change in the saturation is approximately 0.5. The oil saturation remains almost unchanged in regions far from the wells, which suggests the presence of bypassed oil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 xxiii 5.12 Vertical displacement field (cm) at the selected time steps. Pressure drop in the reservoir leads to negative (blue color) vertical displace- ment above the reservoir while pressure increase in the reservoir leads to positive (red color) vertical displacement above the reser- voir. Below the reservoir, the vertical displacement shows opposite behavior as it experiences uplift during reservoir contraction and downward push during reservoir expansion. The asymmetry in the mechanical boundary conditions between the top (fixed traction) and bottom (fixed displacement) boundaries of the model affects the magnitude of vertical displacement. . . . . . . . . . . . . . . . . 121 5.13 Change in the fault pressure (psi) along Fault-1 at the selected time steps. A trend similar to that of Figure 5.10 can be observed. The fault pressure in the reservoir depth interval (appearing as a hori- zontal stripe) declines at early times due to production from both fault blocks. Later, it increases due to injection and then decreases again towards the end of the simulation due to net production. . . . 122 5.14 Effect of spatial distribution and temporal evolution of production on fault stability. Changes in (a) the effective normal traction and (b) the updip shear traction on Fault-1 at t=4500 day with respect to the initial traction values. Liquid production profiles of active producers near Fault-1 are shown in (c) through (g). Surface loca- tions of the wells with respect to the fault are shown in (h). . . . . 124 xxiv 5.15 Effect of spatial distribution and temporal evolution of injection on fault stability. Changes in (a) the effective normal traction (MPa) and (b) the updip shear traction (MPa) on Fault-1 at t=5000 day. Water injection in six wells, shown in (c) through (f) along the fault strike, causes updip shear above the reservoir (red color in (b)) and downdip shear below the reservoir (blue color in (b)). Ten- sile stresses are induced in the overburden and basement, which cause tensile changes in the effective normal traction (red color in (a)). Production dominates injection in the south, which causes compressive changes in the effective normal traction (blue color in (a)) near well 2. Surface locations of the wells with respect to the fault are shown in (g). . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.16 Changes in the Coulomb failure stress (MPa) on Fault-1 at the selected time steps. Red color or positive values indicate induced destabilization. Blue color or negative values indicate induced sta- bilization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.17 Evolution of ΔCFF on selected points within the overburden (b) and reservoir (c) intervals on Fault-1, Fault-3L and Fault-3R. Panel (a) shows the locations of the points in the reservoir interval (blue dots) and the overburden interval (red dots). . . . . . . . . . . . . . 130 xxv 5.18 Comparison of ΔCFF (solid) and reservoir pressure (dotted and dashed) evolutions on (a) Fault-1, (b) Fault-3L and (c) Fault-3R at points within the reservoir interval. The two pressure curves are from two elements connected to the blue dot (inset) on the two sides of the fault. ΔCFF follows the trend of the larger pressure because p f is defined from the larger of the two pressures. Plots (d), (e) and (f) show the characteristic ΔCFF versus Δp rd behavior for Fault- 1, Fault-3L and Fault-3R, respectively. Time evolution is shown through the color of the dots. The color legend on the right shows that the time evolves fromt=0 (black) tot=22250 day (gray). Start of water injection and CO 2 injection are marked in (d). . . . . . . . 131 xxvi Chapter 1 Introduction 1.1 Background We know that injection or production of fluids in a reservoir can perturb its mechanical equilibrium with the surrounding rock by altering the state of stress insideandoutsidethereservoir(Raleighetal.,1976;Yerkes&Castle,1976a;Segall, 1989a). Such stress changes can be significant enough to induce slip on the faults thatarefavorablyorientedandclosetofrictionalfailure. Theeffectofhydrocarbon production (Segall et al., 1994; Sanz et al., 2015; Juanes et al., 2016; v. Thienen- Visser & Fokker, 2017), water production (Gonzalez et al., 2012; Hornbach et al., 2015a; Albano et al., 2017a) and injection (Ellsworth, 2013; Keranen et al., 2013; Gan & Frohlich, 2013; McGarr, 2014; Schoenball et al., 2012; Catalli et al., 2013; Dengetal.,2016;Albanoetal.,2017a;Leietal.,2017;Horton,2012;Hosseinietal., 2018)aswellasgeologicCO 2 sequestration(Morrisetal.,2011;Jha&Juanes,2014; Zhao & Jha, 2019) on fault and fracture activation has attracted much attention recently. Several studies (Cappa & Rutqvist, 2011a,b; Jha & Juanes, 2014; Segall & Lu, 2015b; Juanes et al., 2016) have been conducted with simple physical models of induced seismicity consisting of a confined reservoir, permeable or impermeable faults, and injection or production wells to understand the physical mechanisms behind such events. The contact nonlinearity caused by fault surfaces has beed addressed in some of the previous studies (Chang & Segall, 2016a,b). Most studies have limitations 1 such as they focus in a 2D domain with a flat-layered reservoir and a 1D fault located in the basement rocks with a constant friction coefficient. The seismicity rate is computed a posteriori by substituting the pore pressure and stress changes in Dieterich’s model (Dieterich, 1994), which decouples the fault response from the dynamics of flow and mechanical deformation and also neglects stress interactions amongmultiplehypocentersthatemergesuccessivelyontherupturingfaultsurface (Heimisson, 2019). Meanwhile, most of the previous studies assume a linear elastic infinitesimal deformation response of the rock hosting the fault. We know that deformation response of many rocks, such as clayey sands, ultra-low permeability shales, and carbonates can become viscoelastic or elstoplastic at large effective stresses typ- ically encountered during long-term production or injection (Maury et al., 1996; Bennett et al., 2015; Davis & Selvadurai, 2002). Hard brittle rocks also show inelastic stress-strain behaviors through microcracking activity (Scholz, 1968). Long-term groundwater withdrawal in California’s Central Valley (Ireland et al., 1984; Galloway & Riley, 1999) and oil and water production in Wilmington (Allen, 1968; Allen & Mayuga, 1969) and Gulf of Mexico (Chan & Zoback, 2007) oilfields are examples where both large strains (e.g., > 2%) and plastic deformation are expected to occur. Plastic deformation also occur at relatively small values of the induced stresses if the rock is mechanically weak or ductile, e.g. chalk reservoirs in Ekofisk (Chin & Thomas, 1998), some shales (Sone & Zoback, 2011; Chang et al., 2014), and clay beds or aquitards (Gu & Ran, 2000). Presence of clay particles or interbeds in sandstone aquifers can give rise to a lag in the mechanical response of aquifers under depletion (F.G. Bell, 1986; Shi et al., 2007), which manifests as plasticity. In general, long-term fluid extraction or injection activities can acti- vate nonlinearity in both the material constitutive response and the kinematics of 2 deformation. Use of elastic and infinitesimal deformation models in those cases will predict inaccurate stresses and, therefore, inaccurate onset, magnitude, and dura- tion of the fault slip events. This is especially problematic for smaller magnitude earthquakes or aseismic slip events (Mw<4), which are not felt by humans but may cause damage to surface facilities and leakage of reservoir fluids (or other types of stored fluids such as natural gas or CO 2 ) along the activated fault into freshwater resources. Other manifestations of plastic deformation include drop in the reser- voir permeability, loss of well productivity/injectivity, and wellbore failure (Allen, 1968; Fredrich et al., 2000). Accounting for the above nonlinearities can improve the agreement between model-predicted seismicity and ground deformation (Gonzalez et al., 2012; Bourne & Oates, 2017) and the observed data. Therefore, there is an urgent need to develop simulation methods that resolve coupled flow, geomechanics, and dynamic faulting processes in plastically deforming rocks and demonstrate their application to real reservoirs. Here, we develop a coupled flow-geomechanics-faulting method. Applying the method to faulted reservoir-caprock systems with realistic geological structure, heterogeneous rock properties, multiphase flow and time-dependent and spatially distributed production and injection is another challenge that we tackle. 1.2 Scope of Work In Chapter 2, we present the mathematical problem of coupled multiphase flow, finite strain plasticity, and dynamic fault slip. We discuss the numerical model based on a finite element discretization of the mechanics problem and a finite volume discretization of the flow problem. The model is tested on benchmark consolidation problems. 3 In Chapter 3, we conduct a coupled flow-geomechanics-faulting study that investigates fault behavior during CO 2 injection and oil production. We analyze fault tractions, slip, and reservoir pressure to understand how a fault is reacti- vated in successive steps and how the reactivation process selects the hypocenter locations within and outside the reservoir. We extract characteristic features of induced reactivation in terms of the Coulomb failure stress, the Fault Stress Ratio, and the reservoir pressure perturbation caused by injection or production. In Chapter 4, we applying the finite strain elastoplasticity method to faulted reservoir-caprock systems with oil production. Results from simulations with dif- ferent elastic moduli for reservoir and caprock layers are analyzed to understand the impact of plasticity on characteristics of fault slip in different rock types. In Chapter 5, we apply our modeling and simulation frameworks to a real field. We monitor the evolution of flow and mechanical behavior of the overburden- reservoir-basement storage complex during CO 2 injection and oil production to understand the evolution of pressure, oil saturation and stress in the storage com- plex and to monitor the mechanical stability of reservoir-intersecting faults. 4 Chapter 2 Mathematical Model Consider the deformation of an oil reservoir, confined between caprock and basement and intersected by a fault, due to oil production from a well drilled into the reservoir. See Figure 2.1 for a schematic of the conceptual problem. The porous medium has three phases: the solid skeleton, which is made up of solid grains and pores, and the two fluid phases that occupy the pores (Coussy, 2004). Now, let’s consider the deformation of a representative elementary volume (REV) V 0 ⊂ R 3 of the reservoir (Figure 2.2). Γ 0 is the boundary of the REV. We denote the reference configuration of the volume by the solid particle position vectorX∈R 3 andthefluidparticlepositionvectorsX o andX w , wheresubscripts Reservoir Fault well Sandstone Hypocenter Overburden Basement Caprock |Δσ| , |Δp| > 0 |Δσ| > 0 |Δσ| > 0 Tectonic compression Clay interbeds Figure 2.1: Model problem of well-induced deformation in a faulted reservoir- caprock system. Well operation causes a change in the pore pressure (|Δp| > 0) and a change in the state of stress (|Δσ|> 0) inside the permeable reservoir and a change in the state of stress in caprock and basement. If the stress change is large enough, it can induce plastic/permanent deformation in the rock. If the fault is critically-stressed (depends on fault’s dip, cohesion and friction, and tectonic compression), the stress change can induce fault slip at the hypocenter. 5 reference conguration deformed conguration water oil intermediate conguration Figure 2.2: Kinematics of a porous body saturated with two fluids. ϕ is the solid deformation mapping. Mappings of the two fluid phases are not shown for brevity. Multiplicative elastoplasticity is implied by the decomposition of the deformation gradient F into elastic F e and plastic F p parts and the stress-free intermediate configuration. o and w refer to oil and water phases. Let ϕ : V 0 × [0,t]→ R 3 be the solid deformation mapping, where t denotes time. After deformation, the current or deformed configuration is given by the position vectorx, which is occupied by all three phases simultaneously. The deformed volume is V within the boundary Γ. The solid phase deformation vector isu =x−X, and the solid phase deformation gradient tensor is F =∇ X x = 1 +∇ X u, where∇ X denotes partial derivatives with respect toX and 1 is the rank-2 identity tensor. We will use∇ =F −T ∇ X to denote partial derivatives with respect to current coordinatesx. The superscript −T indicates transpose of inverse. Let M = M w +M o be the total fluid mass in the reference configuration per unit bulk volume in the reference configuration, where M w = ρ w φs w J is the water mass content, M o = ρ o φs o J is the oil mass content, φ is the true or Eulerian porosity defined as the pore volume in the current configuration over the bulk volume in the current configuration, ρ w andρ o are water and oil densities in the current configuration,s w ands o are water and oil saturations in the current configuration with respect to the current pore volume, and J = detF = dV/dV 0 is the Jacobian of solid deformation. The saturations 6 must sum to 1 by definition: s w +s o = 1. The Lagrangian porosity (current pore volume per unit reference bulk volume) is defined as Φ =φJ. As commonly considered in the literature (Sanavia et al., 2002), we take the solid motion as the reference motion and describe the fluid motion relative to the solid motion. We have the following strain and deformation tensors: E = (1/2)(F T F− 1) = (1/2) (∇ X u +u∇ X +∇ X u·u∇ X ), (2.1) ε = (1/2) (∇u +u∇) = sym[∇u], (2.2) b =FF T , (2.3) C =F T F (2.4) whereE is the Green-Lagrange strain tensor,ε is the infinitesimal strain tensor, b is the left Cauchy-Green deformation tensor and C is the right Cauchy-Green deformation tensor. Here, sym denotes the symmetric operator. The linearized volumetric strain in the current configuration is ε v = tr[ε] = logJ. This yields J = exp (ε v ), which is approximated asJ = 1 +ε v in the infinitesimal deformation theory. As fluid is produced from or injected into the reservoir, principles of conservation of mass and momentum of the solid and fluid components govern the evolution of state variables of the system such as stress, strain, fluid mass content, and fluid pressure. 2.1 Governing equations 2.1.1 Balance laws Mass conservation We assume quasi-static mechanical equilibrium. The mass content of the solid phase is conserved. This will be utilized later to obtain the 7 porosity evolution equation. The conservation equation of the water phase mass during deformation, written in the reference configuration of the porous medium, reads as . M w +J∇·w w +Jq w = 0, in V 0 (2.5) where . ( ) denotes the material time derivative,w w =ρ w v w is the water mass flux, q w isthewatermassproductiontermduetowells,v w =− (kk rw /μ w )·(∇p w −ρ w g) is the water volumetric flux or the Darcy velocity (obtained from the conservation of fluid momentum after neglecting inertia, drag, and viscous forces), p w is the excess water pressure (above the initial pore pressure),k is the rank-2 permeability tensor, μ w is the dynamic viscosity of the water, k rw is the relative permeability of water in presence of oil, and g is the acceleration vector due to gravity. The flux terms are in the current configuration. The permeability transforms as k = (1/J)Fk 0 F T , wherek 0 is the permeability tensor in the reference configuration. The water mass content in the current configuration ism w =M w /J. By definition, the water density in the reference configuration is ρ 0 w =Jρ w . Similar to Eq. (2.5), we have another conservation equation for the oil phase mass M o that introduces the Darcy flux of oil in terms of ρ o , oil viscosity μ o , oil relative permeability k ro , and oil phase pressurep o . Since oil and water are immiscible phases, the difference between their pressures is the capillary pressure P c = p o − p w , which is often modeled via a saturation-dependent capillary pressure relation, also known as the retention curve in air-water systems studied in soil hydrology. A non-zero capillary pressure can affect the elastic stiffness of the porous medium. Momentumconservation Giventhepossibilityoflargedeformationsandrota- tions in the finite strain theory, and the need to track the displacement of the solid 8 skeleton from its reference configuration, it is customary to write the conservation of linear momentum in the reference configuration as ∇ X ·P +Jρ b g = 0, ∀X∈V 0 (2.6) where P = JσF −T is the first Piola-Kirchhoff (PK) stress tensor (asymmetric tensor) defined in terms of the total Cauchy stress tensor σ (symmetric). σ is a measure of force acting on an area in the current configuration per unit area in the current configuration. P is a measure of force acting on an area in the current configuration per unit area in the reference configuration. Above, ρ b = (1−φ)ρ s +φ(ρ o s o +ρ w s w ) is the bulk density in the current configuration, and ρ s is the solid grain density. The second PK stress tensor, which measures force in the reference configuration per unit area in the reference configuration (symmetric tensor), is defined asS =JF −1 σF −T , which shows how the stress in the current configuration can be pulled back to the reference configuration. This also implies P =FS. WhenF' 1 andJ' 1,S can be replaced byσ, andE can be replaced by ε such that we recover the infinitesimal strain model. The Kirchhoff stress tensor in the current configuration is defined asτ = Jσ =FSF T . To facilitate constitutive modeling of fluid-saturated rocks, which is discussed in detail below, we define the effective stresses: the effective Kirchhoff stress τ 0 = Jσ 0 and the effective second PK stressS 0 =JF −1 σ 0 F −T . 2.1.2 Poroplastic deformation Following the poroelastoplastic formulation proposed earlier (Simo, 1992; Borja & Alarcon, 1995; Armero, 1999), we use a multiplicative decomposition ofF and 9 an additive decomposition of total fluid massM into their elastic and plastic parts: F =F e F p , M =M e +M p , (2.7) where the superscripts e and p denote elastic and plastic parts, respectively. The elastic part of left Cauchy-Green tensor is b e = F e F e,T . The elastic part of the right Cauchy-Green tensor isC e =F e,T F e and the plastic part isC p =F p,T F p . The volumetric plastic strainε p v =ε v −ε e v = logJ p withJ p = detF p . For isochoric plastic flow, a common assumption in plasticity, J p = 1 and J = detF e . 2.1.3 Fault slip problem To model seismic/aseismic slip on the fault, we treat the fault as a strong dis- continuity, i.e. the displacement field is discontinuous across the fault, as opposed to the weak discontinuity which implies only a discontinuous strain field. The case of a strong discontinuity emerging inside an initially continuous finite element has been treated in the finite strain theory (Borja, 2002; Callari & Armero, 2004) by enhancing the deformation gradient in the nonconforming elements with an additional component related to the displacement jump and enhancing the finite element test functions with a Heaviside function. Here, we consider the conform- ing case where the fault slip is constrained to lie along the edges between adjacent finite elements. The conforming case is practical in many geological settings where geometry of the rupture surface of a seismogenic fault is known from previous earthquakes, well logs, or basin modeling studies. The approach requires construc- tion of a mesh conforming to the fault geometry, which can be time-consuming in a 3D domain containing multiple intersecting or curved faults. However, simulation of the conforming case is expected to be faster than that of the nonconforming case 10 which has to solve for the rupture path as well. We treat the fault as a frictional interface Γ f (reference configuration) across which the displacement field can be discontinuous to accommodate seismic/aseismic slip. We define the fault slip vec- tord f in the current configuration as the displacement discontinuity or jump from the ‘−’ side (Γ + ) to the ‘+’ side (Γ − ) of the fault (e.g. the footwall and hanging wall sides of a dip-slip fault) as u + −u − =d f on Γ f . (2.8) We define ˆ n f as the unit normal vector on the fault (reference configuration) pointing from the -ve side to +ve side. The normal vector can vary in space for a fault that is either naturally arcuate or becomes curved due to reservoir expansion or contraction induced by fluid injection or production. Mechanical equilibrium requires continuity of the total traction vector across the fault: T + f = (FS) + ˆ n f = (FS) − ˆ n f =T − f . The effective fault traction vector is a projection of the effective second PK tensor stress tensor on the fault: l f = FS 0 ˆ n f . We define the local fault coordinate system using three mutually orthogonal unit vectors, (ˆ s 1 , ˆ s 2 , ˆ n f ), which are positive in the left-lateral, updip and opening directions of the fault. l f can be split into a shear traction vector and a normal traction vector as (τ s ,σ 0 n ˆ n f ) whereτ s = (l f · ˆ s 1 )ˆ s 1 + (l f · ˆ s 2 )ˆ s 2 lies on the 2D fault surface and σ 0 n =l f · ˆ n f . Similarly,d f can be split into the slip vectord s = (d f · ˆ s 1 )ˆ s 1 + (d f · ˆ s 2 )ˆ s 2 , which points in the rake direction, and the fault opening d n . 2.2 Constitutive relations Minimization of the free energy of a system, as the system evolves dynamically, is often used to extract the constitutive relations applicable to the system. Let ψ 11 be the free energy (per unit reference volume) of the elastoplastic system defined in terms of three internal variables: b e , M e , and a strain-like internal variable ξ which tracks the isotropic hardening response during plastic deformation. The rate of change of free energy can be written in terms of the stress power (τ :d, S : . E, orP : . F), the fluid accumulation power, the power of faulting produced by the fault traction vector in realizing the slip across the fault, and the internal dissipation: Z V 0 . ψdV 0 = Z V 0 τ :d +χ w . M w +χ o . M o −D dV 0 + Z Γ 0 f T f · . d f dΓ 0 (2.9) whereD = C e S : L p +h . ξ +χ w . M p w +χ o . M p o −P c Φ . s w ≥ 0 is the local internal dissipation (per unit reference volume) due to incremental plastic deformation, hardening, plastic fluid content, and capillarity (due to the surface energy asso- ciated with fluid interfaces in a multiphase system). Here, L p = . F p F p,−1 is the plastic distortion rate, andh is an internal stress-like variable related to the inter- nal strain-like variable ξ as h =−∂ψ/∂ξ. Above, we do not include dissipation due to Darcian fluid flow, which occurs along the fluid pressure gradient, inD. For reversible processes,D = 0. The last term above is the rate of work or power of faulting . W f , which can be approximated (Scholz, 2002) as . W f ≈ T avg f × . d avg f × Γ slip f , where T avg f = (T aft f + T bef f )/2 is the average fault stress magnitude, . d avg f = ( . d aft f + . d bef f )/2 is the average faultslipvelocitymagnitude, Γ slip f isthefaultsliparea, andthesuperscriptsbef and aft indicatebeforeandaftertheslipevent.W f canbedecomposedintothreeparts: the heat released during slip, energy spent during creation of new material surfaces e.g. in the fault process zone, and the seismic radiation energy E s . Ignoring the surfacecreationenergy, whichhasbeenfoundtobe3-4ordersofmagnitudesmaller 12 than the other terms, E s =η seis W f , where η seis = (T aft f −T bef f )/(T aft f +T bef f ) is the seismic efficiency. Based on data from many earthquakes, it has been found that η seis ≤ 0.06 (McGarr, 1999). The moment magnitude of the slip event can be obtained from the Hanks-Kanamori scale (Hanks & Kanamori, 1979), M w = 2 3 log 10 GW f T avg f ! − 10.7 (2.10) It is also important to note that for a given drop in the free energy in Eq. (2.9), ifD increasesduetoplasticdissipation, theenergyavailableforfaulting,W f , decreases. Our results in Section 4.3.2 provide further evidence for this hypothesis. Poromechanics The kinematic or geometric nonlinearity can arise either in the form of large displacement/large rotation and large strain or large displace- ment/largerotationandsmallstrain. Althoughbothcasesrequiretheuseofsecond Piola-Kirchhoff stress and Green-Lagrange strain tensors in the equilibrium and constitutive equations, modeling the former requires hyperelastic constitutive rela- tions with large strain constitutive tensors. The latter behavior can be modeled by simply substituting σ 0 and ε with S 0 and E, respectively. The small strain constitutive tensor does not change during this substitution (Bathe, 1995). In the finite strain theory, the requirement of frame indifference during material response means the material constitutive models are expressed in terms of the Lie derivative of the effective Kirchhoff stress,L v τ 0 . The relation between the total and effective stresses in the current configuration becomes L v τ 0 =L v τ +b . p e 1, (2.11) 13 whereL v τ = F . SF T is the Lie derivative of the total Kirchhoff stress,L v τ 0 = F . S 0 F T is the Lie derivative of the effective Kirchhoff stress,b is the Biot coefficient of the rock fully saturated with water or oil, and p e is the saturation-weighted equivalent pore pressure defined as p e = s w p w +s o p o (Borja, 2006). The Biot coefficient is a tensor in anisotropic media (Bubshait & Jha, 2019; Zhao & Borja, 2020). We assume tensile stresses are positive. In the reference configuration Eq. (2.11) becomes . S 0 = . S +b . p e C −1 (2.12) The effective stresses are responsible for the deformation of the material. In the current configuration, the rate form of this relation is L v τ 0 =c :d, (2.13) where J −1 c is the rank-4 Eulerian elasticity tensor in the current configuration, d = sym[l] is the rate of deformation tensor, and l = ∇(∂u/∂t) = . FF −1 is the gradient of Eulerian velocity vector ∂u/∂t. The stress-strain relation in the reference configuration becomes . S 0 =C : . E where . E = (1/2) . C =F T ·d·F is the pull-back rule, andC is reference configuration’s elastic constitutive tensor related to c as c ijkl = 2F iA F jB F kC F lD C ABCD or F (C : . C)F T = c :d. The incremental form of the pull-back rule provides a relation between the finite and infinitesimal strains asδE =F T δεF. For a homogeneous isotropic linear elastic behavior, the elasticity tensor can be split as follows c = k dr − 2G 3 1⊗ 1 + 2GI =k dr 1⊗ 1 + 2GI 0 (2.14) 14 where k dr = (1/9)1 : c : 1 is the drained bulk modulus, I is the rank-4 identity tensor,I 0 =I− 1 3 1⊗1 is the rank-4 deviatoric operator tensor, andG is the shear modulus. We assume constant elastic moduli and therefore no coupling between volumetric and deviatoric responses (Borja, 2004). Also, we assume that c is independent ofs w andP c . The Biot coefficient can be expressed asb = 1−(k dr /k s ), wherek s is the solid grain modulus. Applying the trace operator to Eq. (2.11), we have L v τ 0 : 1 =L v τ : 1 + 3b . p e . (2.15) Using Eq. (2.13) and noting tr[d] =d : 1 = . ε v , we obtain (c :d) : 1 = 3k dr . ε v =F . SF T : 1 + 3b . p e , (2.16) which resembles the effective stress relation in the infinitesimal formulation once we define,L v τ v = 1 3 F . SF T : 1 andL v τ 0 v = k dr . ε v . The rate of deformation tensor can be split into volumetric and deviatoric parts: d = (1/3) . ε v 1 + . e, where e is the deviatoric strain tensor. The Green-Lagrange strain tensor can be split as . E = (1/3) . ε v F T F +F T . eF. In large displacement, large rotation, but small strain analysis (strain magni- tude less than a couple percentages), material models developed for the infinitesi- mal strain analysis can be used directly, provided thatσ 0 is substituted with the second PK stress tensorS 0 , andε is substituted with the Green-Lagrange strain tensor E (Bathe, 1995). This is possible because S 0 and E, which are energy- conjugates, are invariant under rigid body motions (e.g., rotations of the grid cells of a reservoir/caprock layer deforming under fluid production). We apply the trace operator to Eq. (2.12) to write tr[ . S 0 ] = tr[ . S] +b . p e tr[C −1 ]. Defining . S 0 v = tr[ . S 0 ], 15 . P e = . p e tr[C −1 ], and splitting the strain tensor as . E = (1/3) . E v 1 + . E d and the stress tensor as . S = . S v 1 + . s, we have . S 0 v =K dr . E v = . S v +b . P e , . s = 2G . E d (2.17) whereK dr is defined fromC, ands andE d are deviatoric stress and strain tensors. The first effective stress invariant can be defined asI 0 1 = 3S 0 v and the second invari- ant of the deviatoric stress can be defined asJ 2 = (1/2)s :s. Using Eq. (2.12), we define the fault tractions as follows: σ 0 n = FS +bp f C −1 : ˆ n f ⊗ ˆ n f , τ s =l f −σ 0 n ˆ n f (2.18) where p f is the fault pressure modeled from the two pressures on the two sides of the fault (Jha & Juanes, 2014; Yang & Juanes, 2018). Multiphase flow To express Eq. (2.5) in terms of the water phase pressure as the primary variable, we differentiate M w =ρ w s w Φ and apply the chain rule: . M w =ρ w (Φ(c w . p w s w + . s w ) + . Φs w ), (2.19) where we introduced the constitutive relation defining the water compressibility in the current configuration as c w = (1/ρ w )dρ w /dp w , which allows the water density to be expressed in terms of the water pressure, e.g. ρ w = e cwpw . This defines the chemical potential of water as χ w = R (1/ρ w )dp w = R pw 0 e −cwpw dp w . The . M o term in the oil mass conservation equation is similarly expressed in terms of the oil phase pressurep o by using the oil compressibilityc o . We also define the saturation- weighted Lagrangian fluid density: ρ 0 f = (ρ 0 w s w +ρ 0 o s o ). 16 We assume pressure-dependent fluid viscosities μ w and μ o . Additional consti- tutive relations are required fork rw ,k ro andP c . We use Corey-type relation for the relative permeabilities and van Genuchten-type relation for the capillary pressure, both of which are expressed as analytical functions of s w (the wetting phase satu- ration). We discuss the capillary relation in more detail in the flow rules section below. The constitutive equation of poroelasticity states that the fluid mass increment has two contributions, the volumetric expansion of the skeleton and the increment in fluid pressure (Coussy, 1995), . M w ρ 0 w = (b w . ε v +N ww . p w +N wo . p o ), . M o ρ 0 o = (b o . ε v +N wo . p w +N oo . p o ) (2.20) whereN = [N ww ,N wo ;N wo ,N oo ] =M −1 b is the 2× 2 inverse Biot modulus matrix in the current configuration. Denoting the Biot modulus matrix in the reference configuration as M 0 b , we have M b = J −1 M 0 b . b w = bs w and b o = bs o are the fluid-specific Biot coefficients. Components of the N matrix can be expressed in terms of rock and fluid properties (Lewis & Schrefler, 1998; Jha & Juanes, 2014). For example, for water-saturated conditions (s w = 1) or single-phase flow, M b = (φc w + (b−φ 0 )/k s ) −1 . Fault slip The normal force on the fault is always compressive except when the two sides loose contact at which time the normal force becomes zero. Along with the no-interpenetration condition, this implies the following Kuhn-Tucker conditions: σ 0 n ≤ 0, d n ≥ 0, σ 0 n d n = 0. We define the stick-slip behavior of the fault in terms of σ 0 n ,|τ s |, and the fault friction stress τ c −μ f σ 0 n , where μ f is the 17 coefficient of fault friction and τ c is the fault cohesion. Here,|·| is the magnitude operator or Euclidean norm of a vector. We use the Mohr-Coulomb (M-C) stick- slip condition to determine when a point on the fault slips. The condition is, f s =|τ s |−τ c +μ f σ 0 n ≤ 0, (2.21) i.e., there is no slip until the shear traction magnitude is less than the friction stress, otherwise there is slip and relaxation of fault stresses. f s is also known as the Coulomb Failure Function (CFF) in the induced seismicity literature (Jha & Juanes, 2014). Tracking the change in CFF, as injection or production continues in the reservoir, is a standard method of assessing and mitigating the induced seismicity hazard (Gonzalez et al., 2012; Zhao & Jha, 2019). Initial conditions of a simulation are chosen such that all points on the surface of the fault satisfy the above condition. As the reservoir pressure and stress change with the well operation,|τ s | and σ 0 n change on the fault. If the changes are destabilizing, e.g. increasing shear and/or decreasing compression, the condition is violated and the fault slips. We define the hypocenter as the fault node where the M-C condition is violated earliest in time. We use the slip-weakening relation to capture the rapid drop in μ f immediately after the onset of slip, which is commonly observed in earthquakes (Scholz, 2002). Slip weakening allows the fault friction coefficient to decrease from its static value μ s to a dynamic value μ d over the critical slip distance d c . μ f =          μ s − (μ s −μ d ) |ds| dc , |d s |≤d c , μ d , |d s |>d c . (2.22) The critical slip distance exerts control over the spatial discretization size of the model. This control must be honored to ensure the convergence of a numerical 18 simulation with the discretized fault to its continuum limit. The control should be such that it honors the physical dependence of slip of an individual element of the discretized fault on the stresses of the elements in its neighborhood. As per (Rice, 1980, 1993), the control is the following: the fault mesh discretization size must be smaller than the size of the process zone (the zone of frictional weakening near the rupture front), which affects rupture nucleation, propagation and arrest. The process zone size, which we denote by h ∗ , is related to the ratio of the shear modulus and the maximum stress drop over the critical slip distance. Denoting the numerical mesh discretization size ash m , the convention commonly followed in computational earthquake mechanics models based on slip weakening is h m < h ∗ for all fault mesh elements. Now, we need an expression forh ∗ . For faults obeying the linear slip-weakening model,h ∗ can be expressed as (Uenishi & Rice, 2003; Ma & Elbanna, 2015), h ∗ = Gd c (μ s −μ d )σ 0 n,max . (2.23) whereσ 0 n,max is the maximum effective normal traction over all the elements of the fault mesh. Taking the maximum value is necessary to ensure h m <h ∗ for all the mesh elements. Denoting τ s = μ s σ 0 n and τ d = μ d σ 0 n , τ s −τ d becomes the drop in faultstrengthand (τ s −τ d )/d c becomestheslipweakeningrate, whichisdenotedby W in (Uenishi & Rice, 2003). Another estimate forh ∗ was derived (Rice, 1980) as (9π/16)(1−ν)Gd c /(τ s −τ d ), which isO(1) equivalent to the estimate in Eq. (2.23). Inourstudy, weassumethatthefaultpermeabilitydoesnotchangeduetofault slip. We believe this assumption is justified given our objective of discriminating between production- and injection-induced fault slip based on how stresses evolve before the seismicity. Moreover, recent studies (Cappa & Rutqvist, 2011a,b; Maz- zoldi et al., 2012; Rutqvist et al., 2013) suggest that the shear-induced change in 19 fault permeability has a negligible impact on the size and the moment magnitude of the slip event. Poroplasticity Assuming an isotropic response and using τ 0 =F e S 0 F e,T , the hyperelastic constitutive relation can be written in terms ofb e and eitherτ 0 orS 0 as τ 0 = 2 ∂ψ ∂b e b e , S 0 = 2F e,−1 ∂ψ ∂b e F e (2.24) Tocompletelyspecifyaporo-elastoplasticproblem, threecomponentsareneeded–a yieldcondition, flowrulesbasedonaplasticpotential, andahardeningorsoftening rule. The yield condition specifies the state of stress at which plastic flow initiates, and the potential function is used to determine the magnitude and direction of plastic strain and other internal state variables such as the plastic porosity. The yield condition is generally written in the following form f p (S,p e ,P c ,h) =f p (τ 0 ,P c ,h)≤ 0, (2.25) The potential function for drained plastic flow is given in the form g p (S,p e ,P c ,h) =g p (τ 0 ,P c ,h)≤ 0. (2.26) A non-associative model (f p 6=g p ) ensures that the direction of plastic strain rate vector is not perpendicular to the yield surface, and energy dissipation from plas- ticity is always positive. Non-associated models can address the issue of unreason- able amounts of dilatation that are sometimes predicted by the associated model (f p = g p ), which is based on the assumption of maximum plastic dissipation. In multiphase geomechanics, non-associativity often arises from the commonly-used models for the capillary pressure-saturation relation. 20 2.2.1 Flow rules To determine the amount and direction of plastic deformation (plastic flow rule) in large displacement/rotation–large strain, first we look at the material time derivative ofb e : . b e =lb e +b e l +L v b e (2.27) wherel =l e +l p withl e = . F e F e,−1 andl p =F e L p F e,−1 . The Lie derivative of the strain is L v b e =F dC p,−1 dt F T =F d dt h F −1 b e F −T i F T =−2F e sym[L p ]F e,T . (2.28) Maximization of the dissipationD resulting from Eq. (2.9) provides us with the flowrulesfortheplasticvolumetricstrainandplasticmasscontent(Armero,1999). The former is − 1 2 (L v b e )b e,−1 = . λ ∂g p ∂τ 0 (2.29) where . λ is the plastic multiplier and ∂g p /∂τ 0 is the gradient tensor of the plastic potential function with respect to the effective Kirchhoff stress. The incremental plastic volumetric strain can be computed as . ε p v = . J p J p = . F p :F p,-T = . F p F p,−1 : 1 = tr[L p ] (2.30) The elastoplastic tangent modulus c ep can be obtained from . τ 0 = c : (d− . λ∂g p /∂τ 0 ) =c ep :d. Associative plasticity models guarantee that the elastoplastic tangent tensors are symmetric. 21 Flow rule for plastic fluid mass increments can also be obtained from the dis- sipation maximization principle. For the water phase, we obtain . M p w =ρ 0 w . λ ∂g p ∂p w =−bs w ρ 0 w . λ ∂g p ∂τ 0 : 1 ! (2.31) where it is the assumed that the fluid compressibility over a time increment is small. Substituting sequentially from Eq. (2.29), Eq. (2.28), and Eq. (2.30), we obtain . M p w =bs w ρ 0 w 1 2 (L v b e )b e,−1 : 1 =bs w ρ 0 w . ε p v , . ε p v = . λ·tr " ∂g p ∂τ 0 # (2.32) Flow rule for . s w can be written in terms of the plastic multiplier and the gradient of the yield function with respect toP c , . s w =− . λ(∂f p /∂P c ), which can be written as . s w = (∂s w /∂P c ) . P c . Another candidate formulation is . s w = ∂s w ∂J . J + ∂s w ∂P c   . P c + . J J P c   , (2.33) which makes the capillary pressure a function of deformation: ε v and φ. As mentioned above, we use a van Genuchten-type relation that qualifies as a non- associative flow rule for . s w . However, following (Coussy, 1995), we assume that the effect of deformation on the capillary relation is negligible by neglecting the . J terms. This makes our relative permeability relations, k rw and k ro vs. s w , also independent of deformation. The evolution of plastic multiplier is subject to the Kuhn-Tucker load- ing/unloading conditions: . λ≥ 0, f p ≤ 0, . λf p = 0, (2.34) 22 and the consistency condition, . f p (S,p e ,P c ,h) = 0. (2.35) based on which we obtain . λ = (∂f p /∂τ 0 ) :c :l (∂f p /∂τ 0 ) :c : (∂g p /∂τ 0 ) . (2.36) For large displacement/rotation–small strain, we use the plastic flow rule written in terms ofS 0 : . ε p = . λ ∂g p ∂S 0 , (2.37) and the plastic volumetric strain is calculated from its trace. The elastoplastic tangent modulus C ep is defined as . S 0 =C : ( . E− . ε p ) =C ep : . E. Finally, the hardening rule can be written as . h = . λ∂g p /∂ξ. The hardening rule describes how the yield condition and flow rule are modified during plastic flow. For no hardening and perfectly plastic material, the yield condition and flow rule remain constant during plastic flow. 2.2.2 Drucker-Prager model We consider the Drucker-Prager (D-P) elastoplastic model (Benallal et al., 2008) to describe the behavior of reservoir and caprock under flow-induced load- ing. The D-P yield surface, which is a smooth approximation of the M-C yield surface for continuum modeling of yielding in 3D, has been used in the literature 23 to model plastic deformation of soils and rocks. For large displacement/rotation– small strain, the model can be written in terms of S 0 andE (Bathe, 1995). We define a yield condition for perfect plasticity and with capillarity as follows f p (S,P c ,p e ) =α f (S v +βP e ) + q J 2 −γ−ζP c ≤ 0, (2.38) and the plastic potential is given by g p (S,p e ) =α g (S v +βP e ) + q J 2 , (2.39) where α f , α g , β, γ, and ζ are model constants. Here, β accounts for the plastic volumetric strain of the solid matrix (φ 0 ≤ β≤ 1) and can be assumed equal to the Biot coefficient (Coussy, 1995). A ζ > 0 value can extend the yield surface to the tensile region and allow tensile failure. To calculate increments in plastic strain and plastic fluid mass, we need gradients of functions f p andg p . Note that df p = ∂f p ∂S :dS + ∂f p ∂p e dp e . (2.40) Substituting Eq. (2.39) in Eq. (2.37), we obtain the rates of increment in plastic strain and fluid mass of slightly compressible fluids, . ε p = . λ α g 3 1 + 1 2 √ J 2 s ! , . M p = . λα g βρ 0 f (2.41) 2.3 Numerical model Below we present the numerical solution model for the problem described in the previous section. First, we present the weak forms of the mechanics and flow 24 sub-problems while accounting for finite strain plasticity and fault slip. Next, we linearize the equations that need to be solved to simulate the dynamics of fluid flow, finite strain geomechanics, and fault slip processes in the domain. Then we describe our sequential iterative solution strategy which entails solving the flow problem in a Newton loop while freezing the variations in the total volumetric stress and the volumetric plastic strain, followed by solving the mechanics problem in another Newton loop while freezing the fluid pressure (and saturations in case of multiphase flow). 2.3.1 Weak Forms The weak form of Eq. (2.6) in Lagrangian coordinates is obtained by following the standard variational recipe: take the dot product of the strong form with a test functionη belonging to a suitable space of vector functions (zero for the material points on the prescribed displacement boundary of the body), integrate over the reference volume, use integration by parts, and apply the divergence theorem. The result is (Simo & Hughes, 1998), Z V 0 P :∇ X ηdV 0 = Z Γ 0 σ T·ηdΓ 0 σ + Z V 0 ρ 0 b g·ηdV 0 (2.42) where Γ 0 σ is the Neumann (traction) boundary in the reference configuration,T = P ˆ n is the prescribed traction vector on that boundary, and ˆ n is the outward unit normal to the boundary. Prescribed displacement and prescribed traction boundaries do not overlap and their union equals the complete boundary. For time-dependent tectonic boundary conditions, Eq. (2.42) can be written in the rate form with . P = ( . FS + . SF ) in place ofP. Since we consider constant boundary conditions and use a sequential iterative solution scheme where pressure is frozen 25 during mechanics solve (details below), the weak form in Eq. (2.42) is sufficient. SinceS is symmetric, we haveP :∇ X η =S : sym[∇ X η·F T ]. Defining the virtual strain asδE = sym[∇ X η·F T ], where the variationδ is taken about the current configuration, we can write the weak form as Z V 0 (S 0 −bp e C −1 ) :δE dV 0 = Z Γ 0 σ T·ηdΓ 0 σ + Z V 0 ρ 0 b g·ηdV 0 , (2.43) The left hand side above is the internal virtual work and the right hand side is the external virtual work. The traction continuity equation across the fault,T + f =T − f , contributes two surface integral terms to the internal virtual work, corresponding to tractions on Γ + and Γ − sides of the fault (Jha & Juanes, 2014): Z V 0 \Γ f (S 0 −bp e C −1 ) :δE dV 0 + Z Γ 0 + η T (l f −bp f F −T ˆ n f )dΓ 0 (2.44) − Z Γ 0 − η T (l f −bp f F −T ˆ n f )dΓ 0 − Z Γ 0 σ T·ηdΓ 0 σ − Z V 0 ρ 0 b g·ηdV 0 =R u = 0, where p f is the fault pressure related to the pressures on the two sides of the fault (Jha & Juanes, 2014), e.g. footwall and hanging wall block pressures in case ofadip-slipfault, andR u isthemechanicsresidual. WeusetheLagrangemultiplier method to impose contact forces at the fault that are required to honor the stick- slip condition. As a result, the Lagrange multipliers acquire the physical meaning of the fault tractionsl f defined in Eq. (2.18). Also, they become unknowns of the problem, similar to the displacements. We use piecewise linear displacement test functions with nodal displacement degrees of freedom (d.o.f.) U and piecewise linear fault traction test functions with traction d.o.f. L f on the Lagrange (fault) nodes to discretize the mechanics weak form. In this study, we use a trilinear hexahedral finite element mesh with bilinear quadrilaterals on the 2D fault surface. 26 The residual is evaluated at each integration point in each element. The weak form of the fault slip equation, Eq. (2.8), is Z Γ + η·u + dΓ− Z Γ − η·u − dΓ− Z Γ f η·d f dΓ =R l = 0, (2.45) where R l is the fault residual evaluated at the Lagrange (fault) nodes. The fault slip is discretized using the same test functions as the fault traction. We denote the slip d.o.f. byD f . The weak form of the multiphase flow problem is obtained by applying the same variational recipe to water and oil mass conservation equations. For water, Eq. (2.5) gives Z V 0 . M w η p dV 0 − Z Γ 0 w Wη p dΓ 0 w − Z V 0 Jw w ·∇η p dV 0 + Z V 0 Jq w η p dV 0 =R w = 0, (2.46) where η p is the pressure or saturation test function (zero on the fixed pressure boundary), W = JF −1 w w · ˆ n is the prescribed normal mass flux on the flux boundary Γ 0 w in the reference configuration, and R w denotes the residual for the water equation. The internal flux term is expressed in terms of the product of the fluid phase mobility kk rw /μ w and the phase pressure gradient p w using Darcy’s relation. Similarly, we have the weak form of the oil phase, which we write in the residual form as R o = 0. We use the finite volume method to discretize the weak forms of the flow problem using piecewise constant η p (∇η p = 0) and element-centeredd.o.f. fortheoilphasepressureP o andthewaterphasesaturation S w . P w is obtained from the capillary pressure relation, and S o is obtained from the saturation constraint. We use backward Euler time integration method to discretize the time derivative of fluid mass. We evaluateR o andR w for each mesh element. 27 2.3.2 Mechanics model We use the Total Lagrangian method (Bathe, 1995) to linearize the mechanical equilibriumequations. Denotingthebodyconfigurationatt n+1 asx n+1 =x n +du, the incremental decomposition relations are, S 0,n+1 =S 0,n +dS 0 , E n+1 =E n +dE, u n+1 =u n +du (2.47) The deformation gradient at t n+1 can be written as F n+1 = ∂x n+1 ∂X = ∂x n ∂X + ∂du ∂x n ∂x n ∂X = [1 +∇ x ndu]F n =f n+1 F n (2.48) wheref n+1 = 1+ ∂du/∂x n is known as the relative deformation gradient (Simo &Hughes,1998;Borja,2006). Fromthedefinitionoftheincrementalstrain(Bathe, 1995), dE = sym h ∇ X du + (∇ X u n )(du∇ X ) i + 1 2 (∇ X du)(du∇ X ) = dE L +dE NL (2.49) wheredE L is the linear part of the incremental strain (linear indu) anddE NL is the nonlinear part. Noting that δE n = 0, the integrand in the stiffness term on the left hand side of Eq. (2.44) can be linearized as follows S 0,n+1 :δE n+1 = dS 0 :δdE +S 0,n :δdE NL +S 0,n :δdE L (2.50) ≈ ∂S 0,n ∂E n :dE L ! :δdE L +S 0,n :δdE NL +S 0,n :δdE L 28 wheresecondandhigherorderdutermshavebeendroppedfromtheTaylorexpan- sion of the first term, the second term is already linear becauseδdE NL is linear in du, and the third term is independent ofdu. Linearization of the mechanics weak form gave rise to the algorithmic tangent operator, ∂S 0,n /∂E n , which we denote by D n . It replaces the constitutive elastoplastic tensor (c ep in large deformation– large strain or C ep in large deformation–small strain) and ensures that Newton’s method for the mechanics solver can converge quadratically when close to the solution (Simo & Hughes, 1998; Borja, 2013). In Appendix, we provide the expres- sion for D for our Drucker-Prager elastoplastic material. Discretization results in dE L =B n L dU anddE NL =B NL dU, whereB L andB NL are linear and nonlin- ear strain-displacement matrices. The linear strain-displacement matrix contains products of the shape function derivative and u n and hence has a superscript n. Strain variations are discretized similarly. The discrete linear system of the mechanics problem now has two types of unknowns: the displacement d.o.f. U at the displacement nodes and the fault traction d.o.f. L f at the fault nodes. We solve the system using Newton’s method. Eqs. (2.47) and (2.50) are written with (k) and (k + 1) substituting for n and n + 1, i.e. [U,L f ] (k+1) = [U,L f ] (k) + [dU,dL f ], where (k) is the Newton iteration counter at t n+1 and the correction vector is the solution to the system of linear equations:    K C T l C l 0    (k)   dU dL f    =−    R u R l    (k) , (2.51) where the stiffness matrixK and the direction cosine matrixC l , which converts from the global coordinate system to the fault coordinate system, are obtained via element-by-element assembly of the nodal contributions from the displacement 29 nodes and the fault nodes. R u = [R u ] andR l = [R l ] are vectors of residuals. The element stiffness matrix for element i is obtained as follows: K (k) i = ∂R (k) u,i ∂U (k) i = Z V 0 B (k),T L,i D (k) i B (k) L,i +B (k),T NL,i S 0,(k) i B (k) NL,i dV 0 (2.52) which has a linear part with the elastoplastic tensor and a nonlinear part with the stress tensorS 0,(k) i . The third term in Eq. (2.50) becomesS 0,(k) :δE (k) , which is known and does not contribute to the stiffness matrix but appears as B (k) L,i S 0,(k) i in R (k) u,i . The pore pressure term in Eq. (2.44) does not contribute to the stiff- ness matrix because we freeze the pressure during the mechanics solve as per the sequential iterative solution scheme mentioned earlier. In Newton’s method, an error in the residual implies convergence to an incorrect solution and an error in the system Jacobian implies a slower rate of convergence. Since S 0 is inside R u , the stress path must be evaluated accurately during the simulation otherwise the simulation will diverge away from the true solution. Eq. (2.51) is a saddle point problem resulting from the use of the Lagrange multiplier method to solve the fault contact problem. To solve this system, we use the field split preconditioner option in PETSc (Balay et al., 2019; Aagaard et al., 2013), which allows us to use two different preconditioners for the elasticity and the fault sub-blocks. We use an alge- braic multigrid preconditioner for the former and an approximation of the Schur complement of K for the latter (Aagaard et al., 2013). The Schur complement is−C l K −1 C T l and it is approximated with a sum of two terms corresponding to displacement nodes on the positive and negative sides of the fault (Aagaard et al., 2013). Since only fault side nodes are involved, the Schur approximation is fast to compute. 30 Poroplastic model Large displacement/rotation–large strain First, we describe the model for large displacement, large rotations, and large strain, which may be more applicable to reservoirs with small aspect ratios (similar dimensions in different principal stressdirections)andbrittlerocks. TheleftCauchy-Greentensoratt n iscalculated asb e,n =F e,n F e,T,n . Then Eq. (2.24) is used to evaluateτ 0,n . An elastic predictor– plastic corrector scheme, which is a type of return-mapping algorithm (Armero, 1999; Borja, 2006), is used to solve the elastoplastic problem and determineb e,n+1 , τ 0,n+1 , andS 0,n+1 . The predictor step gives the trial state as b e,tr =f n+1 b e,n f n+1,T , τ 0,tr = 2b e,tr (∂ψ/∂b e,tr ), M p,tr =M p,n . (2.53) whereP c is assumed constant because we freeze pore pressures during the mechan- ics solve. If the yield condition is satisfied with the trial state, i.e. f tr p < 0, we have zero plastic flow during this time step andb e,n+1 =b e,tr , M p,n+1 =M p,tr . If the condition is not satisfied, we have plastic flow and the corrector step calculates the correction to be added to b e,tr . This correction, as per Eq. (2.27), is the Lie derivative ofb e,n+1 : L v b e,n+1 =−2 . λ 3 X i=1 ∂g ∂τ 0 i θ e,n+1 i m i , wherem i =n i ⊗n i , (2.54) θ e,n+1 i are the three principal values (eigenvalues) ofb e,n+1 obtained from the spec- tral decomposition: b e,n+1 = P i θ e,n+1 i m i ,τ 0 i are the principal values ofτ 0,n+1 , and n i are the principal directions or eigenvectors (same for both the stress and strain tensors). The logarithmic principal stretches i , which are physically meaningful and measurable quantities, are related to θ i as i = 1 2 logθ i . Integrating Eq. (2.54) 31 in its principal form over the time step dt with integration limits (t n , e,tr i ) and (t n+1 , e,n+1 i ) gives us the corrected values, e,n+1 i = e,tr i −dλ· ∂g p ∂τ 0 i n+1 , i = 1, 2, 3, (2.55) where dλ is calculated using the discrete consistency condition f n+1 p = 0. b e,n+1 andτ 0,n+1 are now updated. The plastic mass is updated as (Armero, 1999), M p,n+1 =M p,tr +bρ 0 f log J n+1 J e,n J e,n+1 J n ! (2.56) The updated stress is used to re-evaluate the yield condition, which determines whether an additional iteration of the predictor-corrector scheme is necessary. Large displacement/rotation–small strain Next, we describe the model for large displacements/large rotations, but small strains (Bathe, 1995). This may be more appropriate for laterally extensive (large aspect ratio) reservoirs and some ductile rocks deforming under the well operation. Here, the small strain feature can be exploited to replace the iterative predictor–corrector scheme with analytical expressions for the plastic multiplier and the elastoplastic stress, which change within the Newton loop of the mechanics solver. This provides computational efficiency. From Eq. (2.41), the volumetric and deviatoric parts of the plastic strain tensor are dε p v =α g dλ, de p =dλ 1 2 q J n+1 2 s n+1 , (2.57) and Eq. (2.41) provides the incremental fluid mass due to plastic strain as dM p =ρ 0 f α g βdλ (2.58) 32 Given the plastic yielding at t n+1 , we need to update the elastoplastic stress. First, we obtain an explicit algebraic expression for the plastic multiplier. Using the volumetric-deviatoric stress split from Eq. (2.17) and the incremental form in Eq. (2.47), we write: s n+1 = 2G (dE d −de p ) +s n S n+1 v =K dr (dE v −dε p v )−bdP e +S n v (2.59) Substituting Eq. (2.57) into Eq. (2.59),   1 + Gdλ q J n+1 2   s n+1 = 2GdE d +s n , (2.60) S n+1 v =K dr (dE v −α g dλ)−b P n+1 e −P n e +S n v . (2.61) Taking the scalar product of both sides of Eq. (2.60),   1 + Gdλ q J n+1 2   2 s n+1 :s n+1 = 4G 2 dE d :dE d + 4Gs n :dE d +s n :s n (2.62) Using the definition of J n+1 2 and J n 2 ,   dλ + q J n+1 2 G   2 = 2dE d :dE d + 2 G s n :dE d + J n 2 G 2 . (2.63) Defining the right hand side of the above equation as (r n+1 ) 2 , we obtain the plastic multiplier as dλ =r n+1 − q J n+1 2 G (2.64) 33 SubstitutingS n+1 v from Eq. (2.61) andJ n+1 2 from Eq. (2.64) in the discrete consis- tency conditionf n+1 p = 0 provides us the final expression for the plastic multiplier dλ = α f (K dr dE v +P n+1 e (β−b) +S n v +bP n e ) +r n+1 G−γ K dr α f α g +G . (2.65) Assumingβ =b means the role of pore pressure is limited to knowing the effective stress atnth time step,S n v +bP n e . Once we have the plastic multiplier, we calculate the stress and strain as follows. Substituting J n+1 2 from Eq. (2.64) into Eq. (2.57) results in s n+1 =de p 2G(r n+1 −dλ) dλ , (2.66) which when substituted into Eq. (2.60) gives us the expression for the deviatoric plastic strain increment de p = dλ r n+1 de + s n 2G (2.67) dε p v can be updated using Eq. (2.57) and S v can be updated using Eq. (2.61) (recall that P n+1 e is available from the solution of the flow problem at the current sequential iteration). For large positive values of S v , e.g. when the rock is under tension such as during hydraulic fracturing,dλ can be too large making it difficult to project the stress back to the yield surface. In fact, ifdλ>r n+1 , then Eq. (2.64) implies q J n+1 2 < 0, which is not physical. We address this issue by calculating dλ as the minimum ofr n+1 (equivalent to assumingJ n+1 2 = 0) and the value obtained by Eq. (2.65). 34 Fault slip model The fault tractionsL f are additionally constrained by the M-C slip condition, Eq. (2.21). When the condition is violated at a fault node, the fault slip D f is incremented. This is achieved by solving the mechanics problem in a loop: 1. The system of equations Eq. (2.51) is solved, which accounts for the nonlin- earities due to finite strain and plasticity, as described above. 2. The fault tractionsL (k) f are used to check the slip condition Eq. (2.21) and perturbation to the tractions are calculated for the fault nodes violating the condition. 3. Two local equilibrium equations, extracted from the top row of Eq. (2.51), are solved at positive and negative side fault nodes to calculate adjustment inU +,(k) andU −,(k) corresponding to the traction perturbations. The local equilibrium equations are smaller in size than the global system and hence solved using a direct solver. 4. Adjustments in the displacements are used in Eq. (2.8) to compute the slip increment. 5. If the slip increment is below a threshold, exit the loop. Otherwise, go back to step 1 above with the updated slip inR l . See Eqs. (76)-(78) in (Jha & Juanes, 2014) for the detailed algorithm. 2.3.3 Flow model The global mass conservation statements for oil and water (Eq. (2.46)) phases are written for each element in the mesh. The surface integral term with outward 35 normal flux W becomes a summation of outward normal fluxes over all the faces of an element. The discrete fluxW ij for faceij between elementi and elementj is approximated using a nonlinear two-point flux approximation (TPFA) (Nikitin et al., 2014), which expresses the flux in terms of the element pressuresP i andP j , the face transmissibility, and the fluid phase mobility ρ N k rN /μ N . The face transmis- sibility depends on the permeabilitiesk i andk j , element sizes, angles between the face normal vector and the vector connecting the element barycenters as well as on the pressures of the neighboring elements. The dependence on pressure renders the flux approximation scheme nonlinear and applicable to the full tensor permeability field that arises during a finite strain simulation, even when the initial permeability field is isotropic. The flow problem is nonlinear due to pressure-dependent fluid compressibility and viscosity and, in case of multiphase flow, saturation-dependent relative perme- ability and capillary pressure relations. We use Newton’s method to linearize the flowproblemintermsofthevectorofelement-centeredpressureandandsaturation increments. For two-phase oil-water system, denoting the oil pressure increment as dP o and the water saturation increment as dS w , the linear system is    ∂Ro ∂Po ∂Ro ∂Sw ∂Rw ∂Po ∂Rw ∂Sw    (k)   dP o dS w    (k) =−    R o R w    (k) , (2.68) where k is the Newton iteration counter. See (Jha & Juanes, 2014) for details. Nonlinear TPFA gives rise to pressure-dependent transmissibilities, as mentioned above. We assume the transmissibilities to be frozen within a time step, which introduces sparsity in the Jacobian matrix of Eq. 2.68 at the cost of a small reduc- tion in the convergence rate. Following the sequential iterative solution scheme mentioned earlier, we freeze the incremental total volumetric stress (dS v ) and the 36 incremental volumetric plastic strain (dε p v ) to their respective values at the pre- vious sequential iteration. Therefore, their contributions to the Jacobian matrix blocks is zero. 2.3.4 Sequential iterative solution strategy Given the solutions at time stepn, we use a sequential iterative scheme to solve the coupled flow-geomechanics problem at time stepn+1. This approach, which is similartoourapproachintheinfinitesimalstrainmodel(Jha&Juanes,2014), isas follows. At sequential iteration s + 1 of time stepn + 1, we solve the discrete flow system for element-centered pressures and saturations (P n+1,s+1 o ,S n+1,s+1 w ) using Newton’s method, while fixing dS v = S n+1,s v −S n v and dε p v = ε p,n+1,s v −ε p,n v . We use a block GMRES solver with Block Incomplete LU (BILU) preconditioner to solve the system. Next, we solve the discrete mechanical system for nodal displace- mentd.o.f. U n+1,s+1 , element-averagedinternalvariablesε p,n+1,s+1 v andM p,n+1,s+1 , and fault-nodal traction d.o.f. L n+1,s+1 f and slip d.o.f D n+1,s+1 f while assuming that the flow solution is fixed at (P n+1,s+1 o ,S n+1,s+1 w ). The mechanics problem is strongly nonlinear. The nonlinearity due to finite strain and plasticity is addressed via Newton’s method and an iterative return mapping algorithm in case of large strains, and the nonlinearity due to fault slip is addressed via an inner loop over the mechanics problem. We use PETSc’s field split AMG preconditioner with a GMRES linear solver (Balay et al., 2019; Aagaard et al., 2012). Convergence of the entire solution is checked by comparing the norm of the solutions at s and s + 1 iterations. If not converged, the sequential loop continues for another iteration over the two solvers–flow and geomechanics–until a converged solution, (U n+1 ,L n+1 f ,D n+1 f ,P n+1 o ,S n+1 w ), is achieved (Jha & Juanes, 2014). Depending on the strength of coupling between the flow and mechanical processes and the 37 strength of nonlinearities within each problem, the number of sequential iterations required for convergence at a time step varies with time. 2.3.5 Simulator We built our in-house simulator, Pylith-GPRS. The simulator has been described in detail (Jha & Juanes, 2014) and applied to a real oilfield (Juanes et al., 2016), a gas field (Jha et al., 2015), an aquifer (Panda et al., 2018), and synthetic CO 2 sequestration cases (Castineira et al., 2015; Jagalur-Mohan et al., 2018). The simulator solves the coupled system of nonlinear partial differential equations that model the physics of quasi-static mechanical equilibrium of the porous skeleton and conservation of mass for each fluid component (oil, water, and gas). This is a fully nonlinear, two-way coupled, multiphase geomechanics formula- tion (Coussy, 2004; Jha & Juanes, 2014). Fluid flow affects mechanical equilibrium through pressure-induced changes in the effective stress and fluid density-induced changes in the gravitational body force. Mechanical deformation affects fluid flow through the rate of volumetric stress, which changes the pore volume. We use an unconditionally stable scheme, the fixed stress scheme (Kim et al., 2011; Jha & Juanes, 2014; Castelletto et al., 2015), that solves the coupled system of dis- cretized equations sequentially and iteratively till convergence at each time step. A flowchart summarizing the interaction among the models is provided in Figure 2.3. 38 Seismic and geological information Geomechanical model Displacement and stress Fluid pressures and saturations Reservoir flow model Reservoir flow properties, multiphase fluid properties MRST GPRS Pylith Figure 2.3: Flow chart showing the workflow and connections between different modules. Available seismic and geological information are used to construct a geomechanical model and populate it with poromechanical properties. Reservoir flow properties (porosity and permeability) and fluid properties (fluid density, vis- cosity and compressibility) are obtained from an older uncoupled model and trans- ferred to the geomechanical model using the Matlab Reservoir Simulation Toolbox (MRST (Lie, 2016)). The flow-geomechanical simulator, PyLith-GPRS, outputs 3D fields of the effective stress tensor, the displacement vector, fluid phase pres- sures and saturations. 2.4 Benchmark Problem 2.4.1 The Terzaghi problem WeverifyournumericalsimulatorbysolvingtheTerzaghi’sproblem. Theprob- lem deals with a laterally constrained column subjected to a uniform compressive traction on the top surface (Figure 2.4a).The surfaces are set as no-flow boundaries expect the top surface. Model properties are listed as following: length of 50 m, compression of 2 MPa, Young’s modulus of 100 MPa, Poisson ratio of 0.3, Biot coefficient of 1, porosity of 0.2, permeability of 9.869× 10 −17 m 2 , fluid viscosity of 1 cp. We select the point on top surface to plot pressure and displacement versus time (Figure 2.4b, c). Our numerical simulation results agrees with the analytical solution. Using the same model, we then compare the results from finite strain and infinitesimal strain simulation under the cases of small and large deformation. For small deformation case, we set Young’s modulus as 100 MPa, while we set Young’s modulus as 25 MPa for large deformation case. Other properties remain the same. 39 2 a 0 100 200 300 400 0 0.2 0.4 0.6 0.8 1 Numerical solution Analytical solution 0 100 200 300 400 0 0.1 0.2 0.3 0.4 c b Figure 2.4: Terzaghi problem. (a) Model geometry with boundary conditions. The observation point is shown in a black circle. (b) Comparison of dimension- less pressure vs. time from the numerical simulation and analytical solution. (c) Comparison of displacement vs. time from the numerical simulation and analytical solution. 0 100 200 300 400 10 -1 10 0 0 100 200 300 400 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 a b 0 0.02 0.04 0.06 0.08 0.1 -10 -5 0 5 10 -3 Finite strain simulation Infinitesimal strain simulation Figure 2.5: Comparison of finite strain simulation and infinitesimal strain simula- tion under the cases of small and large deformation. (a) Displacement vs. time (b) Dimensionless pressure vs. time. As shown in Figure 2.5, as Young’s modulus decreases, displacement is higher and pressure drop is less. Because large deformation assumption, the displacement in finite strain simulation is higher than in infinitesimal strain simulation. Under the same compression, pressure drop is less for finite strain simulation. Also, the differences between finite strain and infinitesimal strain simulation is larger when we have larger deformation. 40 2.4.2 The Mandel problem We solve Mandel’s biaxial consolidation problem to test the accuracy of the mechanics-to-flow coupling in our coupled geomechanics simulator. See the model in Figure 2.6a. A uniform and constant compression of 2 MPa is applied on the left boundary of the model. The right and bottom boundaries are fixed and no- flow, while the top boundary is a traction-free drained boundary with the excess pressure set atp=0. We choose the following model properties: model dimensions L x ×L y = 50× 10 m 2 , compression 2 MPa, drained Young’s modulus E = 100 MPa, drained Poisson ratio ν = 0.3, Biot coefficient b = 0.9, initial porosity φ(t = 0) = 0.2, initial permeabilityk = 9.869× 10 −14 m 2 , and fluid viscosityμ = 1 cp. We can express these moduli in terms of the moduli used in the numerical model above: E = 9K dr G/(3K dr + G) and ν = (3K dr − 2G)/(2(3K dr + G)). We select the point near the right bottom corner as our observation point to investigate the evolution of pressure and displacement with time (Figure 2.6b, c). The primary difference between Mandel and Terzaghi problems is that the flow problem is coupled to the mechanics problem in the former and uncoupled in the latter. The mechanics-to-flow coupling of the Mandel problem leads to the well- known Mandel-Cryer effect (Wang, 2000) where the pore pressure (Figure 2.6b) and the vertical displacement (Figure 2.6c), away from the drained boundary, rises at early times, before continuing to decline due to drainage. To understand the effect of plasticity, we conduct two simulations with two different material models: a Drucker-Prager (DP) material and an elastic (EL) material. Material properties are as follows: drained Young’s modulus E = 100 MPa,drainedPoisson’sratioν = 0.3,Biot’scoefficientb = 0.9,α f andα g = 0.0562, γ = 1.5 MPa. We plot the evolution of vertical stress, strain, and stress invariants in Figure 2.7, and the evolution of vertical displacement and pressure in Figure 2.8. 41 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 1.2 Numerical solution Analytical solution 0 0.02 0.04 0.06 0.08 0.1 -10 -8 -6 -4 -2 0 2 10 -3 a b c p= 0 No- flow = 2 MPa σ 0 x y No- flow No- flow A A’ Figure 2.6: The Mandel problem. (a) Model geometry with boundary conditions. The observation point A is shown in a black circle. The observation line A’A is shown in dash line. (b) Comparison of dimensionless pressure vs. time from the numerical simulation and analytical solution. (c) Comparison of displacement vs. time from the numerical simulation and analytical solution. During early times, both the total compressive stress and the pore pressure rise due to the Mandel-Cryer effect. The first invariant I 0 1 is nearly constant at 4 MPa as the second invariant √ J 2 is rising from 1.15 to 1.27 MPa. Figure 2.8a and b also show that the deformation and pressure evolve almost identically during early times. During intermediate times, once the pressure depletion effect of drainage reachestheobservationpoint,I 0 1 startstodecreasefaster. √ J 2 continuestoincrease throughout the drainage. When the plastic yield criterion is satisfied (stress path touching the straight blue line), plastic flow begins and the material deforms along the yield line as long as the deviatoric stresses are large enough. Plastic strain begins and the total strain is larger in the plastic simulation compared to the elastic simulation. Plastic dilatancy leads to a higher consolidation (higher u y ) and a lower excess pressure in the plastic material compared to the elastic material (Figure 2.8), which agrees with the observations reported in literature (Armero, 1999). At later times, the stress state returns to the elastic region (below the yield 42 b 0 0.005 0.01 0.015 0.02 1.4 1.6 1.8 2 2.2 DP EL 1 2 3 4 1.15 1.2 1.25 1.3 1.35 1.4 DP EL Yield line 0 0.005 0.01 0.015 0.02 -1 -0.5 0 0.5 1 1.5 2 10 -3 a b c Figure 2.7: Time evolution of the (a) vertical stress, (b) vertical strain and (c) evolution of the stress state at the observation point A in Fig. 2.6a. Vertical plastic strain (blue dot line) starts to increase and the elastic simulation (EL, dash line) shows higher vertical stress and lower vertical strain than Drucker-Prager plastic simulation (DP, solid line) after stresses reach the yield line (blue line). 0 0.005 0.01 0.015 0.02 0.8 0.85 0.9 0.95 1 1.05 1.1 0 0.005 0.01 0.015 0.02 -1 0 1 2 3 10 -3 DP EL a b 0 5 10 -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 0 1 2 3 4 5 6 10 -4 c 0 0.005 0.01 0.015 0.02 -1 0 1 2 3 10 -3 DP EL Figure 2.8: Time evolution of (a) vertical displacement and (b) dimensionless pressure at the observation point A. Elastic simulation (EL, dash line) shows lower vertical displacement and higher pressure than Drucker-Prager plastic simulation (DP, solid line) after plastic failure initiates. (c) Profile of vertical strain along AA’ at t=0.02 day. Vertical strain decreases from point A to A’. Total vertical strain from DP simulation (solid line) is higher than that from EL simulation (dash line) in the cells that have plastic failure (dot line). line Figure 2.7c) because I 0 1 decreases faster due to pressure depletion than the rate of increase in √ J 2 . 2.4.3 The rigid footing problem Theproblemdealswithawater-saturatedrocklayersubjectedtoastripfooting load on the top boundary that increases with time. Because of the axisymmetry of 43 Time (day) 6 8 10 12 14 16 # sequential iterations ε p >0.01 ε p >0.02 ε p >0.04 c -0.1 0 5 4 3 2 1 0 Finite strain Infinitesimal strain E=254 MPa E=127 MPa No-flow No-flow No-flow a b x y water -0.05 0 0.05 0.1 0.15 0.2 Figure 2.9: The rigid footing problem. (a) Model geometry with boundary con- ditions. The observation point, shown as a red dot, is below the footing. (b) Comparison of finite strain and infinitesimal strain simulations for two types of rocks: weak (E = 127 MPa) and strong (E = 254 MPa). (c) Number of sequen- tial iterations for convergence between flow and mechanics problems in finite and infinitesimal strain simulations for weak rock (E = 127 MPa) the problem, we only need to consider the right half of the domain (Figure 2.9a). All the boundaries are set to no-flow boundaries except the top boundary, which allows water to drain as the layer compacts. A load F (Newton per meter) that increases in time at a constant rate is applied over the footing length L. Model propertiesareasfollows: L = 5m, compressionincreasesfrom 0withaloadingrate of 165Pa/sec,drainedPoisson’sratioν = 0.3,Biotcoefficientb = 1,initialporosity φ(t = 0) = 0.2, homogeneous isotropic initial permeability (k(t = 0) = k1) with k = 9.869× 10 −13 m 2 (1 darcy), fluid viscosity μ w = 1 cp, friction angle of 30 degree, dilation angle of 0, and γ = 0.127 MPa. We select a point below the top boundary as our observation point where we plot compression versus displacement (Figure 2.9b). Using the same model, we compare the results from finite strain and infinites- imal strain simulations for two types of rocks: weak (E = 127 MPa) and strong (E = 254 MPa), while keeping other parameters the same. Deformation is smaller in the stiffer rock, as shown in Figure 2.9b. The larger consolidation experienced by the weaker rock is correctly modeled in the finite strain simulation compared to 44 theinfinitesimalstrain, aresultthatagreeswithliterature(seeFig.10in(Carteret al., 1979)). For the same footing load, the compaction is larger in the finite strain simulation. Also, the difference between finite and infinitesimal curves increases as theconsolidationprocesscontinues. InFigure 2.9c, weshowthedifferencebetween the number of sequential iterations for convergence between flow and mechanics problems in the finite and infinitesimal strain simulations. As the plastic strain magnitude increases, both simulations need more iterations to converge. Since the plastic deformation is delayed in the infinitesimal strain model, its number of iterations lag behind that of the finite strain model and suddenly jump at t = 0.2 day. We plot the spatial profiles of pressure, displacement and plastic strain magni- tude from finite strain and infinitesimal strain simulations for the case with strong (E = 254 MPa) rock as shown in Figure 2.10. Pressure increases under com- pression except the drained top boundary where pressure remains unchanged. We observe negative displacement u y in the area under compression and smaller posi- tive displacement in other area. Plastic failure initiates on the compression bound- ary and propagates downwards. The border between the compression boundary and the flow boundary has the highest plastic strain because of high S xy in the finite strain simulation or σ xy in the infinitesimal strain simulation. This serves as a barrier and creates a failure zone. The comparison between finite strain and infinitesimal strain simulations shows that the finite strain model can capture the higher overpressure, larger displacement and plastic strain magnitude resulting from the footing load in this case. 45 0.0e+00 5.4e-02 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 plastic_strain Magnitude 0.0e+00 5.4e-02 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 plastic_strain Magnitude a b c 0 + e 0 . 5 - 0 + e 7 . 1 6 + e 5 . 4 - 6 + e 4 - 6 + e 5 . 3 - 6 + e 3 - 6 + e 5 . 2 - 6 + e 2 - 6 + e 5 . 1 - 6 + e 1 - 0 0 0 0 5 - 0 0.04 0 + e 0 . 5 - 0 + e 7 . 1 6 + e 5 . 4 - 6 + e 4 - 6 + e 5 . 3 - 6 + e 3 - 6 + e 5 . 2 - 6 + e 2 - 6 + e 5 . 1 - 6 + e 1 - 0 0 0 0 5 - Pressure (MPa) 0.08 -0.46 -0.21 (m) 0.04 0 0.027 Plastic strain 0.054 d e f -4.6e-01 3.8e-02 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 displacement Y -4.6e-01 3.8e-02 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 displacement Y Figure 2.10: Results from finite strain (a to c) and infinitesimal strain (d to f) simulations. (a) Overpressure, (b) vertical displacement (exaggerated by a factor of 10), and (c) plastic strain magnitude at t = 0.2 day. 46 Chapter 3 Coupled Multiphase Flow and Geomechanics with Elastic materials 3.1 Introduction We know that injection or production of fluids in a reservoir can perturb its mechanical equilibrium with the surrounding rock by altering the state of stress inside and outside the reservoir (Raleigh et al., 1976; Yerkes & Castle, 1976a; Segall, 1989a). Such stress changes can be significant enough to induce slip on the faults that are favorably oriented and close to frictional failure. The effect of hydrocarbon production (Segall et al., 1994; Sanz et al., 2015; Juanes et al., 2016) or water production (Gonzalez et al., 2012; Hornbach et al., 2015a; Albano et al., 2017a) and injection (Ellsworth, 2013; Keranen et al., 2013; Gan & Frohlich, 2013; McGarr,2014;Schoenballetal.,2012;Catallietal.,2013;Dengetal.,2016;Albano et al., 2017a; Lei et al., 2017) on fault and fracture activation has attracted much attention recently. Several studies (Cappa & Rutqvist, 2011a,b; Jha & Juanes, 2014; Segall & Lu, 2015b; Juanes et al., 2016) have been conducted with simple physical models of induced seismicity consisting of a confined reservoir, permeable or impermeable faults, and injection or production wells to understand the physi- cal mechanisms behind such events. In those studies, an increase in the Coulomb 47 Failure Stress (CFS) (Reasenberg & Simpson, 1992; King et al., 1994; Harris, 1998) on a fault, as a result of production or injection-induced changes in the average reservoir pressure and the poroelastic stress, is used to quantify the induced activa- tion process. Besides CFS, the stress drop (the difference between initial stress and dynamic friction strength) has been identified as a key parameter controlling rup- ture arrest and thus the size of induced earthquakes (Galis et al., 2017; Garagash & Germanovich, 2012). However, in the field, it may be difficult to discriminate between production-induced and injection-induced fault activation based on the pressure and stress data available. More research is needed to understand how pressure and stress evolve differently for production-induced and injection-induced activation. Some of the previous studies (Chang & Segall, 2016a,b) have focused on esti- mating the induced seismicity rate of faults located in the basement below the reservoir. Most studies have been conducted in a 2D domain with a flat-layered reservoir and a 1D fault with a constant friction coefficient. The seismicity rate is computed a posteriori by substituting the pore pressure and stress changes in Dieterich’s model (Dieterich, 1994), which decouples the fault response from the dynamics of flow and mechanical deformation and also neglects stress interactions amongmultiplehypocentersthatemergesuccessivelyontherupturingfaultsurface (Heimisson, 2019) . There are several limitations arising out of these assumptions. First, the fault under question can be a through-going reservoir-intersecting fault, instead of a basement fault, and since the magnitude and direction of induced stresses (normal and shear) depend on the spatial configuration of the interacting bodies, the seismicity rate of an intersecting fault can be different from the base- ment fault. For example, since the basement fault does not intersect the reservoir, the fault surface area exposed to the change in pore pressure and induced stresses 48 is expected to be smaller than the surface area in case of an intersecting fault. Second, the fault rupture behavior in 3D can be fundamentally different from its 2D approximation (Jha & Juanes, 2014). For example, in a 2D approximation of the rupture process, the rupture can propagate only updip and downdip along the fault, whereas in a 3D domain the rupture can propagate laterally within a layer before propagating into caprock and basement above and below the reservoir. Third observation is that the use of constant friction faults in induced seismicity models can limit their real-world applicability because the friction coefficient can change rapidly during slip due to slip weakening processes (Marone, 1998). A change in friction coefficient can lead to direct and indirect changes in the shear and normal stresses acting on the fault and, therefore, in subsequent seismicity. Here, in this chapter, we conduct a coupled flow-geomechanics-faulting study thatinvestigatestheaboveissuesduringCO 2 injection-andoilproduction-induced fault reactivation. We analyze fault tractions, slip, and reservoir pressure to under- stand how a fault is reactivated in successive steps and how the reactivation pro- cess selects the hypocenter locations within and outside the reservoir. We explore the role of wells parameters (injection vs. production, gas vs. liquid) so we can understand how the parameters can be manipulated to mitigate the hazard. We extract characteristic features of induced reactivation in terms of the Coulomb fail- ure stress, the Fault Stress Ratio, and the reservoir pressure perturbation caused by injection or production. We use the characteristics to understand and summa- rize the evolution of fault stability during production and injection. We use it to explain the earlier onset of injection-induced slip compared to production-induced slip when other parameters are the same between the two cases. We present a real-world application of the methodology to the Farnsworth oilfield, where CO 2 is 49 injected and oil is produced, to resolve an ambiguity in the fault stability behavior reported in an earlier publication. 3.2 Physical model We consider a dome-shaped reservoir confined between caprock and basement and intersected by a normal fault dipping at θ=80 ◦ (Fig. 3.2) that is critically stressed and impermeable to fluid flow. Our model is inspired from a previously published model used in multiple induced seismicity studies (Cappa & Rutqvist, 2011a,b; Jha & Juanes, 2014; Vilarrasa et al., 2016). The model setup resembles the Snøhvit offshore field (Hansen et al., 2013) in Norway, where 1600 kilotons of CO 2 was injected during 2007–2012 in the 75–110 m thick Tubaen aquifer situated in a fault block that is sealed by faults on the north and south sides of the aquifer. Natural gas was also produced from this field during 2007–2013. Our model also resembles the St. Johns Dome natural CO 2 reservoir across the Arizona-New Mexico border (Figure 3.1), which has been studied to understand the leakage of CO 2 along and across faults that can slip (Miocic et al., 2019). Weconsidersingle-wellinjectionofsupercriticalCO 2 inawater-saturatedreser- voir or production of oil from an oil reservoir. The well is located in the hanging wall block of the normal fault, at the top of the reservoir dome structure. The geomechanical system is initially in equilibrium with the imposed normal faulting stress regime. We assume that production and injection-induced deformations are infinitesimal, and quasi-static equilibrium conditions hold. In other words, we do not consider inertial effects and seismic wave propagation in this study. Although 50 St. Johns Dome Springerville Volcanic Field CWF Travertines GWC today Paleo GWC Buttes Fault Reservoir W E magmatic CO 2 CO production wells 2 magma leaking CO 2 Figure 3.1: A sketch of the faults and reservoirs in the St. Johns Dome natural CO 2 reservoir located in a faulted anticline. CO 2 migrates from magma across the Coyote Wash Fault (CWF) to the reservoir and the surface, where travertines are formed. This is Figure 1(b) from (Miocic et al., 2019). Used with permission. Licensing information is available at (CC, 2019). this limits our study to production or injection-induced slow slip or creep on reser- voir faults and interseismic strain accumulation and aseismic slip on large seis- mogenic faults, the assumption is valid for many oil, gas and geothermal fields. As shown in a recent study (Wang et al., 2020), flow-induced events are mostly slow slip episodes with the peak slip velocity < 4 micron/sec. Hence, the aseismic process provides an alternate mechanism for induced seismicity (Guglielmi et al., 2015; Wei et al., 2015), and understanding aseismic slip is important for a better assessment of the induced seismicity hazard. Compared to earlier work, we focus on reservoir-intersecting faults, allow for both along-strike and along-dip slip, and capture the dynamic slip weakening process. 3.2.1 Mesh initialization and boundary conditions The computational domain is a 3D box with the dimensions of 2 km× 2 km× 2 km in x, y and z directions respectively. See Fig. 3.2. We use Trelis (CSimSoft, 2016) to create a finite element mesh with hexahedral elements that honors the structure and stratigraphy of the faulted reservoir. The finite element mesh is 51 Well Fault z x b x (km) a Well Fault Tectonic compression 0.5 2.5 y (km) 80 ◦ z (km) Vertical section with fault Reservoir 2 0 0 2 y z Fault surface c Reservoir Well 3D model y Reservoir 376 m Figure 3.2: (a) Physical model with a well injecting or producing from a reservoir (gray color) confined between caprock and basement (white color) and intersected by a fault. (b) A vertical cross section of the model showing the reservoir, well and fault. (c) The fault surface viewed along the x-axis direction. Projection of the well on the fault is shown. Position of a line where detailed analysis of flow-induced stress and slip will be carried out in the Results section below. processed to calculate the element volumes, element centroid coordinates, and transmissibilities of the element faces. The model is initialized under a normal faulting stress regime with the overburden stress (σ zz (t = 0)) being the maximum principal compression stress and the horizontal tectonic stress along thex direction (σ xx (t = 0)) being the minimum principal compression stress. The overburden stress varies linearly with depth z as σ zz (t = 0) =−ρ b gz, where ρ b is the bulk 52 density calculated as ρ b =φρ f + (1−φ)ρ s based on the values given in Table 3.1. For the oil production case, ρ b = 2114.62 kg/m 3 and for the gas injection case, ρ b = 2134 kg/m 3 . The minimum principal stress is set equal to 0.7 times the overburden to create a normal faulting stress regime. Since the top of the model is at z = 500 m, the initial stress at the top boundary is balanced by applying a compression on the top boundary that is equal in magnitude to the weight of a 500 m thick rock layer with the given density. Similarly, σ xx (t = 0) is balanced by applyingaboundarycompressionequalto0.7oftheoverburdenalongthex = 2000 m boundary. Boundaries normal to the y-axis are considered fixed displacement boundaries and the y-direction initial stress is considered equal to the overburden, i.e. σ yy (t = 0) = σ zz (t = 0). These initial and boundary stresses create a state of mechanical equilibrium that results in zero displacements and strains initially. All boundaries are no-flux boundaries. For the gas injection case, we assume that the domain is initially saturated with water, and for the oil production case, we assume the reservoir is saturated with oil at irreducible water saturation. The initial pressures are in hydrostatic equilibrium. A well produces from, or injects into, the top three layers of the reservoir at a constant flow rate. The supercritical CO 2 injection rate is 849.5 cubic meter per day or 1680 ton per day at surface conditions. The rate is chosen to induce slip in a relatively short span of three months while keeping the maximum overpressure below 0.8 times the minimum principal stress to avoid hydraulic fracturing around the well. The oil production rate is 795 cubic meter per day at surface conditions, which is representative of the production from a single fault block in medium-sized oilfields such as the Long Beach Unit of Wilmington oilfield (Yang & Eacmen, 2003). The simulation duration is chosen such that the maximum change in the average reservoir pressure in both cases is similar (approximately 9 MPa). 53 Parameters Value Drained Young’s modulus, E (GPa) 10 Drained Poisson’s ratio, ν 0.25 Solid density, ρ s (kg/m 3 ) 2260 Oil density at surface conditions, ρ o (kg/m 3 ) 806.26 Water density at surface conditions, ρ w (kg/m 3 ) 1000.0 Water compressibility at reservoir conditions, c w (Pa −1 ) 4.79E-10 Water viscosity at reservoir conditions, μ w (cp) 0.4 sCO 2 density at surface conditions, ρ g (kg/m 3 ) 1.76 Biot coefficient, α 0.8 Aquifer permeability (md) 100 Caprock permeability (md) 0.0001 Fault permeability (md) 0.0001 Initial Eulerian porosity, φ 0.1 Static coefficient of friction, μ s 0.6 Dynamic coefficient of friction, μ d 0.2 Critical slip distance, d c (m) 0.1 Fault cohesive strength (kPa) 0 Initial horizontal-to-vertical total stress ratio, r tec 0.7 Well perforation top coordinates (x,y,z) (m) (949.4, 971.3, 1376.7) Well perforation bottom coordinates (x,y,z) (m) (941.1, 972.2, 1440.6) Number of grid cells in the x-direction, N x 39 Number of grid cells in the y-direction, N z 40 Number of grid cells in the z-direction, N z 43 Average grid cell size (m) 50 Table 3.1: Input parameters. 3.2.2 Model parameters We assume a homogeneous isotropic poroelastic medium with the properties given in Table 3.1. These values are representative of a sandstone aquifer at 1.5 km depth. We use a homogeneous permeability field, so we can focus on the mechanistic differences between production and injection-induced fault slip that is independent of any case-specific permeability field. The relative permeability, compressibility and viscosity of the oil phase in the production case simulation is 54 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (psi) 0 1000 2000 3000 1 1.002 1.004 1.006 1.008 1.01 0.39 0.395 0.4 0.405 0.41 (res. vol/surf. vol) (cp) Figure 3.3: Left: Relative permeabilities of oil (k ro ) and water (k rw ) as functions of the water saturationS w . Right: The formation volume factor of oil (B o ) and the oil viscosity (μ o ) as functions of the oil phase pressure p o . The oil compressibility is defined asc o = (1/B o )dB o /dp o . These properties affect the evolution of pressure and, therefore, that of induced stress. shown in Fig. 3.3. As oil production continues,S w andp o decrease, which changes the oil relative mobility and compressibility both of which control the reservoir pressure evolution and, therefore, the induced stress. The multiphase properties (relative permeabilities and capillary pressure) of the sCO 2 -water system and the physical properties of the sCO 2 phase (compressibility and viscosity) are shown in Fig. 3.4. The fault friction values satisfy Eq. (2.23). Substituting G = 4 GPa (based on E andν values above) and the maximum effective normal stress ofσ 0 n,max = 5 MPa (based on maximum value of σ 0 n in Fig. 3.8) in Eq. (2.23), we obtain h ∗ SW = 200 m. The largest element size in our simulation model is less than 55 m. Therefore, the condition of h m < h ∗ , required for fault slip convergence, is satisfied in our study. Similar to some of the previous studies (Cappa & Rutqvist, 2012; Buijze et al., 2019), we assume that fault cohesion is zero in our simulation. Uncertainties in the cohesive strength can affect estimates of fault slip at a given site and must be reduced via lab and field-scale measurements. Also, the value of d c used in 55 (psi) p g B g 0 12 24 (psi) 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (res. vol/surf. vol) 0 1000 2000 3000 4000 5000 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.02 0.04 0.06 0.08 0.1 (cp) B g a b Figure 3.4: (a) Relative permeability of water k rw increases and relative perme- ability of sCO 2 k rg decreases with increasing water saturation S w . The capillary pressureP gw =p g −p w decreases withS w because the rock is assumed to be water- wet. (b) The sCO 2 formation volume factor B g , which is the ratio of the sCO 2 volume at reservoir conditions to the volume at surface conditions, decreases with increasing sCO 2 pressure p g , due to sCO 2 compressibility: c g = (1/B g )dB g /dp g . this study falls within the range of critical slip weakening distance reported in multiple studies (Scholz, 1988; Ide & Takeo, 1997; Olsen et al., 1997; Mikumo et al., 2003; Tinti et al., 2009; Ma & Elbanna, 2015). For quasi-static slip, d c values are expected to be smaller and can be obtained from rate-and-state analysis of friction experiments (Rubin & Ampuero, 2012). However, given our focus on analyzing the evolution of stresses and pressure prior to the onset of slip, choosing a different d c value uniformly across all the simulations will not affect the main conclusions of the study. 3.3 Results We discuss the results of the two case studies shown in Fig. 3.5. We analyze the spatial profiles and time evolution of fault pressure, stresses and slip in the two cases to identify how fault activation characteristics differ between the two cases. To understand the spatiotemporal complexity of the fault activation process, we 56 Gas injection Δp > 0 Δσ n < 0 Oil production Δσ n < 0 Δp < 0 Δσ n > 0 Δσ n > 0 Fault Reservoir Caprock Figure 3.5: Two cases of the model problem studied in this paper. An ellipsoidal reservoir is confined within a low permeability caprock and intersected by a sealing fault. Injection-inducedexpansionofthereservoirappliescompressiononthefault, i.e. Δσ n < 0, and production-induced contraction applies tension, i.e. Δσ n > 0. analyze its evolution at multiple locations on the fault: within the reservoir, above and below the reservoir, and on the reservoir boundaries. We also analyze the profiles of fault shear tractionτ, effective normal tractionσ 0 n , and slip along depth on a line on the fault that is close to the well. 3.3.1 Vertical and lateral positions of the hypocenter Production-induced contraction and injection-induced pressurization of the reservoir lead to fault slip in the two cases. The magnitude of induced shear is high along the top and bottom surfaces of the reservoir due to relatively high pressure differences across the surfaces. Slip is expected to nucleate at a point on these surfaces where induced shear and tectonic shear are acting in the same direction (Segall et al., 1994). We will discuss the induced shear first. Production causes pressure depletion, which causes volumetric contraction of the reservoir. Contraction applies downward pull on the reservoir top surface and upward pull on the reservoir bottom surface. This is production-induced shear. Since the boundary conditions in Fig. 3.2(a) are such that all points on the fault experi- ence downdip tectonic shear, both the production-induced shear and the tectonic 57 0 1 2 3 4 5 -1 5 -5 0 0.5 0.4 0.3 0.2 0.1 0 0.6 0.7 MPa Updip Slip -0.5 -0.4 -0.3 -0.2 -0.1 0 MPa cm MPa 200 m Figure 3.6: Oil production-induced changes in the effective normal traction Δσ 0 n , updip shear traction Δτ, fault pressure Δp f , and updip slip plotted over the fault surface. The plots are at t = 800 day, which is 200 days after the first event time. Slip nucleates in the top reservoir layer (asterisk denotes hypocenter) and propagates laterally along the layer. shear are acting downdip on the top surface. As a result, one of the points on the reservoir top is most likely to slip compared to all other locations on the fault. This is the reason production-induced slip nucleates at the point with a circle sym- bol shown in Fig. 3.8 on the reservoir top. By a similar logic, injection-induced slip is expected to nucleate at a point on the reservoir bottom where the reser- voir expansion-induced shear is aligned with the tectonic shear. This hypocenter location is shown with the diamond symbol in Fig. 3.8. If the stress boundary conditions change to a reverse faulting stress regime, the hypocenter is expected 58 0 2 6 4 8 6 -6 0 6 4 2 0 8 10 MPa Downdip Slip 1 3 5 7 9 MPa cm MPa -4 -2 2 4 200 m Figure 3.7: sCO 2 injection-induced changes in the effective normal traction Δσ 0 n , updip shear traction Δτ, fault pressure Δp f , and downdip slip plotted over the fault surface. The plots are at t=264 day, which is 169 days after the first event time. Slip nucleates in the bottom reservoir layer (asterisk denotes hypocenter), propagates laterally in the reservoir, then updip in the reservoir and downdip in the basement. to be on the bottom surface for production-induced seismicity case and on the top surface for injection-induced seismicity case. If the fault has an offset, e.g. if the hanging wall and footwall blocks are shifted vertically, induced stresses can concentrate on the top and/or bottom of each fault compartments (Buijze et al., 2019), which will affect the location of the hypocenter. The lateral position of the hypocenter on the high-shear surface is determined by the geologic structure (dip, depth, and reservoir thickness), reservoir properties 59 (permeability, porosity and compressibility), and well parameters (injection vs. production, boundaries of the perforation interval relative to the reservoir bound- aries, and fluid type). The dip angle and depth control the magnitude of shear traction as per Eq. (2.18). Reservoir properties control the pressure change at any point along the surface. Indeed, the injection hypocenter is located at the deepest point on the reservoir bottom surface that is closest to the well. This is the point with highest shear magnitude, and it is offset from the production hypocenter in the y direction. In case of strong permeability heterogeneity, the hypocenter may be located near the point with the maximum change in reservoir pressure where the Coulomb failure criterion is met earliest. In that case, the hypocenter may be located within the reservoir instead of being on the reservoir boundary. 3.3.2 Stress path analysis In Fig. 3.8(a) and (b), we examine the time evolution of fault stresses at six carefully selected points in the shear stress vs. effective normal stress space. See the inset figure in (b) for point locations. A curve in this space is called a stress path (Parry, 2005). Injection-induced expansion or production-induced contrac- tion of the reservoir leads to changes in the shear traction. Changes in the effective normal traction is given by the interplay between the fault pressure and total nor- mal stress (Eq. (2.18)), both of which evolve with injection or production. The six points are selected to answer two questions: how do stresses evolve at the hypocen- ters (diamond and circle symbols in Fig. 3.8) and how do stresses propagate from the hypocenter into caprock and basement (squares, triangle and cross symbols). The first question is important in order to increase the accuracy of forecasting the 60 time, magnitude and location of induced events. The second question is impor- tant to understand aftershock generation and mitigation of caprock fracturing, gas leakage, and basement faulting hazards. Injection-induced overpressure leads to an increase in σ 0 n on the fault at the hypocenter and at the base of the well completion, which are shown as diamond and cross symbols in the figure. Reservoir expansion leads to a decrease in σ 0 n on the fault at two points on the top boundary, which are denoted by circle and plus symbols. To satisfy mechanical equilibrium, compression within the reservoir interval is balanced by tension on the fault outside the reservoir interval. Hence, σ 0 n at the caprock point and the basement point increase. The magnitude is smaller at the basement point because it is farther from the well compared to the caprock point. During production, changes in the total stress dominate the fault stability at all six points. Reservoir contraction leads to tensile changes in σ 0 n on points above the reservoir (circle and plus), on a point within reservoir (cross), and on a point below the reservoir (diamond). Tensile changes within the reservoir interval is balanced by compressive changes in the caprock and basement, as shown by the triangle and square symbols. When the stress path of a point reaches the static friction failure line,μ s = 0.6, in Fig. 3.8, slip is initiated at that point as per the Mohr-Coulomb criterion. After slip, the friction coefficient decreases with the slip magnitude as per the slip- weakening law (Eq. (2.22)) until it reaches the dynamic friction line,μ d = 0.2. For production (Fig. 3.8(b)), this happens at and above the reservoir top boundary, with the first event hypocenter being the point at the reservoir top boundary that is nearest to the well (the circle symbol). For injection (Fig. 3.8(a)), it happens at and below the reservoir bottom boundary, with the first event hypocenter being 61 the point at the bottom boundary with the highest shear, which is at the deepest location on the boundary (diamond symbol in the inset figure). Analyzing the stress paths of points other than the first event hypocenter reveals mechanisms responsible for stress migration and aftershock generation. During production (Fig. 3.8(b)), the point on reservoir top (plus symbol) that is deeper and has a higher compression than the first event hypocenter becomes an aftershock location. The points below the production interval, at the reservoir bottom, and below the reservoir (cross, diamond, and square) move away from failure, i.e they arrest the growth of slip. During injection (Fig. 3.8(a)), the points on and above the reservoir top boundary (plus, circle, and triangle) arrest slip because they are being pushed updip, opposite to the tectonic downdip shear. Interestingly, the stress paths of circle and plus symbols change their direction (starts moving down) when the hypocenter slips. This happens because the slip at the hypocenter applies downdip pull on the reservoir top. The point below the reservoir is an aftershock location. In contrast to production, we see that a point inside the reservoir (cross symbol) also becomes an aftershock location. This is because the injection interval is limited to only the top three layers of the reser- voir, and the bottom of the third layer is experiencing downdip shear. The implies that, for aftershock prediction in field, we need to ascertain the true (and dynam- ically changing) boundaries of the injection/production interval in a well through well logging (flowmeter and temperature surveys). It is not sufficient to know the stratigraphic reservoir boundaries, which are static in time. The first event time, or t s1 , in production and injection cases are 600 day and 95 day, respectively. Why the injection-induced event occurs at an earlier time compared to production is one of the questions we aim to answer through characteristic curves of induced fault slip (Section 3.3.5). As we will see below, 62 it is related to how fast the fault is destabilized due to τ and σ 0 n evolving with time. The rate of fault slip depends on the size and aspect ratio of the asperity and on the background stress (Galis et al., 2019). The slip rate in our production and injection cases are of the order of 10 −9 m/s and 10 −8 m/s, respectively, which are low and in agreement with the fluid-induced slow slip events observed in a recent study (Wang et al., 2020). The assumption of instantaneous mechanical equilibrium by the quasi-static mechanics simulator and the use of a relatively large value ofd c in our model lead to a large nucleation size and promote slow slip events in both production and injection cases. 3.3.3 Production-induced fault slip Production-induced contraction of the reservoir leads to a decrease in the effec- tive compression on the fault and increase in the magnitude of dip shear, which lead to fault activation. The reservoir layers shallower than the producing layers are pulled downdip, and the reservoir layers deeper than the producing layers are pulled updip (Fig. 3.9). Fig. 3.10 shows the depth profiles of fault slip and change in tractions along the fault at three time steps–immediately before failure, 100 days after failure, and 200 days after failure. We showed the spatial distribution of pressure, tractions and slip on the fault in Fig. 3.6. Production around a fault is usually associated with mechanical stabilization of the fault when the fault expe- riences similar pressure drop on both sides, and if the across-fault permeability is high enough to remove any pressure jump. In our model, the across-fault perme- ability is low (0.0001 md) to create a confined reservoir, which results into two different pressures on the two sides of the fault. The maximum of the two pres- sures, hence the fault pressure, is on the footwall side of the fault since the hanging wall side is undergoing production. Since the fault is almost impermeable, there 63 σ n τ (MPa) (MPa) μ d = 0 .2 μ s = 0 .6 time Oil production sCO 2 injection -10 -8 -6 -4 -2 0 -4 -2 0 -10 -8 -6 -4 -2 0 -4 -2 0 Inj Hypo Well Prd Hypo Fault surface 0 2 4 0 2 4 0 2 4 6 8 0 2 4 6 8 (MPa) Downdip Slip (cm) (a) (b) (c) (d) −τ μ s = 0 .6 Updip Slip (cm) Figure 3.8: Stress paths (a,b) and shear stress drop (c,d) at six points on the fault as CO 2 injection or oil production continues. Negative σ 0 n indicates compression. Negativeτ indicates downdip shear. μ s andμ d lines are the Mohr-Coulomb failure lines corresponding to the static and dynamic coefficients of friction as per the slip-weakening model. The fault surface inset in (b) shows the locations of the points and the projection of the well. (c) and (d) show the drop in shear stress with increasing slip, which is related to the amount of seismic/aseismic energy released. is negligible change in the fault pressure over the simulation period. Therefore, the effect of pressure is not dominant in this case and the change in total normal stress, Δσ n , dominates the change in effective normal traction Δσ 0 n in Eq. (2.18). Here, the change is defined with respect to the initial value. 3.3.4 Injection-induced fault slip Initially, the reservoir is fully saturated with water and has a hydrostatic pres- sure distribution. sCO 2 is injected with a constant flow rate of 30 million standard 64 -10 -5 -2.5 0 0.88 0.8775 0.875 MPa -0.11 0.04 0 -0.04 c d meter a b -0.08 0.11 0.08 -7.5 S o Figure 3.9: 3D fields of (a) overpressure Δp, (b) oil saturation S o , (c) horizontal displacement U x , and (d) vertical displacement U z resulting from 1000 days of oil production. Well location is shown with a black line in (c) and (d). Δp is almost zero (yellow color) in the footwall block due to the low permeability of the fault. S o decreases due to oil production. Cut-away views of U x and U z fields are shown using a slice along the fault surface and a y = 1000 m slice. Production-induced contraction of the reservoir causes negative (blue color) U x displacements on the fault in the reservoir interval, negative U z in the overburden above the reservoir, and positive (red color) U z in the basement. cubic feet per day (MMSCF/day) for 264 days. Overpressure, sCO 2 accumula- tion, and deformations are shown in Fig. 3.11. Increase in the reservoir pressure decreases the effective compression on the fault and injection-induced reservoir expansion causes an increase in the magnitude of dip shear, both of which lead to fault slip. The change in shear traction induced by injection is opposite in 65 -3 03 0.5 1 1.5 2 Depth along fault (km) -3 03 -3 03 -1 -0.5 0 t = 600 day t = 700 day t = 800 day CFF Slip (cm) (MPa) (MPa) (MPa) Injection Reservoir Figure 3.10: Depth profiles of changes in tractions and slip on the fault due to oil production. direction to the one induced by production. Although injection-induced expan- sion of the reservoir applies slight compression on the fault leading to Δσ n < 0, injection-induced overpressure dominates in Eq. (2.18). As a result, the effective compression on the fault decreases, Δσ 0 n > 0, and the fault slips. During slip, the points on the fault below the hypocenter also show an increase in the shear traction (see the depth profiles in Fig. 3.12) as they move downdip due to normal faulting. We plotted the spatial distribution of pressure, tractions and slip on the fault in Fig. 3.7. 3.3.5 Characteristics of fault activation Weareinterestedinidentifyingthedifferencesbetweenproduction-inducedand injection-induced fault activation in terms of quantities that can be measured or inferred in the field. Based on the spatiotemporal behavior of fault stresses and pressure discussed in Section 3.3, we define three quantities—the Coulomb Failure 66 Stress, CFS =|τ|+μ f σ 0 n , theFaultStressRatiodefinedtobetheshearmagnitude- to-effective normal traction ratio, FSR =|τ|/σ 0 n , and the dimensionless reservoir pressure,p rd ≡p res /p res,0 , where subscript 0 indicates the initial value att=0. The reservoir pressure is defined as p res ≡ R Vres p dV/V res , where V res is the volume of the reservoir in the hanging wall block because the well is located there (Fig. 3.2). In larger scale simulation models with larger size of individual mesh elements, the reservoir pressure values may change significantly from one mesh element to the next due to jumps in permeability. In such cases, V res should correspond to the reservoir volume of a mesh element adjacent to the potential hypocenter instead of the entire reservoir. The reservoir pressure is a measured quantity in the field that decreases during net production and increases during net injection when multiple wells are injecting or producing. Measurement of stresses in the field has been a less common practice, but that is changing fast in the petroleum industry given the role of stresses on oil recovery, hydraulic fractures, horizontal well placement, ground deformation, and induced seismicity. Stresses are regularly calculated using information about optimally oriented faults, background stresses, borehole measurements, and seismic, microseismic and geodetic data (McGarr, 1999; Yaghoubi & Zeinali, 2009; Schmitt et al., 2012; Hearn & Burgmann, 2005; Catalli et al., 2013; Vasco et al., 2017). In certain cases with known rock-fluid properties, fault dip and the faulting regime, we can obtain the initial values of CFS and FSR analytically as described below. Initial CFS and FSR Since we are considering a dip-slip fault in the xz plane, we can perform an approximate analysis in the xz plane. The unit normal vector at a point on the fault isn = [n x ,n z ], where n z /n x = cotθ (normal vector pointing downward and 67 right in Fig. 3.2a) and n 2 x +n 2 z = 1. The dip angle 0 ◦ < θ < 90 ◦ is a constant in our model. Given the shear-free nature of the boundary and initial stresses, the total stress tensor at t=0 is σ 0 =    σ xx,0 0 0 σ zz,0   . (3.1) Defining r tec ≡σ xx,0 /σ zz,0 , we have r tec =0.7 in our model. Assuming an equation- of-state for a slightly-compressible fluid, the fluid density can be written as a function of the initial density, the compressibility and the pressure change asρ f = ρ f,0 exp (c f (p−p 0 )). The initial pressure is hydrostatic,p 0 =ρ f,0 gz, whereg is the magnitude of acceleration due to gravity. Assuming homogeneous porosity φ and solid density ρ s , we have σ zz,0 =−g R ρ b dz =−ρ b gz. Therefore, p 0 =− ρ f ρ b σ zz,0 . Given the definition of effective stress and using the definition, r den ≡ ρ f /ρ b , we obtain σ 0 xx,0 =r tec σ zz,0 −αr den σ zz,0 , σ 0 zz,0 =σ zz,0 −αr den σ zz,0 (3.2) We define the effective stress ratio (Sibson, 1985) as r ≡ σ 0 xx,0 /σ 0 zz,0 = (r tec − αr den )/(1−αr den ). Note 0<r tec <1 and 0<r <1 for the normal faulting tectonic regime and r tec >1 and r>1 for the reverse faulting regime. Magnitude of the effective traction vector at t=0 is given by |l 0 | 2 = σ 0 xx,0 n x 2 + σ 0 zz,0 n z 2 (3.3) and the effective normal traction is σ 0 n,0 =σ 0 xx,0 n 2 x +σ 0 zz,0 n 2 z = (rn 2 x +n 2 z )(αr den − 1)ρ b gz = (rn 2 x +n 2 z ) α− 1 r den p 0 , (3.4) 68 i.e., initially, the effective normal traction is linearly dependent on pressure. For the oil production case, r den = 0.38 and for the sCO 2 injection case, r den = 0.47. From above, we can obtain the following ratio of fault tractions, |l 0 | σ 0 n,0 ! 2 = σ 0 xx,0 n x 2 + σ 0 zz,0 n z 2 σ 0 xx,0 n 2 x +σ 0 zz,0 n 2 z 2 = r 2 + cot 2 θ n 2 x (r + cot 2 θ) 2 (3.5) From τ 2 =|l| 2 −σ 02 n , the initial fault stress ratio is FSR 0 = |τ 0 | σ 0 n,0 = (r− 1) cotθ r + cot 2 θ , (3.6) which is negative since r < 1 and almost independent of depth (r may be weekly depth-dependent through the fluid compressibility) despite τ andσ 0 n being depth- dependent. This means it has the same value for different hypocenter locations of the production and injection cases. We know, CFS 0 =|τ 0 | +μ s σ 0 n,0 . (3.7) Substituting Eq. (3.6) and using Eq. (3.4) to expressσ 0 n,0 in terms ofp 0 , we obtain CFS 0 =−p 0 (μ s +FSR 0 ) (r sin 2 θ + cos 2 θ)((1/r den )−α), (3.8) which varies with depth; it has a higher value for the injection hypocenter than the production hypocenter when the well is operating in the hanging wall block of a normal fault because the injection hypocenter is deeper. 69 Changes in CFS and FSR The value of tectonic stress ratio, r tec , may be uncertain around a fault, which means CFS 0 and FSR 0 may be uncertain. As a result, most induced seismicity studies focus on evaluating the change in Coulomb stress resulting from pressure changes associated with fluid production and injection. Therefore, we present our analysis in that framework. We define ΔCFS = CFS(t 2 )−CFS(t 1 ) = Δ|τ| +μ s Δσ 0 n ΔFSR = FSR(t 2 )−FSR(t 1 ) = |τ| σ 0 n ! 2 − |τ| σ 0 n ! 1 , (3.9) where the subscripts 1 and 2 refer to t 1 and t 2 , and we are considering the period before slip, i.e. t s1 > t 2 > t 1 , so μ f = μ s . It is common to assume t 1 = 0 so the changes are with respect to the initial value. For t > 0, it is difficult to find an explicit analytical relation between CFS and FSR except for geometrically simple cases with spatially uniform material properties and constant boundary conditions. Approximating our reservoir as an isotropic ellipsoidal inclusion within a caprock (see Figure 3.5), the total strain inside the reservoir caused by a uniform pressure change of Δp res can be decomposed into a poroelastic strainε e and a ‘stress-free’ or eigenstrainε ∗ as (Wang, 2000) ε =ε e +ε ∗ = 1 2G Δσ− ν 1 +ν Δσ kk I + α 3K dr Δp res I (3.10) where Δσ kk (summation over repeated indices) is the change in mean stress. ε is the constrained strain of the reservoir inside the caprock andε ∗ is the eigenstrain or ‘stress-free’ strain caused solely by Δp res in absence of any constraint from the caprock. AsperEshelby’sequivalentinclusionmethod(Eshelby,1957),ε =S :ε ∗ , 70 whereS istherank-4Eshelbytensor(Mura,1982), whichisrelatedtotheelasticity tensor asS =PC, whereP is the Hill polarization tensor (Hill, 1965). For an ellipsoidal reservoir,S is independent of spatial position. The stress change can be calculated as Δσ ij =C ijkl ε e kl =C ijkl (S klmn ε ∗ mn −ε ∗ kl ) (3.11) Substitutingε ∗ from Eq. (3.10), we obtain (Segall & Fitzgerald, 1998) Δσ = −αΔp res Σ, (3.12) Σ ij = δ ij − 1− 2ν 1 +ν S ijkk − ν 1 +ν S kknn δ ij (3.13) where Σ is the rank-2 tensor of reservoir geometry (aspect ratios). The changes in tractions on the reservoir-intersecting fault, at points within the reservoir interval and prior to any slip event, can be calculated as Δ|τ| = α|Δp res | sinθ cosθ(Σ xx − Σ zz ), Δσ 0 n = Δp f −αΔp res (Σ xx sin 2 θ + Σ zz cos 2 θ) (3.14) whereS xzkk = 0 has been assumed. The traction changes are with respect to the initial traction values at the given point. For our ellipsoidal-shaped reser- voir in the hanging wall block, the ratio of major and minor axes lengths is approximately 600/50 = 12, for which Σ xx ∼ 0.93(1− 2ν)/(1− ν) = 0.623 and Σ zz ∼ 0.13(1− 2ν)/(1− ν) = 0.087. During injection, Δp res > 0 and Δp f ∼ Δp res . During production, Δp res < 0 and|Δp f ||Δp res | because the fault pressure is defined from the footwall side pressure, which is higher and therefore leads to a more accurate estimate of slip onset time than the hanging wall side pressure. During production, the footwall side pressure drops slightly, 71 even for a sealing fault, because the fault is being pulled towards the hanging wall block, which causes an undrained expansion of the footwall volume adja- cent to the fault. Using the undrained deformation theory to approximate this pressure change, Δp f =−BΔσ kk /3, where B = αM/K u is Skempton’s coeffi- cient (Wang, 2000) defined in terms of the Biot modulus M and the undrained bulk modulus K u = K dr +α 2 M. Δσ kk can be approximated from Eq. (3.12) as −αΔp res (Σ xx + Σ zz ). Hence, Δp f ∼ (αB/3)(Σ xx + Σ zz )Δp res . Substituting these valuesinEq.(3.14)yields Δ|τ|linearlyproportionalto|Δp res |duringbothinjection and production, Δσ 0 n linearly proportional to Δp res during injection, Δσ 0 n weakly proportional to Δp res during production, and Δ|τ| > 0, Δσ 0 n > 0. These results are confirmed in the high-resolution simulations as shown in Figs. 3.8 and 3.13. The changes in FSR and CFS become ΔFSR = |τ 0 | + Δ|τ| σ 0 n,0 + Δσ 0 n − |τ 0 | σ 0 n,0 = Δ|τ|− Δσ 0 n FSR 0 σ 0 n (3.15) ΔCFS = Δ|τ| +μ s Δσ 0 n =σ 0 n ΔFSR + (μ s +FSR 0 )Δσ 0 n (3.16) Based on the proportionality between the fault tractions and the reservoir pressure (Eq. (3.14)), the above relations imply: (1) ΔCFS is linearly proportional to Δp res during both production and injection, and (2) ΔFSR is nonlinearly proportional to 72 Δp res during injection and linearly proportional to Δp res during production. The relation between CFS and FSR can be approximated as ΔCFS = β 1 ΔFSR +β 2 Δp res ΔFSR +β 3 Δp res (3.17) β 1 = σ 0 n,0 < 0 (3.18) β 2 = (β 4 −α(Σ xx sin 2 θ + Σ zz cos 2 θ)) (3.19) β 3 = (μ s +FSR 0 )β 2 (3.20) Δ|τ| Δσ 0 n = |Δp res | Δp res α sinθ cosθ(Σ xx − Σ zz ) β 2 , (3.21) where for points within the reservoir, β 4 = 1,β 2 > 0,β 3 > 0 during injection and β 4 = αB 3 (Σ xx + Σ zz ),β 2 < 0,β 3 < 0 during production. Given the magnitudes of coefficients, the first term on the right hand side of Eq. (3.17) dominates over the other two terms. This implies that d(ΔCFS)/d(ΔFSR) is negative and its magnitudedependsontheinitialpressurep 0 asperEq.(3.4). Eqs.(3.21)and(3.15) imply ΔFSR< 0 and ΔCFS> 0 within the reservoir. Now, we analyze the fault points outside the reservoir, e.g. in the caprock and basement. Here, the change in total stress is Δσ c xx =αΔp res 1−2ν 1−ν − Σ xx , which is opposite in sign and smaller in magnitude compared to Δσ xx in Eq. (3.12) because 1−2ν 1−ν > Σ xx . Inotherwords,thecaprockandbasementexperiencecompressiondur- ing production and tension during injection. As a result, β 2 changes in Eq. (3.19). The pressure change in caprock/basement is also different because of the perme- ability barrier. Low permeability of caprock/basement means the pressure change can be approximated using the undrained deformation theory: Δp c =−BΔσ c kk /3, which is much smaller in magnitude compared to Δp res . Therefore, the change in total normal traction determines the change in the effective normal traction. This is confirmed by Fig. 3.8: σ 0 n increases at square (basement) and triangle (caprock) 73 in Fig. 3.8a and decreases at square and triangle in Fig. 3.8b. For points on the reservoir boundary, the result can be found by combining the results within and outside the reservoir. Three types of characteristics Wedefine Δp rd ≡|1−p rd |asthedimensionlessreservoirpressureperturbation, which increases regardless of production or injection. We propose three charac- teristics of induced fault slip: ΔCFS vs. Δp rd , ΔFSR vs. Δp rd , and ΔCFS vs. ΔFSR. In Fig. 3.14, we plot the first two types at the hypocenters of produc- tion and injection cases. Before slip, ΔCFS vs. Δp rd characteristics are almost linear for both production and injection cases. However, the slopes are different because of the difference in the dependency of Δσ 0 n on Δp res (see Eq. (3.14)). The CFS characteristic is linear because the shear and effective normal tractions are linear functions of the reservoir pressure (see Fig. 3.13 and the discussion below Eq. (3.16)). When the shear stress reaches the friction strength, i.e.|τ| =−μ f σ 0 n (zero cohesion assumed), the fault slips, and we have CFS(t) = 0 for t≥t s1 . This value of CFS is maintained as the well continues to produce or inject. Hence, for t≥ t s1 , ΔCFS(t) =−CFS 0 , which has different values for production and injection cases because their hypocenter depths are different:−CFS 0 is larger for injection because the hypocenter is deeper. The hypocenter depths are sensitive to the tectonic stress regime as explained above. These observations suggest that given ΔCFS vs. Δp rd data, we may not be able to determine whether the impend- ing seismicity is associated with production or injection because the characteristic behavior is linear in both cases. Another difficulty is that ΔCFS is sensitive to errors in the hypocenter depth estimation (Catalli et al., 2013). Uncertainty in 74 the hypocenter depths (Dahm et al., 2007; Cesca et al., 2011) can create uncer- tainty in the estimated ΔCFS values and, therefore, in our forecasts of the induced seismicity hazard. On the other hand, the ΔFSR vs. Δp rd characteristic at the hypocenter (Fig. 3.14b, gray region) is approximately linear for production and nonlinear for injection before slip occurs. After slip, the FSR characteristic curve for injection shows a more rapid shear relaxation, compared to production, due to a higher CFS(t=t s1 ) at the slip onset time or the first event time. As explained below Eq. (3.14), the linearity in the production case is a result of the linearity between Δ|τ| and Δp res and the weak dependency between Δσ 0 n and Δp res . The nonlin- earity in the injection case emerges from the fact that both Δ|τ| and Δσ 0 n are strongly dependent on Δp res . These observations have two implications. (1) Given ΔFSR vs. Δp rd data, either from measurements (Vasco et al., 2017) or numerical simulations, we can visually check whether the FSR characteristic is linear or non- linear and use that knowledge to determine whether to control production wells or injection wells to avoid an impending fault slip event. This can augment the traffic light system proposed earlier (Mignan et al., 2017) to mitigate the induced seismicity hazard. It can be challenging to identify the nonlinear trend if the measured/estimated stress data is noisy or if the simulation used to generate the data is only one-way coupled, i.e. does not account for the feedback from the induced stresses into the pore pressure evolution. (2) Compared to CFS, FSR is less sensitive to errors in the hypocenter depth estimation because FSR is a ratio of tractions. This makes it a more robust estimator of the impending slip. Next, weplotthethirdcharacteristic, ΔFSRasafunctionof ΔCFSinFig.3.15, to extract additional features that might be diagnostic/predictive of impending 75 slip. The region bounded by ΔCFS > 0 and ΔFSR < 0 marks the fault desta- bilization region. We clearly see how the destabilization and stabilization trends dependonthelocationofthepointonthefault(withinoroutsidethereservoir)and the type of well operation (injection or production). There are several interesting and novel results. The characteristic curves at the injection hypocenter (diamond) and within the reservoir (cross) are nonlinear compared to the linear characteristic curve at the aftershock location (square). The analysis below Eq. (3.17) explains this in terms of the reservoir pressure evolution. p res is increasing within the reser- voir during injection, and therefored(CFS)/d(FSR) increases for the diamond and cross points, which causes their characteristics to bend upward compared to the straight line characteristic of the square symbol where the pressure is almost con- stant. In case of production (Fig. 3.15b), points within the reservoir experience large changes in p res , which explains the nonlinearity of the CFS-FSR characteris- tic at the cross point. The pressure at points on the top reservoir boundary (circle and plus) is almost constant, and the characteristic is linear. 3.3.6 The slip onset time The approximate linearity of the ΔCFS vs. Δp rd characteristic suggests a constant value of dCFS dp rd . Therefore, we can estimate the slip onset time or the first event time as t s1 = −CFS 0 dCFS dp rd dp rd dt , (3.22) which is valid for points experiencing destabilization, i.e. ΔCFS > 0. Above, CFS 0 can be estimated using Eq. (3.8) and the average reservoir pressure evolution dp rd dt , which is positive for both production and injection, is known in the field from wellhead pressure measurements or history match simulations. The pressure 76 evolution in our model is almost linear for both types of induced slip, which implies a constant value of dp rd dt . Using these values in Eq. (3.22), we obtain t s1 = 560 day for the production-induced slip case and t s1 = 87 day for the injection-induced slip case. The error in these values is less than 7% compared to the simulated failure times of 600 day and 95 day in the respective high-resolution simulations. A higher value of dCFS/dp rd during injection compared to production answers the question of why injection-induced seismicity occurs earlier than production- induced seismicity. Given the linearity of the ΔFSR vs. Δp rd characteristic for production-induced slip, it can also be used to estimate the time as t s1 =−(μ s + FSR 0 )/ d(|τ|/σ 0 n ) dp rd dp rd dt , which yields 564 day for the production case. This is within 7% of the simulated time. We do not estimate the slip time for the injection case using the FSR characteristic due to its non-linearity. 3.4 A real-world application To test our methodology on a real-world case, we apply it to the dataset of the Farnsworth oilfield in Texas (Zhao & Jha, 2019), which is a part of the Southwest Regional Partnership on Carbon Sequestration established by the US Department of Energy. The Farnsworth field is producing oil, gas and water from the Morrow-B sandstone reservoir since year 1963. Waterflooding started in 1964 to increase the oilrecovery. CO 2 injectionstartedinDecember2010aspartofthetertiaryrecovery and carbon sequestration. More detail can be found in Chapter 5. Fig. 3.16 shows the evolution of ΔFSR as a function of ΔCFS on three faults in the field. The evolutions are plotted on fault points located at the top of the reservoir. We observe that the general trend of the ΔFSR vs. ΔCFS agrees with our prediction in Fig. 3.15. There is one interesting observation in Fig. 3.16(a) that requires 77 further discussion. This relates to the gray-colored cloud of points corresponding to ΔCFS∼−0.4 MPa during year 2013-2016 (21,300-22,250 day), which is during the CO 2 injection period. Based on a relatively constant ΔCFS, one may interpret this as a period of neutral stability, i.e. CO 2 injection is neither destabilizing nor stabilizing the fault, as concluded in an earlier study (Zhao & Jha, 2019). However, in this study, we see that ΔFSR is increasing during this period, which suggests fault stabilization. We show that FSR is able to resolve the ambiguity in the CFS response regarding the fault stability behavior. The reason lies in the difference between the definitions of CFS and FSR, as explained next. Consider Eq. (3.9) with t 1 = 21, 300 day and t 2 = 22, 250 day for a point on Fault-1. During this time period, a decrease in shear (Δ|τ| < 0) can get nullified by an increase in the fault pressure, Δσ 0 n = Δσ n + Δp f > 0 resulting in ΔCFS∼ 0. For example, the shear magnitude on Fault-1 decreases by Δ|τ| =−0.05 MPa and the average reservoir pressure increases due to injection by Δp res = 0.1 MPa causing ΔCFS∼−0.05 + 0.6(Δσ n + 0.1)∼ 0, which suggests neutral stability, since Δσ n due to the compression applied on the fault by the expanding reservoir is negligibly small. For the same values, ΔFSR∼ 0.0092− 0.0077 = 0.0015 MPa, which is positive and suggests stabilization. 3.5 Conclusions To address the fault reactivation hazard in oilfields and other subsurface oper- ations where both production and injection of fluids are ongoing, it is imperative to find distinctive characteristics of production and injection-induced reactivation. Using high-resolution simulations and theories of Coulomb failure and poroelas- ticity, we identify two such characteristics–the change in Coulomb Failure Stress 78 (CFS) and the change in Fault Stress Ratio (FSR) as functions of the average pres- sure perturbation applied to the reservoir. The characteristics can be plotted at selected points along reservoir surfaces across which large pressure differences are anticipated or along boundaries experiencing large shear due to subsidence/uplift, because such points are likely hypocenters of induced seismicity. For the simplified cases of oil production or CO 2 injection in a homogeneous, isotropic, linear elas- tic reservoir intersected by a normal fault, we obtained the characteristic curves for the two types of induced activation. We demonstrated the application of the proposed method to a real-world oilfield where CO 2 injection and oil production occurs. Similar to the widespread success of "type curves" for pressure transient analy- sis and well testing in the petroleum engineering and hydrology domains (Ferris et al., 1962; Bourdet, 2002), we envision that standard type curves of CFS and FSR can be constructed for different types of reservoir-fault configuration (reservoir bounding fault, reservoir-intersecting fault, basement fault), and different types of faulting regime (normal, reverse). Characteristics drawn using field measure- ments can be compared against the standard type curves to determine the type of impending seismicity as well as to forecast the slip onset time in certain cases. The type curves can also be used in a diagnostic mode during post-mortem analysis of any seismic/aseismic fault reactivation event. That is, given the onset time and the reservoir-fault geometric configuration, quantitative comparison between the field characteristics and the standard type curves can be used to identify the wells (producers or injectors) to control. 79 0.8 0.6 0.4 0.2 0 c d meter b -0.11 0.04 0 -0.04 -0.08 0.11 0.08 0 5 7.5 10 MPa 2.5 a Figure3.11: 3Dfieldsof(a)overpressure Δp, (b)sCO 2 saturationS g , (c)horizontal displacementU x , and (d) vertical displacementU z resulting from 264 days of sCO 2 injection. Well location is shown with a black line in (c) and (d). Δp is almost zero (blue color) in the footwall block due to the low permeability of the fault. The sCO 2 saturationS g builds up near the crest of the dome due to the buoyancy and mobility of the injected sCO 2 . Cut-away views ofU x andU z fields are shown using a slice along the fault surface and ay = 1000 m slice. Injection-induced expansion of the reservoir causes positive (red color) U x displacements on the fault in the reservoir interval, positive U z in the overburden above the reservoir, and negative (blue color) U z in the basement. 80 -10 0 10 0.5 1 1.5 2 Depth along fault (km) t = 95 day t = 164 day t = 264 day -5 0 5 -5 0 5 CFF -10 -5 0 Slip (cm) (MPa) (MPa) (MPa) Injection Reservoir Figure 3.12: Depth profiles of changes in tractions and slip on the fault due to sCO 2 injection. -10 -5 05 10 -2 0 2 4 6 8 Traction (MPa) Δp re s Δσ n Δτ Δτ Δσ n Injection Production (MPa) t s1 t s1 Figure 3.13: Evolutions of the effective normal traction, σ 0 n (dash lines), and the updip shear traction,τ (solid lines), as functions of the change in average reservoir pressure, Δp res , for production and injection cases. Evolutions are shown at the respective hypocenter locations shown in Fig. 3.8 (circle for production and dia- mond for injection). t s1 marks the slip onset times for the two cases. The linear behavior in the tractions before t s1 is predicted by Eq. (3.14). 81 0 0.2 0.4 0.6 -1 0 1 2 3 0 0.2 0.4 0.6 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 Inj Inj Prod Prod ΔF SR ΔCF S t s1 t s1 (a) (b) Figure3.14: Characteristiccurvesofinducedfaultslip. (a)ThechangeinCoulomb failure stress ΔCFS and (b) the change in fault stress ratio ΔFSR as functions of the dimensionless reservoir pressure perturbation Δp rd . The evolutions are plotted at the hypocenters of production-induced (circle) and injection-induced (diamond) slip. Gray region indicates the region of induced fault destabilization. (a) In the gray region, before slip, ΔCFS at the event hypocenters evolve approximately linearly. (b) In the gray region, before slip, ΔFSR shows nonlinear behavior for injection and approximately linear behavior for production. For both injection- and production-induced slip, ΔFSR(t s1 ) = −(μ s + FSR 0 ), which is a constant value because FSR 0 is independent of the hypocenter depth. 82 -1 0123 -0.4 -0.2 0 0.2 ΔCF S ΔF SR Oil production (a) (b) sCO 2 injection -1 0123 -0.4 -0.2 0 0.2 time time ΔCF S Inj Hypo Prd Hypo t s1 t s1 time time Figure 3.15: ΔFSR vs. ΔCFS characteristic curves of induced slip for (a) super- criticalCO 2 injectionand(b)oilproduction. Thesymbols–square, triangles, circle, cross and plus–denote different locations on the fault (see the inset figure). Gray region is the region of induced fault destabilization. The curves start at (0, 0) and evolve with time along a straight line of negative slope, with deviations from the line indicating how fast the pressure is changing in the vicinity of that point. ΔCFS after failure is equal to the negative of initial CFS, which is different for points at different depths. 83 -0.4 -0.2 0 -2 0 2 4 6 8 10 10 -3 (MPa) -0.4 -0.2 0 -2 0 2 4 6 8 10 10 -3 0 0.5 1 1.5 2 2.5 3 10 4 (MPa) -0.4 -0.2 0 -2 0 2 4 6 8 10 10 -3 (MPa) CO inj CO inj CO inj water inj water inj water inj time (day) Fault-1 Fault-3L Fault-3R Fault-3L Fault-3R Fault-1 Fault-6 (a) (b) (c) (d) time Figure 3.16: (a), (b), and (c) show the co-evolution of FSR and CFS on three faults in the Farnsworth oilfield. The grayscale of the dots, with the color bar on the right, and the arrows indicate the direction of time. The evolutions are plotted on fault points at the top of the oil reservoir, which are identified as blue dots in (d). The ellipse in (a) highlights the cloud of points during 21300-22350 days when ΔCFS is relatively constant. 84 Chapter 4 Coupled Multiphase Flow and Geomechanics with Elastoplastic materials 4.1 Introduction Production and injection activities in faulted and stress-sensitive reservoirs can raise concerns of fault destabilization (Raleigh et al., 1976; Yerkes & Castle, 1976a; Segall, 1989a). Coupled numerical modeling of fluid flow and mechanical deforma- tion in such reservoirs is an active area of research for the purpose of assessing the hazards of fault reactivation and induced seismicity during oil production (Segall, 1985; Segall et al., 1994; Minkoff et al., 2003a; Chan & Zoback, 2007; Juanes et al., 2016), gas production (Sanz et al., 2015; v. Thienen-Visser & Fokker, 2017), water extraction (Gonzalez et al., 2012), wastewater disposal (Keranen et al., 2013; Hor- ton, 2012; Hosseini et al., 2018), geologic CO 2 sequestration (Morris et al., 2011; Jha & Juanes, 2014; Zhao & Jha, 2019), and hydraulic fracturing (Deng et al., 2016). However, most of these studies assume a linear elastic infinitesimal defor- mation response of the rock hosting the fault. We know that deformation response of many rocks, such as clayey sands, ultra-low permeability shales, and carbonates can become viscoelastic or elstoplastic at large effective stresses typically encoun- tered during long-term production or injection (Maury et al., 1996; Bennett et al., 85 2015; Davis & Selvadurai, 2002). Hard brittlerocks also show inelastic stress-strain behaviors through microcracking activity (Scholz, 1968). Long-term groundwater withdrawal in California’s Central Valley (Ireland et al., 1984; Galloway & Riley, 1999) and oil and water production in Wilmington (Allen, 1968; Allen & Mayuga, 1969) and Gulf of Mexico (Chan & Zoback, 2007) oilfields are examples where both large strains (e.g., > 2%) and plastic deformation are expected to occur. Plastic deformation also occur at relatively small values of the induced stresses if the rock is mechanically weak or ductile, e.g. chalk reservoirs in Ekofisk (Chin & Thomas, 1998), some shales (Sone & Zoback, 2011; Chang et al., 2014), and clay beds or aquitards (Gu & Ran, 2000). Presence of clay particles or interbeds in sandstone aquifers can give rise to a lag in the mechanical response of aquifers under deple- tion (F.G. Bell, 1986; Shi et al., 2007), which manifests as plasticity. In general, long-term fluid extraction or injection activities can activate nonlinearity in both thematerialconstitutiveresponseandthekinematicsofdeformation. Useofelastic and infinitesimal deformation models in those cases will predict inaccurate stresses and, therefore, inaccurate onset, magnitude, and duration of the fault slip events. This is especially problematic for smaller magnitude earthquakes or aseismic slip events (Mw<4), which are not felt by humans but may cause damage to surface facilities and leakage of reservoir fluids (or other types of stored fluids such as natural gas or CO 2 ) along the activated fault into freshwater resources. Other manifestations of plastic deformation include drop in the reservoir permeability, loss of well productivity/injectivity, and wellbore failure (Allen, 1968; Fredrich et al., 2000). Once a rock fails plastically, its stress path can deviate significantly from that of an elastically deforming rock in reservoirs that are producing or injecting for decades. Plastic dissipation also changes the mechanical energy budget such that 86 energy available for brittle shear or tensile failure events can be expected to be smaller (Scholz, 2002). Given that fault slip is an irreversible deformation process and it is coupled to stresses in the host rock and driven by the elastic strain energy accumulated during production/injection, we hypothesize that the changes induced by plasticity in the stress path of the host rock can change the onset, location, magnitude and duration of the induced slip events. The present study tests this hypothesis. Accounting for the above changes can improve the agreement between model- predicted seismicity and ground deformation (Gonzalez et al., 2012; Bourne & Oates, 2017) and the observed data. Therefore, there is an urgent need to develop simulation methods that resolve coupled flow, geomechanics, and dynamic faulting processes in plastically deforming rocks and demonstrate their application to real reservoirs. Over the years, multiple coupled flow and elastoplastic geomechanical frameworks have been proposed (Armero & Simo, 1993; Borja & Alarcon, 1995; Armero, 1999; Sanavia et al., 2002; Borja, 2004; Li et al., 2005; Benallal et al., 2008), mostly in absence of fault slip. Here, we build on that foundation to develop such a method. Applying the method to faulted reservoir-caprock systems with realistic geological structure, heterogeneous rock properties, and multiphase flow is another challenge that we tackle. 4.2 Effect of finite strain plasticity on the flow problem It is important to understand the effect of plastic strain on the solution of the flow problem, especially within an iterative solution scheme where the flow and mechanics problems are solved sequentially. Therefore, we can re-write the fluid 87 mass increment Eqs. (2.20) explicitly in terms of the plastic strain. For example, for water, . M w =ρ 0 w   X N=o,w b w b N k dr +N wN ! . p N + b w k dr L v τ v +b w . ε p v   . (4.1) which can be combined with Eq. (2.19) to obtain the porosity evolution equation, 1 J . Φ = b k dr L v τ v + X N=o,w b 2 s N k dr + (b−φ 0 )s N k s ! . p N +b . ε p v (4.2) Substituting Eq. (4.1) in Eq. (2.5) gives us the water phase equation: X N=o,w b w b N k dr +N wN ! . p N + b w k dr L v τ v +b w . ε p v + 1 ρ w ∇·w w + q w ρ w = 0, (4.3) which is analogous to the pressure equation under infinitesimal deformation (Eq. (30) in (Jha & Juanes, 2014)). The oil phase equation is obtained simi- larly. For large displacement/rotation–small strain,L v τ v is substituted with . S v , and k dr is substituted with K dr . This suggests that the fixed stress iterative solu- tion scheme for infinitesimal deformation can be extended to finite strain plasticity by fixing . S v , . ε p v during the solution of the flow problem. 4.3 Analogy between fault slip and plastic yield- ing Mathematically, there are many similarities between fault failure and plastic yielding model (Owen & Hinton, 1980), e.g., the stick-slip condition f s ≤ 0 is analogous to the yield conditionf p ≤ 0, the slip rate . d s is analogous to the plastic strain rate . ε p , and the slip speed| . d s | is analogous to the plastic multiplier . λ. One 88 can also see analogies between the effective stress invariant I 0 1 = 3(S v +βP e ) and the effective normal tractionσ 0 n , √ J 2 and the shear traction|τ s |, the coefficient of internal friction α f and the fault friction coefficient μ f , and the internal cohesion- like parameter γ and the fault cohesion τ c . The analogy arises because plastic failure is usually identified with frictional sliding of solid grains along their contact surface, which is analogous to two sides of a fault sliding under friction. However, the analogy has physical limitations because seismic/aseismic slip may change the chemical composition of the fault zone material such that it is no longer the same as the host rock. In the Results section below, we exploit the analogy to relate the stress path of a material point inI 0 1 vs. √ J 2 space to the stress path of a point on a fault in the σ 0 n vs. |τ s | space. We use this to understand the impact of plastic failure on fault slip onset, location and magnitude. 4.4 Production-induced plastic failure and fault slip Motivated by the setup in Figure 2.1, we use the same 3D geomechanical model of a reservoir-caprock-basement system (Figure 4.1) as in Chapter 2. Similarly, the model is initialized under mechanical and hydrostatic equilibrium with a normal faulting stress regime. We apply a lateral compression on the east boundary (x- positive surface) which is 70% of the overburden to create a normal faulting stress regime. Both the compression and overburden increase linearly with the lithostatic gradient. Initial stresses are in equilibrium with the gravitational and boundary stresses such that the initial displacement is zero. All boundaries are no-flux boundaries. The reservoir is initially saturated with 12% water and 88% oil. The initial reservoir pressure is around 2075 psi. 89 Well Fault z (km) x (km) 1.4 1.6 1.5 0 2 1 b x (km) a Well Fault Tectonic compression Top point 0.5 2.5 y (km) 80 ◦ z (km) Reservoir cross-section (zoom-in) Reservoir 2 0 0 2 Footwall cell Hanging wall cell Figure 4.1: (a) Physical model with a well producing from a reservoir (gray color) confined between caprock and basement (white color) and intersected by a fault. (b) A vertical cross section of the reservoir showing the layered topography, well and fault. Position of the Top point and Hanging wall cell, Footwall cell, where detailed analysis of stress and slip will be carried out in the Results section below. We use the linear slip-weakening friction model to describe quasi-static fric- tional failure on the fault. The friction coefficient decreases from a static value μ s = 0.5 to a dynamic value μ d = 0.2 as the slip magnitude increases from zero to the critical slip distance d c = 0.1 m. The value of d c falls within the range of critical slip weakening distance reported in multiple studies (Scholz, 1989; Ide & Takeo, 1997; Olsen et al., 1997; Mikumo et al., 2003). We use a D-P perfectly plastic model where the yield surface touches the outer vertices of the M-C model, i.e. the inscribed D-P model. We use a friction angle ofφ f = 20 degree, dilatation angle of φ g = 20 degree, and cohesion of c = 1 MPa (Sarris & Papanastasiou, 2013). See the appendix for the relations between (α f ,α g ,γ) and (φ f ,φ g ,c). The model is populated with the following homogeneous isotropic properties: rock bulk density ρ b = 2460 kg/m 3 , initial permeability k(t = 0) = 5 mD, initial porosity φ(t = 0) = 0.1, saturated Biot coefficient b = 0.8, and drained Poisson’s ratio ν = 0.25. The overburden and underburden are hydraulically isolated from the reservoir. Transmissibilities of the fault in the three directions (lateral shear, dip 90 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (psi) 0 1000 2000 3000 1 1.002 1.004 1.006 1.008 1.01 0.39 0.395 0.4 0.405 0.41 (res. vol/surf. vol) (cp) Figure 4.2: (a) Relative permeability of waterk rw (solid) and relative permeability of oil k ro (dash) as functions of the water saturation s w . (b) The oil formation volume factor Bo and the oil viscosity μ o . shear, and normal directions) are set to zero. The water density is ρ w = 1000 kg/m 3 and viscosity is μ w = 0.4 cp at initial reservoir pressure. The oil density at the initial pressure is ρ o = 806.2 kg/m 3 . The oil relative permeability in presence of water k ro , the oil formation volume factor Bo, and the oil viscosity μ o are plot- ted as saturation- and pressure-dependent functions in Figure 4.2. The formation volume factor is defined as the ratio of the fluid density at the atmospheric pressure to the fluid density at reservoir pressure. Bo is related to the oil compressibility as c o =−(1/Bo)dBo/dp o . We consider four cases (listed in Table 4.1) by varying geomechanical proper- ties such as the drained Young’s modulus E and bulk modulus K dr . We produce oil at 2000 bbl/day (318 m 3 /day) for 850 days in the four cases. For each case, we run finite strain and infinitesimal strain simulations. To understand the effect of rock type on fault slip and plastic failure, we compare the spatial profiles and time evolution of reservoir pressure, vertical displacement, plastic magnitude and slip area in four cases. The slip area is important because it determines the seis- mic event magnitude as per Eq. (2.10). We analyze the stress path at different locations: the hanging wall cell, the footwall cell and the fault top point as shown 91 Table 4.1: Four cases of reservoir-caprock system defined in terms of their elastic moduli E (GPa) K dr (GPa) Case 1 Overburden and underburden 26.4 17.6 Reservoir 26.4 17.6 Case 2 Overburden and underburden 26.4 17.6 Reservoir 12.9 8.6 Case 3 Overburden and underburden 12.9 8.6 Reservoir 26.4 17.6 Case 4 Overburden and underburden 12.9 17.6 Reservoir 5.3 3.6 in Figure 4.1 to identify the effect of plastic failure on fault slip. We also analyze the stress path in the four cases to understand the rock type-induced temporal complexity of plastic failure and fault activation process. 4.4.1 Spatial profiles and time evolution We plot 3D fields of Case 1 results at t = 850 day from the finite strain elasto- plastic simulation in Figure 4.3 and Figure 4.4. Due to production, the reservoir pressure decreases and the water saturation increases near the well and propagates in the hanging wall side of the reservoir. Because of the impermeable fault, pres- sure and water saturation remain unchanged on the footwall side. The pressure remains almost constant in caprock and basement because they are hydraulically isolated from the reservoir. To understand the geomechanical signatures of reservoir pressure depletion, we plot vertical displacement of the whole field in Figure 4.4a. Pressure depletion leads to volumetric contraction of the reservoir and causes negative vertical dis- placements above the reservoir and positive displacement below the reservoir. The 92 -1500 -1000 -500 0 Pressure difference(psi) 0 0.5 1 1.5 Saturation difference 10 -3 a b Figure 4.3: 3D fields of (a) pressure change and (b) water saturation change at day 850. Pressure change is almost zero (yellow color) in the footwall block due to the low permeability of the fault. Water saturation increase due to oil production. 3×10 -4 -0.03 0.03 0 0 1.5×10 -4 0.44 (m) Plastic strain x z y z x a b 0 0.22 Slip d s (mm) c Figure 4.4: 3D fields of (a) vertical displacement, (b) plastic strain magnitude in the footwall block, and (c) slip magnitude on the fault surface at day 850. (a) shows a cut-away view of the u z field using a slice along the fault surface and a y = 1000 m slice. maximum downward displacement decreases from 3 cm above the reservoir to 1.2 cm on the top surface. We examine the evolution in stresses vs. average reservoir pressure space of the selected footwall block cell where plastic failure initiates and it’s corresponding hanging wall block cell in Figure 4.5. Production-induced reservoir contraction leads to negative normal stresses for the hanging wall block cells. However, for the footwall block cells, reservoir contraction in the hanging wall block applies 93 pull in the x direction and leads to an increase in S 0 xx . In y and z directions, footwall block cells experience smaller contraction than hanging wall block cells so S 0 yy and S 0 zz decrease slightly. The combination of normal stress changes cause a rapid decrease in effective volumetric stress for hanging wall block cells, and a small decrease in effective volumetric stress for footwall block cells (Figure 4.5d). Reservoir contraction increases the shear stress magnitude in both hanging wall and footwall blocks as shown in Figure 4.5c. The deviatoric stress invariant is dominated by changes in the normal stresses. The increase in S 0 xx overcomes the decrease in S 0 yy and S 0 zz and leads to a large increase in the deviatoric stress in the footwall block. However, for hanging wall block cells, the decrease in S 0 xx balances out with the decrease in S 0 yy and S 0 zz and the deviatoric stress remains almostunchanged(Figure4.5bandd). Thus, plasticfailureinitiatesinthefootwall block and propagates laterally and vertically. The plastic strain magnitude in the footwall block is plotted in Figure 4.4b. Reservoir contraction also applies downward pull on the reservoir top surface and upward pull on the reservoir bottom surface. This is production-induced shear. Since the boundary conditions in Figure 4.1(a) are such that all points on the fault experience downdip tectonic shear, both the production-induced shear and the tectonic shear are acting downdip on the top surface. As a result, one of the points on the reservoir top is most likely to slip compared to all other locations on the fault. This is confirmed by Figure 4.4c which shows that the slip nucleates near the top of reservoir. To compare the results among finite strain and infinitesimal strain simulations, we keep the pressure drop, which is the driving force behind induced seismic- ity, the same between the two simulations (Figure 4.6a). The distributions of 94 -70 -65 -60 -55 -50 6 6.5 7 7.5 8 8.5 9 -70 -65 -60 -55 -50 6 6.5 7 7.5 8 8.5 9 Footwall cell Hanging wall cell Failure line 12.5 13 13.5 14 14.5 0 0.05 0.1 0.15 0.2 0.25 12.5 13 13.5 14 14.5 0 5 10 15 20 25 30 12.5 13 13.5 14 14.5 -25 -20 -15 -10 a b Time d c 12.5 13 13.5 14 14.5 0 0.02 0.04 0.06 12.5 13 13.5 14 14.5 -20 -18 -16 -14 -12 -10 12.5 13 13.5 1 0 5 10 15 20 25 30 Well Fault z (km) x (km) 1.4 1.6 1.5 0 2 1 b x (km) Well Fault Tectonic compression Top point 0.5 2.5 y (km) 80 ◦ z (km) Reservoir cross-section (zoom-in) Reservoir 2 0 0 2 Footwall cell Hanging wall cell Figure 4.5: Case 1 stress components evolution of (a) effective normal stress, (b) and (c) deviatoric stress versus pressure for selected footwall block cell (black line) and corresponding hanging wall block cell (blue line). (d) Stress paths of the (a) footwall block cell (black points) and hanging wall block cell (blue points) from Case 1. displacement, plastic strain and fault slip remain similar in both cases to Fig- ure 4.4, which implies that finite strain plasticity did not impact the location of the hypocenter. However, the magnitudes of these quantities are different in the two simulations which implies that the seismic/aseismic slip onset time and the magnitude are different. This is an important result that can explain why some of the induced seismicity simulation models, which assume infinitesimal strain over decade-long production from a oil/gas reservoir or injection in an aquifer, under- or over-estimate the timing and magnitude of induced events. Compared with the 95 10 11 12 13 14 0 0.2 0.4 0.6 0.8 1 1.2 10 11 12 13 14 0 1 2 3 4 5 6 10 -4 10 11 12 13 14 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 c d Finite strain elastoplastic Infinitesimal strain elastic 10 11 12 13 14 -0.1 -0.08 -0.06 -0.04 -0.02 0 Time Reservoir top point Model top point 0 200 400 600 800 1400 1600 1800 2000 2200 12 12.5 13 13.5 14 14.5 -52.55 -52.5 -52.45 -52.4 -52.35 -52.3 52.25 Case 1 Case 2 Case 3 Case 4 a Time (day) b e Slip Area (x10 m) Time Time Figure 4.6: Time evolution of (a) the average reservoir pressure. Evolution of (b) the vertical displacement at the top of the reservoir and (c) at the top of the model, (d) plastic strain magnitude at the footwall block cell, and (c) slip area on the fault surface versus the average reservoir pressure for the four cases with finite strain and infinitesimal strain simulations. finite strain simulation, as shown in Figure 4.6, the vertical displacement in the infinitesimal simulation is smaller for the same amount of pressure drop because of a lack of plastic yielding. In the infinitesimal simulation, compared to the finite strain elastoplastic simulation, the fault slips earlier, at a smaller pressure drop value, and generates a larger slip area (larger M w ) at the end of production. With different rock types, we observe spatial distributions of pressure, satu- ration, displacement, plastic strain and fault slip that are similar to Figure 4.4. However, the time evolutions are different and revealing. With the same amount of oil produced, when the reservoir rock has a smaller bulk modulus and Young’s modulus, we observe less reservoir pressure drop in Case 2 and Case 4 compared to Case 1 (Figure 4.6a). When the overburden and underburden rocks have smaller bulk and Young’s moduli compared to those of the reservoir (Figure 4.6b Case 3), the reservoir pressure drop is again smaller. To identify the effects of pressure 96 depletion on geomechanical behavior, we plot the evolution of vertical displace- ment and plastic strain magnitude of selected points as well as the fault slip area as functions of the average reservoir pressure in Figure 4.6b to e. A weaker reser- voir with smaller elastic moduli can be compressed more easily. So, with the same amount of pressure decline, we observe larger vertical deformation in Case 2 and Case 4 compared with Case 1 at both the reservoir top and model surfaces. The mechanics-to-flow coupling (stress and strain terms in Eq. (4.3)) reduces the drop in reservoir pressure during production. This mechanism is also known as the compaction drive in petroleum engineering literature (Aziz & Settari, 1979). Com- paring Case 1 and Case 3 shows interesting results of practical significance. We observe that Case 3 vertical displacement is larger than Case 1 vertical displace- ment at the reservoir top and smaller than Case 1 vertical displacement at the model top. The reason is as follows. The caprock in Case 3 has a smaller modulus and it deforms more. However, the downward pull on the caprock decreases as we move away from the reservoir top towards the ground surface because the pressure change is localized inside the reservoir and the induced stresses decay rapidly with distance. With smaller elastic moduli, the plastic failure initiates at a smaller pressure drop as shown in Figure 4.6d. The plastic strain magnitude also increases faster in Case 2 and Case 4 compared to Case 1. In Case 3, a larger portion of the pressure depletion-induced stresses are spent in driving caprock and basement deformations, so we observe a delayed onset of plastic failure and a smaller plastic strain magnitude. The fault slips earlier at a smaller pressure drop in Case 2 compared to Case 1. Since the pressure drop is smaller in Case 3 and Case 4, we donotencounterfaultslipduringthesimulationperiod. Theslipisfurtherdelayed in Case 3 because caprock and basement rocks are weaker (smaller moduli). 97 4.4.2 Effect of plastic failure on fault slip To further understand the relationship between plastic failure and fault slip, we show the results from finite and infinitesimal strain simulations for Case 1. We first examine the time evolution in the space of the deviatoric stress invariant √ J 2 vs. the effective volumetric stress invariant I 0 1 . We consider the selected footwall block cell. As shown in Figure 4.7a, production-induced compaction leads to negative I 0 1 in both finite and infinitesimal strain simulations. Before the onset of plastic failure, √ J 2 increases in the footwall block cell and later reaches the D-P yield surface. In the finite strain simulation, √ J 2 decreases andI 0 1 increases (lower compression) after reaching the yield surface. In infinitesimal strain simulation, √ J 2 continues to increase and I 0 1 continues to decrease (higher compression). In Figure 4.7b, we show the time evolution in the shear traction (τ here refers to|τ s | in Eq. (2.21)) vs. the effective normal traction (σ 0 n ) space of the hypocenter (top point on the fault surface). Reservoir contraction causes tension on the fault surfacewhichincreasesσ 0 n andcausesdowndipshearonthetoppointwhichreduces τ. Initially, the stress paths of both elastoplastic and elastic simulations towards the M-C failure line are similar. However, when plastic failure occurs, since plastic failure lowers the compression, σ 0 n increases rapidly while τ continues to decrease. When plastic failure propagates to the overburden and reservoir boundary, σ 0 n increases at a slower rate and τ continues to decrease. This results in a delayed fault slip and a smaller magnitude seismic event in the elastoplastic simulation compared to the elastic simulation. 4.4.3 Effect of rock type Duetoproduction-inducedcontraction,inducedstressesdependonthepressure drop. We further analyze the evolution of stress as a function of the average 98 -53 -52.5 -52 -51.5 -51 5.5 6 6.5 7 7.5 8 8.5 -10.2 -10 -9.8 -9.6 -5 -4 -3 -2 -1 Time a Time -10.2 -10 -9.8 -9.6 -9.4 -5 -4 -3 -2 -1 Finite strain elastoplastic Infinitesimal strain elastic Failure line b Figure 4.7: Stress paths of the (a) footwall block cell and (b) top point on the fault surface from Case 1 finite strain and infinitesimal strain simulations. The solid line in each plot is the M-C failure line. The arrows indicate the direction of time. The red star marks the onset of plastic failure. reservoir pressure from the finite strain simulation for different rock types at the selected footwall block cell and the fault surface top point. Initially, the model is under equilibrium with the same initial reservoir pressure and initial I 0 1 and √ J 2 in the four cases. Production reduces the reservoir pressure. Before plastic failure begins, production-induced reservoir contraction leads to decrease inI 0 1 and increase in √ J 2 . As shown in Figure 4.8a and b, with the same amount of pressure drop, I 0 1 decreases less while √ J 2 increases more when the reservoir has smaller elastic moduli. This leads to a relatively faster growth in the yield function f p , definedinEq.(2.38),forCases2and4comparedtoCase1,asshowninFigure4.9a. When overburden and basement have smaller moduli, comparing Case 3 and Case 1, the pulling effect due to compaction of the hanging wall reservoir is shared by overburden and basement deformations such that I 0 1 of the footwall cell decreases faster and √ J 2 increases slower. The yield function f p grows relatively slower for Case 3 and it requires additional pressure drop to reach plastic failure, as shown in Figure 4.9a. 99 10 11 12 13 14 -5 -4 -3 -2 -1 10 11 12 13 14 -10.2 -10 -9.8 -9.6 -9.4 10 11 12 13 14 5.8 6 6.2 6.4 6.6 6.8 7 10 11 12 13 14 -52.5 -52 -51.5 -51 a Time 12 12.5 13 13.5 14 14.5 .55 2.5 .45 2.4 Case 1 Case 2 Case 3 Case 4 b c d Figure 4.8: Effect of rock type. Evolution of (a) effective stress invariantI 0 1 and (b) deviatoric stress invariant √ J 2 in the footwall cell, and evolution of (c) effective normal fault traction σ 0 n and (d) shear traction τ at the top point as functions of the average reservoir pressure p res from the four cases in Table 4.1. The arrows indicate the direction of time. We highlight the initiation of plastic failure with red star and fault slip with red pentagon in Case 1. In the √ J 2 vs. I 0 1 space shown in Figure 4.10a, as I 0 1 decreases, √ J 2 increases faster for cases with smaller reservoir moduli. Since the yield condition requires less √ J 2 at higherI 0 1 , cases with smaller reservoir moduli tends to fail at a smaller pressure drop. On the other hand, when caprock and basement moduli are smaller, √ J 2 grows at a smaller rate with decreasingI 0 1 . This means a larger pressure drop is required to reach the D-P failure line. Reservoir is relatively stable. After plastic failure, stresses stays on the failure line. I 0 1 increases and √ J 2 decreases. 100 We show the evolution ofσ 0 n andτ as functions of the average reservoir pressure from finite strain simulation for different rock type cases at the selected top point in Figure 4.8c and d. As pressure drops, contraction-induced tension on the fault surface leads to an increase inσ 0 n for the four cases. σ 0 n increases slightly faster for cases with smaller reservoir moduli (Case 2 and Case 4) while it grows slower for cases with smaller caprock and basement moduli (Case 3). After plastic failure, σ 0 n increases at a higher rate. For Cases 2 and 4, which have smaller reservoir moduli, σ 0 n increases faster than Cases 1 and 3. After plastic failure reaches the reservoir top boundary, σ 0 n continues to increase with a slower rate in the four cases. Reservoir compression causes downdip shear at the top of the reservoir which leads to a decrease in τ at the reservoir top point. With smaller reservoir moduli, the reservoir contracts more for the same value of pressure drop, the top point experiences larger downdip shear, and τ decreases more as shown in Cases 2 and 4. We defined the Coulomb Failure Function (CFF) in Eq. (2.21). In producing reservoirs, the change in τ dominates the evolution of CFF on faults because the fault pressure p f in σ 0 n (Eq. (2.18)) remains either constant or drops slightly. This is because p f , during production, is proportional to the pressure on the non- producing side of the fault, so it can correctly predict the slip onset time (Jha & Juanes, 2014; Yang & Juanes, 2018). We observe that for a given pressure drop, CFF grows faster in cases when the elastic moduli of the reservoir rock is smaller. It leads to an earlier onset of slip in Cases 2 and 4 with a smaller pressure drop. On the other hand, when caprock and basement have smaller moduli, bothσ 0 n and τ show smaller variation for a given drop in the reservoir pressure. Initiation of fault slip requires higher pressure depletion as shown in Figure 4.9. 101 10 11 12 13 14 -1 -0.8 -0.6 -0.4 -0.2 0 a Time 10 11 12 13 14 -4 -3 -2 -1 0 b Time 12 12.5 13 13.5 14 14.5 -52.55 -52.5 -52.45 -52.4 Case 1 Case 2 Case 3 Case 4 Figure 4.9: Effect of rock type on the rate of failure. Evolution of (a) yield function f p and (b) Coulomb Failure Function CFF, as a function of average reservoir pressure p res from four cases. The arrows indicate the direction of time. -53 -52.5 -52 -51.5 -51 6 6.2 6.4 6.6 6.8 7 -10.2 -10 -9.8 -9.6 -5 -4 -3 -2 -1 Time a Time b -10.2 -10 -9.8 -9.6 -9 -5 -4 -3 -2 -1 Case 1 Case 2 Case 3 Case 4 Failure line Figure 4.10: Stress paths of the (a) footwall block cell and (b) top point on the fault surface from the four cases. The solid straight line in (a) is the D-P failure line and in (b) is the M-C failure line. The arrows indicate the direction of time. As shown in Figure 4.10b, before plastic failure, in the cases where reservoir rocks have smaller moduli, τ decreases faster with increasing σ 0 n . Plastic failure leads to a further increase in σ 0 n . Because the decrease in τ is larger than the increase in σ 0 n , Cases 2 and 4 reach M-C line earlier. For weaker caprock and basement rocks, τ decreases at a smaller rate with increasing σ 0 n . The fault is relatively stable. 102 4.5 Conclusion We presented a novel simulation framework of coupled multiphase flow, finite strain elastoplasticity and fault slip in petroleum reservoirs. The framework incor- porates realistic fault constitutive behavior by using a slip-weakening dynamic friction model. We applied the framework to synthetic examples as well as to a production-induced seismicity case. We investigated the effect of plastic failure on fault stability. We studied the impact of variation in the elastic moduli between reservoir and caprock and how this affects production-induced plasticity and slip. The main highlights are as follows: 1. The poroplastic reservoir exhibits larger vertical deformation and delayed slip than the poroelastic reservoir after the same amount of oil production because of plastic yielding. 2. Plastic dissipation releases a part of the mechanical energy generated by production-induced stresses. This reduces the energy available for seis- mic/aseismic slip and results in a delayed slip and a smaller magnitude slip event. 3. Reservoirs with smaller elastic moduli have a smaller pressure depletion for the same amount of production. 4. For the same amount of pressure drop, when the reservoir rock has smaller elastic moduli than those of the caprock, the vertical displacement is larger and plastic failure and fault slip begin earlier. When the caprock moduli are smaller than those of the reservoir, the vertical displacement is larger on the reservoir top surface and smaller on the ground surface (masking effect of 103 the caprock). A higher pressure drop is required to activate plastic failure and fault slip in this case. 104 Chapter 5 Role of Well Operations and Multiphase Geomechanics on Fault Stability 5.1 Introduction It is known that production or injection of fluids in subsurface reservoirs can change stress inside and outside the reservoir. Such changes may perturb mechanical equilibrium and cause slip on the pre-existing faults or trigger earth- quakes (Yerkes & Castle, 1976b; Segall, 1989b; Segall & Fitzgerald, 1998; Evans et al., 2012; Won-Young, 2013; Segall & Lu, 2015a; Won & Paul, 2016; Albano et al., 2017b). Coupling between fluid flow and mechanical deformation in porous media has been used in recent studies for assessment of fault activation and induced seis- micity (Jha & Juanes, 2014; Juanes et al., 2016; Fan et al., 2016; Panda et al., 2018). Various coupled multiphysics numerical simulation frameworks have been proposed for poroelastic media with mechanical discontinuities such as faults and fractures (Mainguy & Longuemare, 2002; Minkoff et al., 2003b; Fan et al., 2016; Garipov et al., 2016). However, it still remains a challenge to apply such modeling and simulation frameworks to real fields with complex structural and stratigraphic features e.g. intersecting faults and layers, petrophysical heterogeneity, and time- dependent and spatially distributed production and injection. There is an urgent 105 need to develop coupled modeling frameworks that can efficiently handle realistic geology, heterogeneous properties, and complex well schedules. Injection and storage of CO 2 in depleted oil reservoirs is an important tech- nology that could reduce greenhouse gas emissions while enhancing the recovery of hydrocarbons (Lackner, 2003; Franklin, 2009; Szulczewski et al., 2012) from the oilfield. However, injection-induced pressurization, simultaneous movement of CO 2 and hydrocarbon phases, and sudden changes in fluid compressibility due to CO 2 dissolution in oil (Ampomah et al., 2016) can lead to changes in the stress within and around the reservoir. This could potentially cause irreversible changes in the poromechanical properties of the rock, opening of fractures, or activation of critically-stressed existing faults. Such changes may affect the oil recovery, cre- ate undesired pathway for CO 2 leakage or induce seismic events (Streit & Hillis, 2004; Rutqvist et al., 2007; Chiaramonte et al., 2008; Rutqvist et al., 2008; Vidal- Gilbert et al., 2009, 2010; Bissell et al., 2011; Preisig & Prévost, 2011; Verdon et al., 2011; Olden et al., 2012; Castelletto et al., 2013; Teatini et al., 2014). Thus, it is important to monitor the evolution of flow and mechanical behavior of the overburden-reservoir-basement storage complex during CO 2 injection and oil pro- duction. Here, we report results from a coupled multiphase flow-geomechanical modeling and monitoring study of the Farnsworth Unit (FWU) oilfield of Texas to understand the evolution of pressure, oil saturation and stress in the storage complex and to monitor the mechanical stability of reservoir-intersecting faults. 5.2 Farnsworth Unit Geology The Farnsworth Unit field is located in Ochiltree County, Texas, which is on the northeast shelf of the Anadarko Basin (Figure 5.1) (El-kaseeh et al., 2017). 106 Figure 5.1: Map of the Farnsworth Unit (FWU) field. Left: The border of FWU and some of the well locations are shown. Right: FWU is located in the northern Texas within the Anadarko Basin (shaded region), and it is shown within a green square. The main source of anthropogenic CO 2 used during injection is highlighted by two star-shaped symbols: the northern star is the Arkalon Ethanol Plant in Liberal, Kansas, and the southern star is the Agrium Fertilizer Plant in Borger, Texas. The Upper Morrow Sandstone, which is the focus of this study, consists of 115 m to 150 m thick Morrow Shale and discontinuous coarse-grained Morrow-B Sand- stone. The Morrow-B Sandstone layers with typical thicknesses between 6.7 m and 15.2 m create the hydrocarbon-bearing reservoir. Conformably overlying the Upper Morrow is the Thirteen Finger Limestone and the lowest member of the Pennsylvanian Atoka Series. The base of Thirteen Finger Limestone consists of a thin coal bed overlain by fine crystalline dolomitic limestones and shales. Morrow Upper Shale and Thirteen Finger limestone serve as the caprock formations in FWU. Faulting in the Farnsworth Unit closely resembles the deformation in the nearby Wichita Megashear zone, which is a major sinistral strike-slip fault system in the foreland of the Amarillo-Wichita Mountain Front in the southern Anadarko Basin(Evans,1979). TheWichitaMegashearexperiencedverticalblockmovement followed by lateral strike slip movement. The vertical fault movement is expressed 107 Figure 5.2: Map view showing the faults within FWU at the Morrow Sandstone depth level. The FWU boundary is shown with thick black lines. Township and range boundaries are shown in a light gray color. Nine faults are interpreted using FWU’s 3D seismic survey. Except Fault-6, all other faults intersect both the Morrow-B Sandstone and Thirteen Finger Limestone. as normal faults with a southward displacement (Ampomah et al., 2016). The left lateral strike slip movement created anticline folding. A total of nine faults were interpreted from 3D surface seismic data in FWU (Figure 5.2). Fault-1 through Fault-5 are interpreted as part of a sinistral wrench fault system. Except Fault-6, all faults intersect both Morrow-B Sandstone and Thirteen Finger Limestone. The Farnsworth Unit was discovered in 1955 and unitized in 1963. Water injection-based secondary recovery commenced in 1964. Between 1955 to 2010, approximately 38 million stock tank barrels (MMstb, 1 stb = 0.159 m 3 ) of oil and 28 thousand standard cubic feet (Mscf, 1 scf = 0.0283 m 3 ) of natural gas was produced from FWU. CO 2 injection started in December 2010 as part of tertiary recovery in the oilfield, and it was implemented through five-spot well patterns. Since 2013, three to five new well patterns have been added each year for CO 2 storage-EOR (Figure 5.1 above). Anthropogenic CO 2 is injected at a rate of about 108 518 tonnes/day or about 0.2 million metric tonnes/year. During October 2010– July 2016, 13 injectors injected more than 0.8 million metric tonnes of CO 2 , which helped to recover additional 1.5 million barrels (MMbbls) of oil from FWU. 5.3 Numerical Model Using available seismic and well logging information in FWU, we create a 3D geological model of reservoir-overburden-underburden system. In this study, we focus on the western part of FWU because CO 2 injection is confined to that region. Threemajorstratigraphichorizonsareconsideredinthemodel: MorrowShaleTop, Morrow-B Sandstone Top and Morrow-B Sandstone Base as shown in Figure 5.3. Morrow shale serves as the caprock while Morrow-B sandstone is the reservoir layer subjected to production and injection. The depth of each horizon is estimated by averaging the well marker depths from logs (Ampomah et al., 2016; Balch et al., 2017). These depths are shown in Figure 5.3 along a cross-section A-A’ defined in Figure 5.2. Three major faults compartmentalizing mechanical and hydraulic response of the oilfield are included in the model. 5.3.1 Mesh The model dimensions are 10 km× 10 km× 1 km in x, y and z directions. The depth interval of the model is fromz=1 km toz=2 km below the ground sur- face. Three major faults are considered: Fault-1, Fault-3 and Fault-6. According to White et al. (2017), Fault-3 is almost vertical with a small dip to the south. We assume the dip angle to be 6 ◦ as shown in Figure 5.3. Similarly, Fault-1 is bending to west with an angle of 6 ◦ . Both Fault-1 and Fault-3 cross through the three horizons. However, Fault-6 is younger i.e. more recent than Morrow-B Sandstone 109 (a) (b) x y z y Fault 1 Fault 3 Fault 6 A A' A' A Morrow Shale Morrow B top Morrow B base Fault 3 Figure 5.3: (a) Top view and (b) side view along cross-section AA’ show the model geometryandthemechanicalboundaryconditions. Fault-1andFault-3arecutting through the model while Fault-6 ends below the reservoir layer. Fault-3 is located near the center of the model and dips towards south with a 6 ◦ angle as shown in (b). Three layers are highlighted in (b): Morrow Shale, Morrow-B Sandstone Top and Morrow-B Sandstone Base. (a) The minimum principal stress is applied on the y-positive surface. The maximum principal stress is applied on the x-positive surface. (b) A normal compression equal to the overburden is applied on the top boundary. and terminates just below the reservoir. We assume Fault-6 to be vertical based on regional geology. Since Fault-1 and Fault-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 by splitting Fault-3 into two smaller faults, Fault-3L and Fault-3R, in the model. The resulting gap between Fault-3L and Fault-1, and between Fault-3R and Fault-1 is approximately 350 m long (it is smaller in the refined mesh described in SI). Fluids can move through the gap to establish partial hydraulic communication across the faults, which is also supported by a previous flow simulation study in FWU (Ampomah et al., 2017a). 110 Conformal 3D meshing of a domain containing multiple arbitrarily-shaped and intersecting surfaces is an computational challenge and an open research problem. According to one estimate (SpaceClaim, 2018), analysts can spend upwards of 75% of their time trying to fix, simplify and work with complex geometric objects prior to meshing. The fault meshing challenge has led to exclusion of many (Juanes et al., 2016) or all (Hornbach et al., 2015b; Lei et al., 2017) intermediate-scale faults, which are potential sources of induced seismicity, from prior modeling stud- ies. We developed a generalized meshing workflow that can honor intersections of multiple curved faults and layers in 3D. Our workflow, which is implemented in Trelis (Trelis, 2017), has the following steps. First, we create non-uniform ratio- nal B-spline (NURBS) representation of fault surfaces and stratigraphic horizons. Then we build a 3D solid geometry of faults/horizons using de-featuring and topo- logical operations of webcut, chop, trim, merge and imprint to create hanging and footwall fault blocks and layers in the reservoir, overburden, and underburden. We remove small-scale features that cause failure of the meshing algorithm. We use a bottom-up meshing approach in which 1D curve objects of a horizon or a fault surface are meshed first followed by meshing of the parent 2D object. Mesh- ing of a selected number of 2D surfaces in a 3D fault block guides the meshing of the 3D block itself, which is followed by meshing of the entire 3D solid vol- ume. We use tetrahedral elements to honor arbitrary 3D geometry of faults and achieve conformance between the mesh and the structural geology of the region. This increases the accuracy of fault stresses predicted to result from CO 2 injection and hydrocarbon production. We use adaptive meshing to generate coarser mesh near the domain boundaries and in hydraulically-inactive regions, and finer mesh within active flow regions e.g. within thin reservoir layers, and around wells and 111 Figure 5.4: Computational domain mesh showing the locations of the fault surfaces and the wells. We use tetrahedral elements with smaller-sized elements within the reservoir layers and next to the fault surfaces and larger-sized elements near the domain boundaries. 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. faults. This reduces the total number of mesh elements and enables faster simula- tions. This model has a total of 33760 tetrahedral elements as shown in Figure 5.4. The average element volume is approximately 1.8× 10 5 m 3 in the reservoir and 3.6×10 6 m 3 outside of the reservoir. To analyze the mesh sensitivity of the results and conclusions, we created a finer mesh of the domain with 261207 elements that is described in SI. Consistent finite element and finite volume mesh files are created from both meshes to provide input to the coupled simulator. 112 5.3.2 Fault model There are two stages in fault modeling: mesh representation and dynamic fail- ureoffaults. Thefiniteelementmeshisprocessedtointroducezerothicknessinter- face elements and cohesive nodes on each fault surface (Aagaard et al., 2013). This allows representing the fault as a surface of displacement discontinuity embedded in the continuum and calculation of the fault slip and traction vector components in the three directions–along strike, along dip, and outward normal (Jha & Juanes, 2014). We solve the nonlinear contact problem subject to the Mohr-Coulomb frictional failure condition at the fault cohesive nodes to model fault slip (Aagaard et al., 2013; Jha & Juanes, 2014). 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 change in mechanical stability induced by production and injection. The likelihood of fault failure is assessed by evaluating the change in Coulomb Failure Function, ΔCFF, which is a function of the changes in pressure and stress around the fault (Reasenberg & Simpson, 1992; Jha & Juanes, 2014): ΔCFF = CFF(t)−CFF(0) = (|τ| +μ f σ 0 n ) t − (|τ| +μ f σ 0 n ) 0 = hq τ 2 lat +τ 2 dip t − q τ 2 lat +τ 2 dip 0 i +μ f [(σ 0 n ) t − (σ 0 n ) 0 ] = Δ|τ| +μ f Δσ 0 n = Δ|τ| +μ f (Δσ n +bΔp f ). (5.1) Above, Δ|τ| is the magnitude of change in the shear traction vector. μ f is the fault friction coefficient. Δσ 0 n = Δσ n +bΔp f is the change in effective normal traction, Δσ n is the change in total normal traction and Δp f is the change in fault pressure. Here, we consider tension to be positive and compression to be negative, so a 113 positive ΔCFF indicates a destabilizing effect on the fault. We consider all shear to be destabilizing and use the magnitude of shear traction in ΔCFF calculation. Useofinterfaceelementstorepresentafaultallowsustoreproducethepressure jump across the fault, which is observed in the field as a difference in the flow rate and pressure in wells located across the faults that compartmentalize the reservoir. Because of the pressure discontinuity, fault stability can be assessed by evaluating the Coulomb failure criterion on both the hanging wall and the footwall sides of the fault separately. Since the failure criterion is met first with the larger of the two pressures, we define the fault pressurep f to be the larger of the two pressures as proposed recently (Jha & Juanes, 2014; Yang & Juanes, 2018; Jagalur-Mohan et al., 2018). The value of the Biot coefficientb is usually high in fault zones. Ifb is assumed to represent an unconsolidated fault zone undergoing rupture, or if a conservative estimate of the induced event trigger time is desired, then b can take values near its upper limit i.e. b=1 (Segall & Fitzgerald, 1998). We assume a uniform static friction coefficient value of μ f =0.6 on all the faults. The friction coefficient value does not affect the results of this study because the faults do not slip, and we focus on the relative change in values of the state variables. 5.3.3 Petrophysical properties The model is populated with poromechanical properties obtained using two representative well logging datasets. Different reservoir properties are assigned alongandacrossthelayerstomodelanisotropicresponse. Theaveragebulkdensity is ρ b =2480 kg/m 3 , the average sonic compressional velocity is V p =3600 m/s, the 114 average sonic shear velocity isV s =1920 m/s. Dynamic values of the bulk modulus and Young’s modulus can be calculated as follows: K =ρ b V 2 p − 4 3 V 2 s , E = ρ b V 2 s (3V 2 p − 4V 2 s ) V 2 p −V 2 s (5.2) 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 (Ampomah et al., 2017b) on to the coupled geome- chanical model. The older model is constructed in Eclipse (Schlumberger, 2017), is uncoupled to mechanics and has been used to track CO 2 migration 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. Quality of the property mapping algorithm is evaluated by comparing the property histograms of the two models (Figure 5.5). We observe that the two histograms have similar values for the average permeability and the permeability variance. During model calibration, also known as history-matching, we increased the permeability values by a factor of three to honor the historical production and injection rates and pressures. The average permeability in the x- and y-directions is approximately 189 mD (1 mD = 9.87× 10 −16 m 2 ), and the average permeability in the z-direction is approximately 144 mD. The permeability field is assumed to 115 -1 0 1 2 34 0 0.02 0.04 0.06 0.08 0.1 0.12 Probability -1 0 1 2 34 0 0.02 0.04 0.06 0.08 0.1 0.12 (a) (b) (d) -1 0 1 2 34 0 0.02 0.04 0.06 0.08 0.1 0.12 Existing Model Geomechanical Model Existing Model Geomechanical Model Figure 5.5: 3D mapping of rock properties from an existing flow-only uncoupled model into the geomechanical model. Histograms of log-normalized permeability in (a) x-direction, (b) y-direction, and (c) z-direction from the uncoupled Eclipse model (blue color) and the coupled geomechanical model (orange color). be constant in time. This is a reasonable assumption in Farnsworth because of the consolidated lithology of the formations and a relatively small change in the average reservoir pressure, which is actively maintained through waterflooding and CO 2 injection. Pore volume is one of the key properties of any reservoir simulation model because it is related to the original hydrocarbon volume, ultimate recovery from primary depletion, and the rate of change in pore pressure due to produc- tion/injection. The pore volume is a product of porosity and the bulk volume of the reservoir. 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. The reservoir bulk volume is determined by the mesh structure. We use a three phase–three component Black-Oil model (Aziz & Settari, 1979) to represent the multiphase flow system during flow simulation. The Black-Oil flow model must be initialized with phase-dependent fluid properties: compressibility, viscosity, and the dissolved phase volume fraction as functions of the fluid phase 116 0 1000 2000 3000 Pressure (psi) 1 1.4 1.8 0 1000 2000 3000 Pressure (psi) 1 1.004 1.008 Bw 0 2000 4000 6000 8000 Pressure (psi) 0 0.004 0.008 Bg (a) (b) (c) Figure 5.6: Formation volume factors for (a) oil, (b) water, (c) gas phases. The bubble point pressure of oil is approximately 1900 psi (1 psi = 6894.76 Pa) below which dissolved gas is released from the oil phase. We assume a constant water compressibility, which corresponds to a constant dB w /dp. We use CO 2 properties for the gas phase instead of the reservoir gas properties. This allows us to honor the CO 2 injection and transport behavior in the field more accurately. pressure. We provide these properties for oil, water and gas phases that are rep- resentative of FWU fluids at reservoir conditions. The Formation Volume Factors of the fluid phases, which determine their compressibilities (Aziz & Settari, 1979), are plotted in Figure 5.6. The multiphase flow model is defined with saturation- dependent relative permeability tables for oil, water and gas phases (Figure 5.7). The relative permeability curves for water (k rw ), for oil in presence of water (k row ), for oil in presence of gas, (k rog ), and for gas (k rg ) are based on lab measurements and calibration of the Eclipse flow model to measured production volumes in the field. We use the segregation model to calculate the three-phase oil relative per- meability (k ro ) as a function ofk row ,k rog , and saturations of water and gas phases. Transmissibilities of the three faults in the three directions (lateral shear, dip shear, and normal directions) are set to zero. Partial hydraulic communication across the faults exists through the gap between the faults as mentioned in the Mesh section above. This honors the fault gouge calculations (Hutton, 2015) and pressureobservationsinthefield. Comparedtothescenariowithconductivefaults, zero-transmissibility faults cause the reservoir pressure changes resulting from fluid 117 0 0.2 0.4 0.6 0.8 1 Sg 0 0.2 0.4 0.6 0.8 1 krog 0 0.2 0.4 0.6 0.8 kro 0 0.2 0.4 0.6 0.8 1 Sw 0 0.2 0.4 0.6 0.8 1 krw 0 0.2 0.4 0.6 0.8 1 krow (a) (b) Figure 5.7: The three-phase (oil, water and gas) relative permeability model in our simulation is based on two two-phase relative permeabilities. (a) Relative permeability of water k rw (dash) and relative permeability of oil in presence of waterk row (solid)asfunctionsofthewatersaturationS w . (b)Relativepermeability of gas k rg (dash) and relative permeability of oil in presence of gas k rog (solid) as functions of the free gas phase saturation S g . The relative permeability of oil k ro is calculated dynamically during the simulation as a function ofk row ,k rog , and the saturations. 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 2 injection pressure and storage capacity. 5.3.4 Initial and boundary conditions The flow-geomechanical model is initialized under mechanical and hydrostatic equilibrium with normal faulting stress regime imposed via appropriate tectonic stresses at the boundaries. Following Ampomah et al. (2016), we set the initial oil saturation to 0.69 and the initial water saturation to 0.31 for the whole model. Based on Figure 5.7, the water phase is immobile at this saturation. The initial pressure follows the hydrostatic pressure gradient as expected based on initializa- tion of the flow model. The initial stress field is set to honor the oblique strike-slip normal faulting stress regime of Fault-3, which strikes east-west. The minimum horizontal principal stress is assumed to be 0.9 times the maximum horizontal principal stress, which we assume to be equal to the vertical stress. The minimum 118 stress is applied on the northern boundary (y-positive surface) and the maximum principal stress is applied on the eastern boundary (x-positive surface). The ver- tical stress is assumed to be lithostatic: σ v =−g R ρ b (z)dz (compression negative), where g is the value of gravitational acceleration. A normal compression equal to the overburden weight is applied on the top boundary at z=1 km (Figure 5.3). 5.3.5 Production and injection history The coupled fluid flow and geomechanical simulation spans approximately 60 years of production and injection from January 1956 to July 2016. There are 76 wells in the field: 45 producers and 31 injectors (Figure 5.4). Production and injection profiles of typical production and injection wells are plotted in Figure 5.8. Cumulativeliquid(oilandwater)productionattheendofthesimulationisapprox- imately 60 MMstb (million stock tank barrel). Water injection starts in year 1967 and the cumulative injection volume is approximately 70 MMstb. Gas injection starts in 2010 and the cumulative gas injection volume is approximately 30 MMscf (million standard cubic feet). Since the objective of this study is to identify the fault stability mechanisms arising out of the coupling between multiphase flow, geomechanics and well operation, we do not focus on matching these numbers on cumulative production and injection volumes exactly in our simulation. 5.4 Results 5.4.1 Field results Evolution of the average pressure in the reservoir is plotted in Figure 5.9(a). Theaveragepressureremainsalmostconstantinoverburdenandbasementbecause theyarehydraulicallyisolatedfromthereservoir. Smallchangesinoverburdenand 119 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 0 200 400 600 800 1000 Qo+Qw (stb/day) 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 0 500 1000 1500 2000 2500 Qw (stb/d) 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 0 1000 2000 3000 4000 5000 Qg (scf/d) (a) (b) (c) Figure 5.8: Typical production and injection rate profiles of FWU wells. (a) The liquid production rate (Qo+Qw, oil+water) of a producer, (b) the water injection rate Qw of a water injector, and (c) the gas injection rate Qg of a gas injector. Production starts at the beginning of the simulation. Water injection starts in year 1995. Gas injection starts in year 2010. Note the variability and discontinuity in the rate profiles, which will be related to fault stability patterns in the analysis below. basement pressures do occur due to undrained deformation caused by changes in the total stress due to bulk volume changes in the reservoir. The reservoir pressure decreases at the beginning of the simulation due to production. Around year 1966, injection begins which leads to a temporary increase in the pressure. Since the total production flow rate is higher than the total injection flow rate, the reservoir pressure decreases thereafter. During 1995-2010, the production rate is curtailed for most of the producers while water injection continues at relatively high rates causing the average pressure to increase. CO 2 injection in year 2013 causes an increase in both production and pressure. The maximum pressure drop in the reservoir during the entire simulation is approximately 200 psi (1 psi = 6894.76 Pa). We select four representative time steps–day 4500 (year 1969), day 5000 (year 1970), day 13250 (year 1996) and day 22250 (year 2016)–to further investigate the result. We show the evolution of vertical displacement at the reservoir boundaries in Figure 5.9(b). The top boundary of the reservoir, Morrow-B Top, experiences downwardandupwardmovements, i.e. negativeandpositiveverticaldisplacement, 120 (b) 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 1940 1990 2040 2090 2140 Average reservoir pressure (psi) 0 1.5 3 4.5 6 Cumulative liquid production (STB) 10 7 22250 day 13250 day 5000 day t = 4500 day (a) Figure5.9: (a)Timeevolutionoftheaveragereservoirpressure(psi)duetoproduc- tion, injection and flow-induced stresses. The maximum pressure drop is approxi- mately 200 psi. Cumulative liquid production is shown by the red dot curve. (b) Time evolution of the vertical displacement at two selected points on Morrow-B Top (blue line) and Morrow-B Base (green dash line). Vertical displacement on the two surfaces have opposite signs because of their opposite motion in response to volumetric contraction/expansion of the reservoir with pressure decrease/increase. Displacement at the reservoir bottom has smaller magnitude due to the zero dis- placement boundary condition at the bottom boundary of the model. due to decreasing and increasing reservoir pressure, respectively. Morrow-B Base experiences movements in the opposite direction for the same pressure changes. The displacement magnitude is smaller at Morrow-B Base because of the zero displacement boundary condition at the bottom boundary of the model. We plot 3D fields of changes in the pressure and oil saturation for the reservoir region, i.e. within the depth interval of 1380-1410 m (Figure 5.10 and Figure 5.11). Att=0 day, the changes in pressure and oil saturation with respect to their initial values are zero. For 0<t<4500 day, production leads to a drop in pressure within the reservoir. Small amount of injection happens in the northeast part of FWU causing an increase in pressure and a decrease in the oil saturation around those injection wells. At t=5000 day, the increase in pressure and the decrease in oil saturation in the northeast part are larger while the pressure continues to decrease 121 t = 4500 day t = 5000 day t = 13250 day t = 22250 day (a) (b) (c) (d) -100 -200 0 100 200 (psi) Figure 5.10: Changes in the reservoir pressure compared to the initial pressure at the selected time steps. Pressure decreases in regions with net production and increases in regions with net injection. The maximum pressure drop is 200 psi approximately and occurs around t=13250 day. The final pressure drop at t=22250 day is 150 psi approximately. in other parts of the reservoir. At t=13250 day, a pressure drop of 200 psi can be observed in the northeast part, while it is approximately 100-150 psi in the remainder of the reservoir. As more wells are brought on-line, the oil saturation continues to drop around those wells. At the end of the simulation, t=22250 day, the average pressure drop in the reservoir is approximately uniform at 150 psi, and the average drop in the oil saturation around the active wells is 0.5 approximately. We want to identify the geomechanical signatures of reservoir compartmental- ization and bypassed oil reserves due to the presence of faults, sub-optimal well configurations, and various EOR processes such as high permeability channeling, viscous fingering, or gravity segregation of CO 2 (Lake et al., 2015). We plot 3D fields of the vertical displacement for the whole model over a depth interval of 1000-2000 m at the selected time steps (Figure 5.12). At t=0 day, the model is at mechanical equilibrium and all displacements are zero. Due to production-induced 122 t = 4500 day t = 5000 day t = 13250 day t = 22250 day (a) (b) (c) (d) Figure 5.11: Changes in the oil saturation with respect to the initial oil saturation at the selected time steps. The oil saturation decreases around both producers and injectors, as expected. The maximum change in the saturation is approximately 0.5. The oil saturation remains almost unchanged in regions far from the wells, which suggests the presence of bypassed oil. pressure depletion, the vertical displacement becomes negative (cause for subsi- dence) near the top of the model around 4500 day. At 5000 day, due to injection- induced pressure increase, the vertical displacement becomes positive (cause for uplift) in the northeastern part of the model where the injectors are located. At t=13250 day, negative vertical displacements can be seen above the reservoir while positive vertical displacements are seen below the reservoir. This is consistent with a volumetric contraction of the reservoir layers (Segall & Fitzgerald, 1998). At the end, the model shows a maximum subsidence of approximately -0.05 cm at the top boundary. The non-trivial shape of the displacement contours reflects the time-dependent and spatially variable flow rates of the FWU wells. 5.4.2 Fault results We choose Fault-1 as a representative fault to show the analysis steps. Fault-1 is a N-S striking fault located in the middle of the model within y=4000 m and 123 t = 4500 day t = 5000 day t = 13250 day t = 22250 day (a) (b) (c) (d) (cm) 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 Figure 5.12: Vertical displacement field (cm) at the selected time steps. Pressure drop in the reservoir leads to negative (blue color) vertical displacement above the reservoir while pressure increase in the reservoir leads to positive (red color) verti- caldisplacementabovethereservoir. Belowthereservoir, theverticaldisplacement shows opposite behavior as it experiences uplift during reservoir contraction and downward push during reservoir expansion. The asymmetry in the mechanical boundary conditions between the top (fixed traction) and bottom (fixed displace- ment) boundaries of the model affects the magnitude of vertical displacement. y=10000 m (Figure 5.3). The fault cuts through the model in the z-direction and is almost vertical. Figure 5.14(h) and Figure 5.15(g) show the wells within a distance of 1000 m from Fault-1. There are 18 such wells, with 12 producers and 6 injectors, which started operating during the first 5000 days of simulation. They are numbered in Figure 5.14(h) and Figure 5.15(g). Most of the producers started earlier while the injectors started later at approximately t=4000 day. Given the importance of the fault pressure p f in determining fault stability (Equation 5.1), we analyze the evolution of fault pressure on Fault-1 across the selected time steps in Figure 5.13. Since there are wells on both sides of Fault-1 causing pressure to change on both sides, the fault pressure within the reservoir depth interval shows a trend with time that is similar to that of the reservoir pressure in Figure 5.10. 124 -100 -50 0 50 100 t = 4500 day t = 5000 day t = 13250 day t = 22250 day (a) (b) (c) (d) y z Reservoir (psi) Figure 5.13: Change in the fault pressure (psi) along Fault-1 at the selected time steps. A trend similar to that of Figure 5.10 can be observed. The fault pressure in the reservoir depth interval (appearing as a horizontal stripe) declines at early times due to production from both fault blocks. Later, it increases due to injection and then decreases again towards the end of the simulation due to net production. Next we analyze the change in tractions on Fault-1. We consider traction changes at t = 4500 and 5000 days to distinguish between production-dominated and injection-dominated effects. There are three tractions on a fault surface: lat- eral or along-strike shear, along-dip shear, and effective normal. The induced lateral shear on Fault-1 is negligible compared to the induced dip shear and effec- tive normal tractions. This is primarily a result of the well locations and flow rates. Reservoir contraction- or expansion-induced lateral shear around a producer or an injector cancel out when the effects are superposed over multiple wells distributed relatively uniformly in the fault blocks, which is the situation in FWU. In con- trast, the dip shear and normal tractions induced by net production or injection from all the wells concentrated within a relatively thin reservoir interval are not negligible. For example, inside the reservoir, production will induce compression near the producer and extension away from the producer (Segall & Fitzgerald, 1998). Therefore, production leads to positive (tensile) normal traction changes inside the reservoir and negative (compressive) normal traction changes in the overburden and basement (Figure 5.14(a)). Along the dip direction, production 125 - 10000 7000 5000 6000 8000 9000 Well 1 Well 7 Well 8 Well 2 Well 3 Well 4 Well 9 Well 6 Well 11 Well 5 Well 10 Well 12 -0.15 -0.075 0.00 0.075 0.15 -0.5 -0.2 0.00 0.2 0.5 y z 300 0 (a) (c) (d) (e) (f) (g) 4500 1 2 3 4 5 6 7 8 9 10 11 12 (b) (MPa) (MPa) Figure 5.14: Effect of spatial distribution and temporal evolution of production on fault stability. Changes in (a) the effective normal traction and (b) the updip shear tractionon Fault-1 att=4500 daywith respectto theinitial traction values. Liquid production profiles of active producers near Fault-1 are shown in (c) through (g). Surface locations of the wells with respect to the fault are shown in (h). will induce downdip shear in overburden and updip shear in the basement (Fig- ure 5.14(b)). The induced traction directions are opposite in case of injection as shown in Figure 5.15. We calculate the change in Coulomb Failure Function (ΔCFF) on the faults. Positive or negative ΔCFF implies a decrease or an increase in the fault stability, respectively. The results for Fault-1 are shown in Figure 5.16. Let’s consider the reservoir domain. At the end of primary depletion of FWU, at t=4500 day, 126 production-induced contraction causes a tensile change in the total normal traction (Δσ n > 0)andanincreaseinthesheartractionmagnitude(Δ|τ|> 0)immediately adjacent to the reservoir. The fault pressure decreases inside the reservoir as shown inFigure5.13. Sincethemagnitudeofthepressurechangeislargerthanthatofthe total traction change, the effective normal traction decreases. Also, the effective normal traction changes are larger than the shear traction changes given the fault dipangle. Therefore, ΔCFFisnegative(stabilizing)att=4500day. Att=5000day, when injection dominates, the fault pressure increases and overcomes the decrease in the total normal traction due to injection-induced expansion. Therefore, the effective normal traction increases and ΔCFF is positive (destabilizing) att=5000 day. As production continues and the fault pressure further decreases, results similar to t=4500 day can be observed at t=13250 day. For t>13250 day, because of water and CO 2 injection, the fault pressure starts to increase but remains less than its initial value (Δp f <0). Therefore, an increase in ΔCFF can be observed, although it remains negative. To understand the difference in evolution of mechanical stability within the reservoir and the overburden, we plot ΔCFF at selected points on the faults in the reservoir and the overburden interval (Figure 5.17). As expected, the variations are larger within the reservoir interval with a range of -0.4 MPa to 0.1 MPa. The ΔCFF range in the overburden is much smaller: -0.017 to 0.0025 MPa. Positive (destabilizing) changes in the Coulomb stress occurs around year 1970 after water injection begins. Then it decreases and becomes negative due to production and pressure depletion. From year 1980 to 1985, the cumulative injection volume is higher than the cumulative production volume. Thus the second peak of ΔCFF appears around year 1985 (Figure 5.17(c)). During 1995-2010, field development slowed down causing a decrease in the field production and a build-up in the field 127 pressure and ΔCFF. Around year 2010, CO 2 injection began and ΔCFF increases accordingly. However, it remains negative. To explain the evolution of ΔCFF and how it relates to pressure evolution, we plot ΔCFFandtheelementpressuretogetherinFigure5.18(a)to(c)forthepoints in the reservoir interval. For the selected point shown in the inset figure, we choose two elements in the reservoir interval from the two sides of the fault, i.e. from the hanging wall and footwall blocks of the fault. Pressures from the two elements are plotted along with ΔCFF at the selected point. ΔCFF shows a general trend similar to that of the two pressures, which is expected because well operations are driving the changes in mechanical stability of the faults. However, there are two non-trivial observations. First, ΔCFF is not following either of the two pressures for the entire duration. Instead, it follows the larger of the two pressures because that defines the fault pressurep f . The fault block with the larger pressure switches between hanging wall and footwall blocks as time evolves due to the variation in well operation with time. Second observation is that during 2011-2016 ΔCFF on Fault-1 decreases while the two pressures are increasing or staying constant. In other words, ΔCFF does not seem to follow the larger pressure during this period, which corresponds to the period of CO 2 injection and tertiary recovery in FWU. The reason is that the change in shear traction dominate ΔCFF during this period. Injection in the top reservoir layer applies downdip shear on points along the boundary between the top and bottom reservoir layers. Since we are interested in extracting characteristic behaviors of fault destabi- lization induced by well operations, we plot ΔCFF as a function of the ‘driving force’ i.e. the reservoir pressure drop. Let Δp rd =|1− (p res /p res,0 )| be the magni- tude of dimensionless change in reservoir pressure due to production or injection. Depending on the model discretization size and inter-well spacing, this could be 128 defined by averaging the pressure change over multiple grid elements in the reser- voir. We plot ΔCFF versus Δp rd for the three faults in Figure 5.18(d) to (f). The characteristic trend is a straight line with a negative slope. The deviations from this trend correspond to the start of water and CO 2 injection (Figure 5.18(d)) and major changes in well operations around the fault (Figure 5.18(e)). During production-dominated period, Δp rd increases and ΔCFF decreases. Hence, points evolving downward along the trend line suggests stabilization on this plot. During injection, the reservoir pressure increases such that the pressure drop compared to the initial pressure, Δp rd , decreases. ΔCFF may increase or decrease depending on the relative magnitudes of fault tractions and pressure (Eq. (5.1)). Points evolving upward along the trend line, i.e. increasing ΔCFF, suggests destabilization. How- ever, injection does not always translate into destabilization as shown by the post- CO 2 -injection gray dots evolving horizontally towards left in Figure 5.18(d), which suggests an approximately constant ΔCFF and neutral stability during injection. A similar analysis for Fault-3L and Fault-3R shows that injection is increasing the risk of induced destabilization on these two faults. Moreover, the destabilization features on the two faults are different (gray dots jump down for Fault-3L and up for Fault-3R) due to different locations of the two faults with respect to the wells (refer to Figure 5.11(c), (d)). 5.5 Discussion and Conclusion We presented results from a coupled flow-geomechanical simulation study of CO 2 storage and enhanced oil recovery (SEOR) in the Farnsworth Unit oilfield of Texas. We demonstrated a novel workflow to incorporate complex geological structure, petrophysical heterogeneity, and multi-well schedule of a real oilfield 129 producing under different recovery mechanisms–primary depletion, waterflooding, and CO 2 flooding–into a rigorous multiphase flow-geomechanics simulation frame- work. We applied the framework to the FWU dataset to predict the geomechanical effects of fluid transport and storage on faults. The main highlights are as follows: 1. Our workflow (Figures 5.5 through 2.3) demonstrates how to use existing uncoupled models to initialize and distribute properties and wells in a cou- pled model. This is important because many oilfields have uncoupled, flow- only simulation models that are calibrated continually with significant effort because of their role in reporting hydrocarbon reserves and finding new drilling locations in the oilfield. Finding an efficient workflow to convert them into coupled models can accelerate the adoption of the coupled model- ing technology. 2. Our workflow shows how to identify and characterize geomechanical signa- tures of flow compartmentalization and bypassed oil due to faults, well con- figurations, and well rate variations. The simulated horizontal and vertical displacement fields (e.g. Figure 5.12) can be used to determine whether the reservoir regions, left unswept during CO 2 or water-flooding, can or cannot be detected with the aid of surface deformation data e.g. InSAR and GPS. 3. We demonstrate that we can maintain an approximately constant ΔCFF (Figure 5.18(d)) on critically-oriented faults, and thereby control their induced seismicity risks, by carefully managing the wells. 130 4000 10000 7000 5000 6000 8000 9000 -0.5 -0.2 0.00 0.2 0.5 (MPa) -0.15 -0.075 0.00 0.075 0.15 y z (MPa) -0.5 -0.2 0.00 0.2 0.5 (MPa) (a) (b) (c) Figure 5.15: Effect of spatial distribution and temporal evolution of injection on fault stability. Changes in (a) the effective normal traction (MPa) and (b) the updip shear traction (MPa) on Fault-1 at t=5000 day. Water injection in six wells, shown in (c) through (f) along the fault strike, causes updip shear above the reservoir(redcolorin(b))anddowndipshearbelowthereservoir(bluecolorin(b)). Tensile stresses are induced in the overburden and basement, which cause tensile changes in the effective normal traction (red color in (a)). Production dominates injection in the south, which causes compressive changes in the effective normal traction (blue color in (a)) near well 2. Surface locations of the wells with respect to the fault are shown in (g). 131 y z t = 4500 day t = 5000 day t = 13250 day t = 22250 day (a) (b) (c) (d) -0.50 -0.25 0.00 0.25 0.50 (MPa) Figure5.16: ChangesintheCoulombfailurestress(MPa)onFault-1attheselected timesteps. Redcolororpositivevaluesindicateinduceddestabilization. Bluecolor or negative values indicate induced stabilization. (a) fault 3L fault 3R fault 1 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 -10 -5 0 5 CFF (MPa) 10 -3 CFF (overburden) fault 1 fault 3L fault 3R Figure 5.17: Evolution of ΔCFF on selected points within the overburden (b) and reservoir (c) intervals on Fault-1, Fault-3L and Fault-3R. Panel (a) shows the locations of the points in the reservoir interval (blue dots) and the overburden interval (red dots). 132 (a) (b) (c) 0 0.05 0.1 -0.6 -0.4 -0.2 0 0.2 0 0.05 0.1 0 0.05 0.1 0 1 2 3 10 4 (d) (e) (f) 196 1900 1935 1970 2005 2040 2075 2110 2145 2180 Pressure (psi) 1960 1970 1980 1990 2000 2010 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 dCFF (MPa) 2180 2110 2040 1970 1900 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 dCFF (MPa) 1900 1935 1970 2005 2040 2075 2110 2145 2180 Pressure (psi) Pressure (psi) Water injection CO2 injection 4500 day 5000 day 22250 day 13250 day Figure 5.18: Comparison of ΔCFF (solid) and reservoir pressure (dotted and dashed) evolutions on (a) Fault-1, (b) Fault-3L and (c) Fault-3R at points within the reservoir interval. The two pressure curves are from two elements connected to the blue dot (inset) on the two sides of the fault. ΔCFF follows the trend of the larger pressure because p f is defined from the larger of the two pressures. Plots (d), (e) and (f) show the characteristic ΔCFF versus Δp rd behavior for Fault-1, Fault-3L and Fault-3R, respectively. Time evolution is shown through the color of the dots. The color legend on the right shows that the time evolves from t=0 (black) to t=22250 day (gray). Start of water injection and CO 2 injection are marked in (d). 133 Reference List Aagaard B, Kientz S, Knepley M, Strand L, Williams C (2012) PyLith User Man- ual, Version 1.8.0 Computational Infrastructure for Geodynamics, University of California, Davis. AagaardBT,KnepleyMG,WilliamsCA(2013) Adomaindecompositionapproach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation. J. Geophys. Res. 118:3059–3079. Albano M, Barba S, Tarabusi G, Saroli M, Stramondo S (2017a) Discriminat- ing between natural and anthropogenic earthquakes: insights from the Emilia Romagna (Italy) 2012 seismic sequence. Nat. Sci. Reports . Albano M, Barba S, Tarabusi G, Saroli M, Stramondo S (2017b) Discriminat- ing between natural and anthropogenic earthquakes: insights from the Emilia Romagna (Italy) 2012 seismic sequence. SCI. REP-UK. . Allen DR (1968) Physical changes of reservoir properties caused by subsidence and repressurizing operations, Wilmington field, California. J. Pet. Tech- nol. 20:23–29. Allen DR, Mayuga MN (1969) The mechanics of compaction and rebound, Wilm- ington Oilfield, long beach, ca, usa In Land Subsidence, Vol. 89, pp. 410–413. IASH-UNESCO, Tokyo. Ampomah W, Balch R, Cather M, Rose-Coss D, Dai Z, Heath J, Dewers T, Mozley P (2016) Evaluation of CO2 storage mechanisms in CO2 enhanced oil recovery sites: Application to Morrow Sandstone Reservoir. Energy & Fuels 30:8545–8555. Ampomah W, Balch RS, Cather M, Rose-Coss D, Gragg E (2017a) Numerical sim- ulation of CO2-EOR and storage potential in the Morrow Formation, Ochiltree County, Texas In SPE Oklahoma City Oil and Gas Symposium, number SPE- 185086-MS, Oklahoma City, Oklahoma. 134 AmpomahW,BalchRS,CatherM,WillR,GundaD,DaiZ,SoltanianMR(2017b) Optimum design of CO2 storage and oil recovery under geological uncertainty. Appl. Energ. 195:80–92. Ampomah W, Balch RS, Ross-Coss D, Hutton A, Cather M, Will RA (2016) An integrated approach for characterizing a sandstone reservoir in the Anadarko Basin. Offshore Technology Conference . ArmeroF(1999) Formulationandfiniteelementimplementationofamultiplicative modelofcoupledporo-plasticityatfinitestrainsunderfullysaturatedconditions. Comput. Meth. Appl. Mech. Eng. 171:205–241. Armero F, Simo JC (1993) A priori stability estimates and unconditionally sta- ble product formula algorithms for nonlinear coupled thermoplasticity. Int. J. Plasticity 9:749–782. Aziz K, Settari A (1979) Petroleum Reservoir Simulation Elsevier, London. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, Eijkhout V, Gropp WD, Karpeyev D, Kaushik D, Knepley MG, May DA, McInnes LC, Mills RT, Munson T, Rupp K, Sanan P, Smith BF, Zampini S, Zhang H (2019) PETSc Web page https://www.mcs.anl.gov/petsc. Balch R, McPherson B, Grigg R (2017) Overview of a large scale carbon capture, utilization, and storage demonstration project in an active oil field in texas, usa. Enrgy. Proced. 114:5874 – 5887. Bathe KJ (1995) Finite Element Procedures Prentice Hall, New Jersey. Benallal A, Botta AS, Venturini WS (2008) Consolidation of elastic–plastic sat- urated porous media by the boundary element method. Comput. Meth. Appl. Mech. Eng. 197:4626–4644. Bennett K, Berla L, Nix W, Borja R (2015) Instrumented nanoindentation and 3d mechanistic modeling of a shale at multiple scales. Acta Geotech 10:1–14. Bissell R, Vasco D, Atbi M, Hamdani M, Okwelegbe M, Goldwater M (2011) A full field simulation of the in Salah gas production and CO2 storage project using a coupled geo–mechanical and thermal fluid flow simulator. Enrgy. Proced. 4:3290–3297. Borja RI (2002) Finite element simulation of strain localization with large defor- mation: capturing strong discontinuity using a Petrov–Galerkin multiscale for- mulation. Comput. Meth. Appl. Mech. Eng. 191:2949–2978. 135 Borja RI (2004) Cam-Clay plasticity. Part V: A mathematical framework for three- phase deformation and strain localization analyses of partially saturated porous media. Comput. Meth. Appl. Mech. Eng. 193:5301–5338. Borja RI (2006) On the mechanical energy and effective stress in saturated and unsaturated porous continua. Int. J. Solids and Struct. 43:1764–1786. Borja RI (2013) Plasticity Modeling & Computation Springer, New York. Borja RI, Alarcon E (1995) A mathematical framework for finite strain elastoplas- ticconsolidationPart1: Balancelaws, variationalformulation, andlinearization. Comput. Meth. Appl. Mech. Eng. 122:145–171. Bourdet D (2002) Well Test Analysis: The use of Advanced Interpretation Models Elsevier, Amsterdam. Bourne SJ, Oates SJ (2017) Extreme threshold failures within a heterogeneous elastic thin sheet and the spatial-temporal development of induced seismicity within the groningen gas field. jgr 122:10,299–10,320. Bubshait A, Jha B (2019) Coupled poromechanics-damage mechanics model- ing of fracturing during injection in brittle rocks. Int. J. Numer. Methods Eng. pp. 1–21. Buijze L, van den Bogert P, Wassing B, Orlic B (2019) Nucleation and arrest of dynamic rupture induced by reservoir depletion. J. Geophys. Res. 124. Callari C, Armero F (2004) Analysis and numerical simulation of strong discontinuities in finite strain poroplasticity. Comput. Meth. Appl. Mech. Eng. 193:2941–2986. Cappa F, Rutqvist J (2011a) Impact of CO 2 geological sequestration on the nucle- ation of earthquakes. Geophys. Res. Lett. 38:L17313. Cappa F, Rutqvist J (2011b) Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO 2 . Int. J. Greenh. Gas Control 5:336–346. Cappa F, Rutqvist J (2012) Seismic rupture and ground accelerations induced by CO 2 injection in the shallow crust. Geophys. J. Int. 190:1784–1789. Carter JP, Booker JR, Small JC (1979) An improved method for calculating water influx. Int. J. Numer. Anal. Methods Geomech. 3:107–129. 136 Castelletto N, Teatini P, Gambolati G, Bossie-Codreanu D, Vincké O, Daniel JM, Battistelli A, Marcolini M, Donda F, Volpi V (2013) Multiphysics mod- eling of CO2 sequestration in a faulted saline formation in Italy. Adv. Water Resour. 62:570–587. Castelletto N, White JA, Tchelepi HA (2015) Accuracy and convergence properties of the fixed-stress iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Methods Geomech. . Castineira D, Jha B, Juanes R (2015) Uncertainty quantification and inverse modeling of fault poromechanics and induced seismicity: Application to a syn- thetic carbon capture and storage (CCS) problem In 50th US Rock Mechan- ics/Geomechanics Symposium, Am. rock Mech. Assoc., number ARMA-2016- 151. Catalli F, Meier MA, Wiemer S (2013) The role of Coulomb stress changes for injection-induced seismicity: The Basel enhanced geothermal system. Geophys. Res. Lett. 40:72—-77. CC (2019) Creative Commons – Attribution 4.0 International – CC BY 4.0 https: //creativecommons.org/licenses/by/4.0/legalcode Accessed: 2019-05-16. Cesca S, Dahm T, Juretzek C, Kuhn D (2011) Rupture process of the 2001 May 7 Mw 4.3 Ekofisk induced earthquake. Geophys. J. Int. 187:407–413. Chan AW, Zoback MD (2007) The Role of Hydrocarbon Production on Land Subsidence and Fault Reactivation in the Louisiana Coastal Zone. Journal of Coastal Research 23:771–786. Chang C, Mallman E, Zoback MD (2014) Time-dependent subsidence associated with drainage-induced compaction in Gulf of Mexico shales bounding a severely depleted gas reservoir. AAPG Bulletin 98:1145–1159. Chang KW, Segall P (2016a) Injection induced seismicity on basement faults including poroelastic stressing. J. Geophys. Res. 121:2708—-2726. Chang KW, Segall P (2016b) Seismicity on basement faults induced by simulta- neous fluid injection–extraction. Pure App. Geophys. . Chiaramonte L, Zoback MD, Friedmann J, Stamp V (2008) Seal integrity and feasibility of CO2 sequestration in the Teapot Dome EOR pilot: geomechanical site characterization. Environ. Geol. 54:1667–1675. Chin LY, Thomas LK (1998) Fully-coupled geomechanics and fluid-flow analysis of wells with stress-dependent permeability In SPE International Conference and Exhibition, p. SPE 48857, Beijing, China. 137 Coussy O (1995) Mechanics of Porous Continua John Wiley and Sons, Chichester, England. Coussy O (2004) Poromechanics John Wiley and Sons, Chichester, England. CSimSoft (2016) Trelis 16.1 User Documentation CSimSoft, LLC. DahmT,KrugerF,StammlerK,KlingeK,KindR,WylegallaK,GrassoJR(2007) The 2004 Mw 4.4 Rotenburg, Northern Germany, Earthquake and its possible relationship with gas recovery. Bull. Seismol. Soc. Am. 97:691–704. Davis RO, Selvadurai APS (2002) PLASTICITY AND GEOMECHANICS Cam- bridge University Press, New York. Deng K, Liu Y, Harrington RM (2016) Poroelastic stress triggering of the Decem- ber 2013 Crooked Lake, Alberta, induced seismicity sequence. Geophys. Res. Lett. 43:8482–8491. Dieterich JH (1994) A constitutive law for rate of earthquake production and its application to earthquake clustering. J. Geophys. Res. 99:2601–2618. El-kaseeh G, Will R, Balch R, Grigg R (2017) Multi–scale seismic measurements for CO2 monitoring in an EOR/CCUS project. Enrgy. Proced. 114:3656–3670. Ellsworth WL (2013) Injection-induced earthquake. Science pp. 1–7. Eshelby JD (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. A 241:376–396. Evans JL (1979) Major structural and stratigraphic features of the anadarko basin In Hyne NJ, editor, Pennsylvanian Sandstones of the Mid-Continent, pp. 97–113. Tulsa Geological Society Special Publication. Evans KF, Zappone A, Kraft T, Deichmann N, Moia F (2012) A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe. Geothermics 41:30–54. Fan Z, Eichhubl P, Gale JFW (2016) Geomechanical analysis of fluid injection and seismic fault slip for the Mw4.8 Timpson, Texas, earthquake sequence. J. Geophys. Res. 121:2798–2812. Ferris JG, Knowles DB, Brown RH, Stallman RW (1962) Theory of Aquifer Tests US Government Printing Office, Washington. F.G. Bell JC MC (1986) A review of the engineering behaviour of soils and rocks withrespecttogroundwater. Geochem. Soc. Lond. Eng. Geol. Spec. Publ.3:1–23. 138 Franklin OM (2009) Onshore geologic storage of CO2. Science 325:1656–1658. Fredrich JT, Arguello JG, Deitrick GL, de Rouffignac EP (2000) Geomechanical modeling of reservoir compaction, surface subsidence, and casing damage at the Belridge diatomite field. SPE Reserv Eval Eng 3:348–359. Galis M, Ampuero JP, Mai PM, Cappa F (2017) Induced seismicity provides insight into why earthquake ruptures stop. Science Advances 3. Galis M, Ampuero J, Mai P, Kristek J (2019) Initiation and arrest of earthquake ruptures due to elongated overstressed regions. Geophys. J. Int. 217. Galloway DL, Riley FS (1999) San Joaquin Valley, California - Largest human alteration of the Earth’s surface In Galloway DL, Jones DR, Ingebritsen SE, editors, Land subsidence in the United States, pp. 23–34. U.S. Geological Survey Circular. Gan W, Frohlich C (2013) Gas injection may have triggered earthquakes in the Cogdell oil field, Texas. Proc. Natl. Acad. Sci. U.S.A. 110:18786––18791. Garagash D, Germanovich L (2012) Nucleation and arrest of dynamic slip on a pressurized fault. J. Geophys. Res. 117. Garipov TT, Karimi-Fard M, Tchelepi HA (2016) Discrete fracture model for coupled flow and geomechanics. Comput. Geosci. 20:149–160. Gonzalez PJ, Tiampo KF, Palano M, Cannavo F, Fernandez J (2012) The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading. Nature Geosci. 5:821–825. Gu XY, Ran QQ (2000) Land subsidence In Proc. 6th Int. Symposium on Land Subsidence, Vol. 11, pp. 355–365. Guglielmi Y, Cappa F, Avouac J, Henry P, Elsworth D (2015) Seismicity triggered by fluid injection–induced aseismic slip. Science 348. Hanks TC, Kanamori H (1979) A moment magnitude scale. J. Geophys. Res. 84:2348–2350. Hansen O, Gilding D, Nazarian B, Osdal B, Ringrose P, Kristoffersen JB (2013) Snøhvit: The history of injecting and storing 1 Mt CO 2 in the fluvial Tubaen Fm. Energy Procedia 37:3565–3573. Harris RA (1998) Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard. J. Geophys. Res. 103:24347–24358. 139 Hearn EH, Burgmann R (2005) The effect of elastic layering on inversions of GPS dataforcoseismicslipandresultingstresschanges: Strike-slipearthquakes. Bull. Seismo. Soc. Am. 95:1637—-1653. Heimisson E (2019) Constitutive law for earthquake production based on rate- and-state friction: Theory and application of interacting sources. J. Geophys. Res. 124. Hill R (1965) Continuum micro-mechanics of elastoplastic polycrystals. J. Mech. Phys. Solids 13:89–101. Hornbach MJ, Deshon HR, Ellsworth WL, Stump BW, Hayward C, Frohlich C, Oldham HR, Olson JE, Magnani MB, Brokaw C, Luetgert JH (2015a) Causal factors for seismicity near Azle, Texas. Nature Comm. 6:1067–1074. Hornbach MJ, DeShon HR, Ellsworth WL, Stump BW, Hayward C, Frohlich C, Oldham HR, Olson JE, Magnani MB, Brokaw C, Luetgert JH (2015b) Causal factors for seismicity near Azle, Texas. Nat. Comm. 6. Horton S (2012) Disposal of hydrofracking waste fluid by injection into subsur- face aquifers triggers earthquake swarm in central Arkansas with potential for damaging earthquake. Seismol. Res. Lett. 83:250–260. Hosseini SM, Goebel THW, Jha B, Aminzadeh F (2018) A probabilistic approach to injection-induced seismicity assessment in the presence and absence of flow boundaries. Geophys. Res. Lett. 45:8182–8189. Hutton A (2015) Geophysical modeling and structural interpretation of a 3d reflec- tion seismic survey in Farnsworth Unit, TX M.S. thesis, New Mexico Institute of Mining and Technology, Department of Geophysics. Ide S, Takeo M (1997) Determination of constitutive relations of fault slip based on seismic wave analysis. J. Geophys. Res. 102:27379–27391. IrelandRL,PolandJF,RileyFS(1984) LandsubsidenceintheSanJoaquinValley, California as of 1980. US Geol. Surv. Prof. Pap. 93:437–I. Jagalur-Mohan J, Jha B, Wang Z, Juanes R, Marzouk Y (2018) Inferring fault frictional and reservoir hydraulic properties from injection-induced seismicity. Geophys. Res. Lett. 45. Jha B, Bottazzi F, Wojcik R, Coccia M, Bechor N, Mclaughlin D, Herring TA, Hager BH, Mantica S, Juanes R (2015) Reservoir characterization in an under- ground gas storage field using joint inversion of flow and geodetic data. Int. J. Numer. Anal. Methods Geomech. 39:1619–1638. 140 Jha B, Juanes R (2014) Coupled multiphase flow and poromechanics: A compu- tational model of pore pressure effects on fault slip and earthquake triggering. Water Resour. Res. 50:3776–3808. Juanes R, Jha B, Hager BH, Shaw JH, Plesch A, Astiz L, Dieterich JH, Frohlich C (2016) Were the may 2012 Emilia-Romagna earthquakes induced? A coupled flow-geomechanics modeling assessment. Geophys. Res. Lett. 43. Keranen KM, Savage HM, Abers GA, Cochran ES (2013) Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 M w 5.7 earthquake sequence. Geology 41:699–702. Kim J, Tchelepi HA, Juanes R (2011) Stability, accuracy and efficiency of sequen- tial methods for coupled flow and geomechanics. Soc. Pet. Eng. J. 16:249–262. King GCP, Stein RS, Lin J (1994) Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am. 84:935–953. Lackner KS (2003) A guide to CO2 sequestration. Science 300:1677–1678. Lake LW, Johns R, Rossen B, Pope G (2015) Fundamentals of Enhanced Oil Recovery Society of Petroleum Engineers, Texas. Lei X, Huang D, Su J, Jiang G, Wang X, Wang H, Guo X, Fu H (2017) Fault reactivation and earthquakes with magnitudes of up to Mw4.7 induced by shale- gas hydraulic fracturing in Sichuan Basin, China. Nat. Sci. Reports 7. Lewis RW, Schrefler BA (1998) The Finite Element Method in the Static and Dynamic deformation and Consolidation of Porous Media Wiley, Chichester, England, second edition. Li X, Liu Z, Lewis RW (2005) Mixed finite element method for coupled thermo- hydro-mechanical process in poro-elastic-plastic media at large strains. Int. J. Numer. Meth. Eng. 64:667–708. Lie KA (2016) An Introduction to Reservoir Simulation Using MATLAB: User guide for the Matlab Reservoir Simulation Toolbox (MRST) SINTEF ICT. Ma X, Elbanna AE (2015) Effect of off-fault low-velocity elastic inclusions on supershear rupture dynamics. Geophys. J. Int. 203:664–677. MainguyM,LonguemareP(2002) Couplingfluidflowandrockmechanics: Formu- lations of the partial coupling between reservoir and geomechanical simulators. Oil Gas SCI. Technol. 57:355–367. 141 Marone C (1998) Laboratory-derived friction laws and their application to seismic faulting. Ann. Rev. Earth Planet. Sci. 26:643–696. Maury V, Piau JM, Halle G (1996) Subsidence induced by water injection in water sensitive reservoir rocks: The example of ekofisk In European Petroleum Conference, p. SPE 36890, Milan, Italy. Mazzoldi A, Rinaldi AP, Borgia A, Rutqvist J (2012) Induced seismicity within geologic carbon sequestration projects: maximum earthquake magnitude and leakage potential. Int. J. Greenh. Gas Control 10:434–442. McGarr A (1999) On relating apparent stress to the stress causing earthquake fault slip. J. Geophys. Res. 104:3003–3011. McGarr A (2014) Maximum magnitude earthquakes induced by fluid injection. J. Geophys. Res. 119:1008—-1019. Mignan A, Broccardo M, Wiemer S (2017) Induced seismicity closed-form traffic light system for actuarial decision-making during deep fluid injections. Sci. Rep. 7:13607. MikumoT,OlsenKB,FukuyamaE,YagiY(2003) Stress-breakdowntimeandslip- weakening distance inferred from slip-velocity functions on earthquake faults. Bull. Seismol. Soc. Am. 93:264–282. Minkoff SE, Stone CM, Bryant S, Peszynska M, Wheeler MF (2003a) Coupled fluid flow and geomechanical deformation modeling. J. Pet. Sci. Eng. 38:37–56. MinkoffSE,StoneCM,BryantS,PeszynskaM,WheelerMF(2003b) Coupledfluid flow and geomechanical deformation modeling. J. Pet. Sci. Eng. 38:37 – 56. Miocic JM, Gilfillan SMV, Frank N, Schroeder-Ritzrau A, Burnside NM, Haszel- dine RS (2019) 420,000 year assessment of fault leakage rates shows geological carbon storage is secure. Scientific Reports 9. Morris JP, Hao Y, Foxall W, McNab W (2011) A study of injection-induced mechanical deformation at the In Salah CO 2 storage project. Int. J. Greenh. Gas Control 5:270–280. Mura T (1982) Micromechanics of Defects in Solids Martinus Nijhoff Publishers, The Hague. Nikitin K, Terekhov K, Vassilevski Y (2014) A Monotone Nonlinear Finite Volume Method for Diffusion Equations and Multiphase Flows. Comp. Geosci. 18:311–324. 142 Olden P, Pickup G, Jin M, Mackay E, Hamilton S, Somerville J, Todd A (2012) Use of rock mechanics laboratory data in geomechanical modeling to increase confidence in CO2 geological storage. Int. J. Greenh. Gas Con. 11:304–315. Olsen KB, Madariaga R, Archuleta RJ (1997) Three-dimensional dynamic simu- lation of the 1992 Landers Earthquake. Science 278:834–838. Owen D, Hinton E (1980) Finite Elements in Plasticity: Theory and Practice Pineridge Press Limited, Swansea. Panda D, Kundu B, Gahalaut VK, Burgmann R, Jha B, Asaithambi R, Yadav RK, Vissa NK, Bansal AK (2018) Seasonal modulation of deep slow-slip and earthquakes on the main himalayan thrust. Nat. Comm. 9:4140. Parry RH (2005) Mohr Circles, Stress Paths and Geotechnics CRC Press, London, 2nd edition. Preisig M, Prévost JH (2011) Coupled multi–phase thermo–poromechanical effects. case study: CO2 injection at in Salah, Algeria. Int. J. Greenh. Gas Con. 5:1055–1064. Raleigh CB, Healy JH, Bredehoeft JD (1976) An experiment in earthquake control at Rangely, Colorado. Science 191:1230–1237. Reasenberg PA, Simpson RW (1992) Response of regional seismicity to the static stress change produced by the Loma Prieta earthquake. Science 255:1687–1690. Rice JR (1980) The mechanics of earthquake rupture. Rice JR (1993) Spatio-temporal complexity of slip on a fault. J. Geophys. Res. 98:9885–9907. Rubin A, Ampuero J (2012) Earthquake nucleation on rate and state faults – aging and slip laws. J. Geophys. Res. 113. Rutqvist J, Birkholzer J, Cappa F, Tsang CF (2007) Estimating maximum sus- tainable injection pressure during geological sequestration of CO2 using cou- pled fluid flow and geomechanical fault–slip analysis. Energ. Convers. Man- age. 48:1798–1807. Rutqvist J, Birkholzer J, Tsang CF (2008) Coupled reservoir—geomechanical anal- ysis of the potential for tensile and shear failure associated with CO2 injection in multilayered reservoir—caprock systems. Int. J. Rock Mech. Min. 45:132–143. Rutqvist J, Rinaldi AP, Cappa F, Moridis GJ (2013) Modeling of fault reactivation andinducedseismicityduringhydraulicfracturingofshale-gasreservoirs. J. Pet. Sci. Eng. 107:31–44. 143 Sanavia L, Schrefler BA, Steinmann P (2002) A formulation for an unsaturated porous medium undergoing large inelastic strains. Comput. Mech. 28:137–151. Sanz PF, Lele SP, Searles KH, Hsu SY, Garzon JL, Burdette JA, Kline WE, Dale BA, Hector PD (2015) Geomechanical analysis to evaluate production-induced fault reactivation at groningen gas field. society of petroleum engineers. Soc. Pet. Eng. p. SPE 174942. Sarris E, Papanastasiou P (2013) Numerical modeling of fluid–driven frac- tures in cohesive poroelastoplastic continuum. Int. J. Numer. Anal. Methods Geomech. 37:1822–1846. Schlumberger (2017) Eclipse 100 Reference Manual 2017.1 Houston. Schmitt DR, Currie CA, Zhang L (2012) Crustal stress determination from bore- holes and rock cores: Fundamental principles. Tectonophysics 580:1–26. Schoenball M, Baujard C, Kohl T, Dorbath L (2012) The role of triggering by static stress transfer during geothermal reservoir stimulation. J. Geophys. Res. 117:B09307. Scholz CH (1968) Microfracturing and the inelastic deformation of rock in com- pression. J. Geophys. Res. 73:1417–1432. Scholz CH (1988) The critical slip distance for seismic faulting. Nature 336:761–763. Scholz CH (1989) Mechanics of faulting. Ann. Rev. Earth Planet. Sci. 17:309–334. Scholz CH (2002) The Mechanics of Earthquakes and Faulting Cambridge Univer- sity Press, Cambridge, UK, 2nd edition. Segall P (1985) Stress and subsidence resulting from subsurface fluid with- drawal in the epicentral region of the 1983 Coalinga earthquake. J. Geophys. Res. 90:6801–6816. Segall P (1989a) Earthquakes triggered by fluid extraction. Geology 17:942–946. Segall P (1989b) Earthquakes triggered by fluid extraction. Geology 17:942. Segall P, Fitzgerald SD (1998) A note on induced stress changes in hydrocarbon and geothermal reservoirs. Tectonophys. 289:117–128. Segall P, Grasso JR, Mossop A (1994) Poroelastic stressing and induced seismicity near the Lacq gas field. J. Geophys. Res. 99:15423–15438. 144 Segall P, Lu S (2015a) Injection–induced seismicity: Poroelastic and earthquake nucleation effects. J. Geophys. Res. 120:5082–5103. Segall P, Lu S (2015b) Injection induced seismicity: Poroelastic and earthquake nucleation effects. J. Geophys. Res. 120:5082—-5103. Shi XQ, Xue YQ, Ye SJ, Wu JC, Zhang Y, Yu J (2007) Characterization of land subsidence induced by groundwater withdrawals in Su-Xi-Chang area, China. Environ. Geol. 52:27. Sibson RH (1985) A note on fault reactivation. J. Struct. Geol. 7:751–754. Simo JC (1992) Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Com- put. Meth. Appl. Mech. Eng. 99:61–112. Simo JC, Hughes TJR (1998) Computational Inelasticity Springer, New York. Sone H, Zoback MD (2011) Visco-plastic properties of shale gas reservoir rocks In 45th US Rock Mechanics Geomechanics Symposium, number ARMA 11-417, San Francisco, CA. SpaceClaim A (2018) De-featuring model for CAE. StreitJE,HillisRR(2004) Estimatingfaultstabilityandsustainablefluidpressures for underground storage of CO2 in porous rock. Energy 29:1445–1456. Szulczewski ML, MacMinn CW, Herzog HJ, Juanes R (2012) Lifetime of carbon capture and storage as a climate–change mitigation technology. P. Natl. Acad. Sci. USA 109:5185–5189. Teatini P, Castelletto N, Gambolati G (2014) 3D geomechanical modeling for CO2 geological storage in faulted formations. a case study in an offshore northern Adriatic reservoir, Italy. Int. J. Greenh. Gas Con. 22:63–76. Tinti E, Cocco M, Fukuyama E, Piatanesi A (2009) Dependence of slip weakening distance (dc) on final slip during dynamic rupture of earthquakes. Geophys. J. Int. 177. Trelis (2017) User Documentation 16.1 CSimSoft. Uenishi K, Rice JR (2003) Universal nucleation length for slip-weakening rupture instability under nonuniform fault loading. J. Geophys. Res. 108. v. Thienen-Visser K, Fokker PA (2017) The future of subsidence modelling: com- paction and subsidence due to gas depletion of the groningen gas field in the netherlands. Netherlands Journal of Geosciences 96:s105—-s116. 145 Vasco D, Harness P, Pride S, Hoversten M (2017) Estimating fluid-induced stress change from observed deformation. Geophys. J. Int. 208:1623–1642. Verdon J, Kendall JM, White D, Angus D (2011) Linking microseismic event observations with geomechanical models to minimise the risks of storing CO2 in geological formations. Earth Planet Sc. Lett. 305:143–152. Vidal-Gilbert S, Nauroy JF, Brosse E (2009) 3D geomechanical modelling for CO2 geologic storage in the Dogger carbonates of the Paris Basin. Int. J. Greenh. Gas Con. 3:288–299. Vidal-Gilbert S, Tenthorey E, Dewhurst D, Ennis-King J, Ruth PV, Hillis R (2010) Geomechanical analysis of the Naylor Field, Otway Basin, Australia: Implica- tions for CO2 injection and storage. Int. J. Greenh. Gas Con. 4:827–839. Vilarrasa V, Makhnenko R, Gheibi S (2016) Geomechanical analysis of the influ- ence of co 2 injection location on fault stability. J. Rock Mech. Geotech. Eng. 8. Wang HF (2000) Theory of Linear Poroelasticity with Applications to Geomechan- ics and Hydrogeology Princeton University Press. Wang L, Kwiatek G, Rybacki E, Bonnelye A, Bohnhoff M, Dresen G (2020) Labo- ratory study on fluid-induced fault slip behavior: The role of fluid pressurization rate. Geophys. Res. Lett. 47. Wei S, Avouac J, Hudnut K, Donnellan A, Parker J, Graves R, Helmberger D, Fielding E, Liu Z, Cappa F, Eneva M (2015) The 2012 brawley swarm triggered by injection-induced aseismic slip. Earth and Planetary Science Letters 422. White MD, Esser RP, McPherson BP, Balch RS, Liu N, Rose PE, Garcia L, Ampomah W (2017) Interpretation of tracer experiments on inverted five–spot well–patterns within the western half of the Farnsworth Unit Oil Field. Enrgy. Proced. 114:7070–7095. Won CK, Paul S (2016) Seismicity on basement faults induced by simultaneous fluid injection–extraction. Pure Appl. Geophys. 173:2621–2636. Won-Young K (2013) Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio. J. Geophys. Res. 118:3506–3518. Yaghoubi AA, Zeinali M (2009) Determination of magnitude and orientation of the in-situ stress from borehole breakout and effect of pore pressure on borehole stability-case study in Cheshmeh Khush oil field of Iran. J. Pet. Sci. Eng. 67. 146 Yang AP, Eacmen JC (2003) Fault identification and modeling for production enhancement In SPE Western Regional AAPG Pacific Section Joint Meeting, number SPE 83495, p. 13607, Long Beach, CA, USA. Yang Z, Juanes R (2018) Two sides of a fault: grain-scale analysis of pore pressure control on fault slip. Phys. Rev. E 97:022906. Yerkes RF, Castle RO (1976a) Seismicity and faulting attributable to fluid extrac- tion. Eng. Geol. 10:151–167. Yerkes RF, Castle RO (1976b) Seismicity and faulting attributable to fluid extrac- tion. Eng. Geol. 10:151–167. Zhao X, Jha B (2019) Role of well operations and multiphase geomechanics in controlling fault stability during CO 2 storage and enhanced oil recovery. J. Geophys. Res. 124. Zhao Y, Borja RI (2020) A continuum framework for coupled solid deforma- tion–fluid flow through anisotropic elastoplastic porous media . Comput. Meth. Appl. Mech. Eng. 369:113–225. 147

Abstract (if available)

Abstract Mitigating the hazard of fault reactivation during water, gas injection and oil production is important for ensuring environmental sustainability of geologic carbon sequestration and hydrocarbon recovery in the long term. In this Dissertation, we present a novel multiphase flow-geomechanics framework that can capture the nonlinearities of material (poroplasticity), kinematic (finite strain) and contact (dynamic fault slip). ❧ Using our high-resolution framework, we compare the fault destabilization effects of gas injection with those of oil production. Focusing on the confined multiphase reservoir located in the hanging wall block of a normal fault, we analyze the fault tractions, slip, and reservoir pressure to understand how a fault is destabilized in successive steps, how hypocenter and aftershock locations are selected, and what role the well type parameter (injection or production) can play in controlling the fault stability. We extract and summarize the characteristic features of induced fault activation in terms of the dependence of Coulomb Failure Stress (CFS) and the Fault Stress Ratio (FSR) on the dimensionless reservoir pressure perturbation caused by injection or production. ❧ By combining the finite strain poroplasticity with Mohr-Coulomb theory of fault failure, we conduct a fully-coupled multiphase flow and geomechanical study of production-induced fault slip in a poroplastic reservoir. We evaluate the impact of plasticity on reservoir pressure, deformation, induced stress, and fault slip by comparing infinitesimal strain elastic and finite strain poro-elastoplastic models. For real-world applications, we consider different scenarios where the reservoir is either mechanically weaker or stronger than the caprock. ❧ We also apply our workflow to incorporate complex geological structure, petrophysical heterogeneity, and multi-well schedule of a real oilfield. Our multiphase geomechanics formulation captures development of a field across different recovery mechanisms, from primary depletion to waterflooding to CO2 flooding, and predicts the geomechanical effects of fluid transport and storage on fault stability. We show that our workflow is able to characterize geomechanical signatures of flow compartmentalization and bypassed oil due to faults, inadequate well pattern configurations, and well rate variations.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

Reactivation of multiple faults in oilfields with injection and production

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Modeling and simulation of complex recovery processes

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

Dynamics of CO₂ leakage through faults

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

Application of data-driven modeling in basin-wide analysis of unconventional resources, including domain expertise

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

PDF

Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating

PDF

Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems

PDF

Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data

PDF

Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

PDF

Control of displacement fronts in porous media by flow rate partitioning

Asset Metadata

Creator Zhao, Xiaoxi (author)

Core Title A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 06/13/2021

Defense Date 12/13/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag finite strain,geomechanics,induced seismicity,multiphase flow,OAI-PMH Harvest,poroplasticity

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jha, Birendra (committee chair), Ershaghi, Iraj (committee member), Sahimi, Muhammad (committee member)

Creator Email xiaoxiz@usc.edu,zhaoxx91@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-414161

Unique identifier UC11668105

Identifier etd-ZhaoXiaoxi-9219.pdf (filename),usctheses-c89-414161 (legacy record id)

Legacy Identifier etd-ZhaoXiaoxi-9219.pdf

Dmrecord 414161

Document Type Dissertation

Rights Zhao, Xiaoxi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA