University of Southern California Dissertations and Theses

Reactivation of multiple faults in oilfields with injection and production

(USC Thesis Other)

Reactivation of multiple faults in oilfields with injection and production

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content REACTIVATION OF MULTIPLE FAULTS IN OILFIELDS WITH INJECTION AND PRODUCTION by Abdulrahman Abdulaziz Bubshait 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) May 2022 Copyright 2022 Abdulrahman Abdulaziz Bubshait Dedication To my parents, Abdulaziz and Noura, for their endless love, devotion and for nurturing my passion. To my best friend and wife, Reem, for her infinite love, unconditional believe, and encouragement. To my daughter, Noura, and son, Abdulaziz, for bringing the greatest joy to my life. To my brother, Nawaf, and sisters, Maitha, Ahad, and Shahad, for their boundless love and for all the joyful days we had together. To my parents-in-law, Mohammed and Amal, for their sincere support. ii Acknowledgements First and foremost, I would like to praise and thank Allah, the most gracious, who has granted me countless blessings, knowledge, and opportunities so that I have been finally able to accomplish the thesis. Apart from my efforts, the success of this thesis depends essentially on the motivation and advice provided by many others. Therefore, I would like to take this opportunity to express my gratitude to the people who have been influential in the successful completion of this thesis. I would like to show my greatest appreciation to Professor Birendra Jha, my advisor, I can’t thank him enough for his tremendous support and help. I feel motivated and encouraged every time attending his meeting. Without his encouragement and guidance, this thesis would not have materialized. My deepest gratitude and acknowledgment to Professor Iraj Ershaghi, the Director of the Petroleum Engineering Program at USC, for his endless support and guidance. My enormous appreciation is addressed to Professor Yehuda Ben-Zion for his valuable inputs as a member of the Ph.D. dissertation committee. I am indebted to Saudi Aramco for sparing no efforts in sponsoring and supporting me during this journey. iii TableofContents Dedication ii Acknowledgements iii ListofTables vi ListofFigures vii Abstract xvi Chapter1: Introduction 1 Chapter2: Methodology 5 2.1 Identifying Study Candidates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Production and Injection Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Seismic Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Focal Mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5 Stress Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.6 Building Coupled Flow and Geomechanics Model . . . . . . . . . . . . . . . . . . . . . . . 10 2.7 Flow and Geomechanics Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.8 Model Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Chapter3: MathematicalModelingandGoverningEquations 15 3.1 Coupled multiphase flow and geomechanics . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Shear failure model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 Modeling Fault Failure Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3.1 Moment Tensor Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.3.2 Potency-Magnitude Scaling Method . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Chapter4: CornerPointMeshingAlgorithm 26 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.4 Result and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Chapter5: SyntheticModel 41 5.1 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.2 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 iv 5.3 Result and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.3.1 Formal Stress Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.3.2 Coupled Flow and Geomechanics Simulation . . . . . . . . . . . . . . . . . . . . . 47 5.3.2.1 Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3.2.2 Case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Chapter6: AzleInducedSeismicEvents 59 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 6.1.1 Role of wastewater injection and water production . . . . . . . . . . . . . . . . . . 60 6.1.2 Interaction and propagation of induced stresses . . . . . . . . . . . . . . . . . . . . 62 6.1.3 Spatiotemporal resolution in multiphysics models . . . . . . . . . . . . . . . . . . . 63 6.2 Geologic structure, fluid volumes and well configuration . . . . . . . . . . . . . . . . . . . 65 6.3 Coupled multiphase flow–geomechanics–fault failure model . . . . . . . . . . . . . . . . . 67 6.3.1 Production and injection history . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.5.1 Pressure and stress distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 6.5.2 Interaction and propagation of stresses along depth and with time . . . . . . . . . 78 6.5.3 Stress rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 6.5.4 Fault stresses and induced failure mechanisms . . . . . . . . . . . . . . . . . . . . . 84 6.5.4.1 Synthetic Fault . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.5.4.2 Antithetic Fault . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.5.4.3 Evolution of principal stresses due to interaction between induced stresses 92 Chapter7: DFW-IrvingIncucedSeismicEvents 96 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.2 Geologic structure, fluid volumes and well configuration . . . . . . . . . . . . . . . . . . . 98 7.3 Coupled multiphase flow–geomechanics–fault failure model . . . . . . . . . . . . . . . . . 99 7.3.1 Bend Arch–Fort Worth Basin model . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.3.2 DFW–Irving model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.3.3 Coupling Bend Arch-Fort Worth Basin and DFW-Irving models . . . . . . . . . . . 103 7.3.4 Model Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.3.5 Production and Injection History . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 7.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.4.1 Bend Arch-Fort Worth Basin model pressure and stress distributions . . . . . . . . 108 7.4.2 DFW-Irving model pressure and stress distributions . . . . . . . . . . . . . . . . . 110 7.4.3 Interaction and propagation of stresses along depth and with time . . . . . . . . . 113 7.4.4 Stress rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.4.5 Fault stresses and induced failure mechanisms . . . . . . . . . . . . . . . . . . . . . 115 7.4.5.1 DFW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 7.4.5.2 Irving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7.4.5.3 Potency-Magnitude Scaling and Moment Tensor comparison . . . . . . 118 Chapter8: Conclusions 122 8.1 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Bibliography 125 Appendix 131 v ListofTables 3.1 Potency-magnitude scaling correlation parameters . . . . . . . . . . . . . . . . . . . . . . . 25 4.1 Meshing Model Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Meshing Model Poroelastic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.1 Synthetic Seismological Model Fault Parameters . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2 Synthetic Model Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.3 Synthetic Model Matrix Poroelastic Properties . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.4 Synthetic Model Fault Poroelastic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.1 Azle Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 7.1 DFW-Irving Model parameters including the initial stress and pressure gradient values. h is the vertical thickness of a layer or a fault.E is Young’s modulus.ν is Poisson’s ratio.b w is Biot coefficient. ψ is internal friction angle. c is cohesion. . . . . . . . . . . . . . . . . . 104 vi ListofFigures 1.1 Why stress changes over space and time? and how to model stress evolution in reservoir? are the two questions that motivated this work . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1 Stress mapping workflow that combines stress inversion with coupled flow-geomechanical simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Focal mechanism for different types of faulting regimes (a) strike slip, (b) normal and (c) reverse. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 (a) Map showing one of Azle 2015 seismic events epicenter (star) and the locations of seismic station (triangles). (b) shows the calculated focal mechanism for the selected event. 9 2.4 Stress inversion results plotted using stereonets. Directions of principal stresses can be read from the locations of the red, green, and blue points with respect to the center which represents the vertical direction, and the perimeter which represents the horizontal direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.1 A schematic (not to scale) showing a horst and graben region. . . . . . . . . . . . . . . . . 28 4.2 A schematic of the physical model showing the reservoir in yellow, synthetic fault in red and antithetic fault in green and the overall geomechanical domain . . . . . . . . . . . . . 29 4.3 Cross sectional views of the domain in the X − Z and Y − Z planes showing the proportional dimensions of the reservoir and the faults compared to the geomechanical domain. It also shows the curvature of the faults and the semi trapezoidal geometry of the reservoir. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Cross sectional views of the domain in theX-Z andY -Z planes showing the mechanical boundary conditions. The reservoir is bounded by no-flow boundaries from all six sides. . 30 4.5 Dual Geomechanics and Flow Coordinates Meshing Routine Workflow. . . . . . . . . . . . 32 4.6 Model Coordinate system and the point of origin. This cube shows all six sides, Front (F), Back (B), Right (R), Left (L), Up (U) and Down (D). . . . . . . . . . . . . . . . . . . . . . . . 33 4.7 Corner point ordering algorithm (STARS User Guide, 2016). . . . . . . . . . . . . . . . . . 34 vii 4.8 Corner point representation resulting from our proposed meshing routine. The reservoir is represented by the flow mesh and the non-reservoir domain (overburden, underburden, and sideburden) is represented by the geomechanical mesh surrounding the reservoir. Coupling between flow and geomechanics is honored within the reservoir mesh. . . . . . . 35 4.9 Curved fault representations by different meshing techniques showing the accuracy in describing curved faults and the principal stress acting on the surfaces of individual fault mesh elements. (a) Classical structured meshing technique. (b) Pillar meshing technique. (c) Corner-point meshing technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.10 Poisson’s ratio distribution in (a) the Pillar Mesh and (b) the Corner Point Mesh. Young’s Modulus in (c) the Pillar Mesh and (d) the Corner Point Mesh. . . . . . . . . . . . . . . . . 36 4.11 Permeability (md) distribution in (A) the Pillar mesh and (B) the Corner point mesh. . . . . 37 4.12 Pressure (psi) profile for the two meshes at different time steps. (a) Pressure profile for pillar mesh at t = 1200 days (b) at t = 10000 days. (c) Pressure profile in the flow grid only for corner point mesh at t = 2000 day and (d) Pressure at t = 10000 day . . . . . . . . . . . 38 4.13 Vertical Displacement (ft.) for the two meshes at different time steps. (a) Vertical displacement for the pillar mesh at t = 1200 days (b) at t = 10000 days. (c) Vertical displacement for the corner point mesh at t = 2000 days and (d) at t = 10000 days. . . . . . 39 4.14 Stress path in the failure space of Von Mises stress vs. the effective mean stress. The failure envelope is shown as a thick dash line. Stress paths of the fault cells for the two mesh types are shown for different time steps. Green triangles represent the stress path for the pillar mesh and the purple triangles represent the extrapolated failure point. Red circles represent the stress path for the corner point mesh and the blue circle represent the extrapolated failure point. Approximating the curved fault geometry (corner point) with a stair-stepped profile (pillar) leads to inaccurate stresses and fault stability predictions in case of the pillar mesh. To allow comparison between the two stress trajectories, the pillar mesh stress values are shifted to match the corner point stress values at t = 0. . . . . . . . 40 5.1 (a) Seismological model representing a normal fault with specified strike, rake, and dip angles. (b) The poroelastic model shows the reservoir in the footwall block (Reservoir (A)) and the hanging wall block (Reservoir (B)). The overburden and underburden layers, shown with a hashed region, have a permeability value lower than that of the reservoir. The mechanical boundary conditions are shown. All boundaries are no-flow boundaries. A well is located in either the footwall reservoir (Case 1) or the hanging wall reservoir (Case 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2 The seismological model showing three seismic events with the respective focal mechanisms on the fault plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.3 (a) Gausian distributions of strike, dip, and rake angles. (b) sample of synthetic focal mechanisms data used for stress inversion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 viii 5.4 Permeability distribution and well location in the poro-elastoplastic model. (a) Case 1 and (b) Case 2 correspond to two different simulation models with two different well locations. 46 5.5 Stress inversion results plotted using stereonets. Directions of principal stresses can be read from the locations of the red, green, and blue points with respect to the center which represents the vertical direction and the perimeter which represents the horizontal direction. The difference between the maximum and minimum principal stress directions for the 1500 events is 90 degree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.6 Pressure at t=1428.6 and t=10000 day for Case 1. . . . . . . . . . . . . . . . . . . . . . . . . 47 5.7 Effective Mean Stress at t=1428.6 and t=10000 day for Case 1. Compression is positive. . . 48 5.8 Effective Plastic Strain at t=1428.6 and t=10000 for Case 1. . . . . . . . . . . . . . . . . . . 49 5.9 Vertical Displacement at t=1428.6 and t=10000 for Case 1. Unit is feet. . . . . . . . . . . . . 50 5.10 von Mises stress at t=1428.6 and t=10000 for Case 1. . . . . . . . . . . . . . . . . . . . . . . 50 5.11 Stress path in the failure space of von Mises stress vs. effective mean stress. The failure envelope is shown as a dash line. Stress paths of four fault cells are shown by plotting the stresses at five timesteps (value of t indicates the timestep number, not actual time in days). Each time step has a unique symbol. At the first timestep (t=1), all four cells have the same von Mises stress (3000 psi) and the effective mean stress (7000 psi) because the initial stresses are uniform in the model (see Table 2). At t=2, the cell next to well perforation moves, and other cells stay at the initial point . . . . . . . . . . . . . . . . . . . 52 5.12 Depth profiles of effective plastic strain at different time steps (time in days shown in the legend) during the simulation. The nucleation point (hypocenter) is at the well perforation depth and failure propagates both updip and downdip, with a larger updip span because of the traction-free top boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.13 Total Effective Plastic Strain for Case 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.14 Pressure at t=1666 and t=10000 for Case 2. Unit is psi . . . . . . . . . . . . . . . . . . . . . 54 5.15 Effective Mean Stress at t=1666 andt=10000 day for Case 2. Unit is psi . . . . . . . . . . . 55 5.16 Effective Plastic Strain at t=1666 and t=10000 day for Case 2 . . . . . . . . . . . . . . . . . 55 5.17 TVertical Displacement at t=1666 day and t=10000 day for Case 2. . . . . . . . . . . . . . . 56 5.18 Von Mises Stress at t=1666 and t=10000 for Case 2. . . . . . . . . . . . . . . . . . . . . . . . 57 5.19 Stress path in the failure space of von Mises stress vs. effective mean stress. The failure envelope is shown as a dash line. Stress paths of four fault cells are shown by plotting the stresses at five timesteps (value of t indicates the timestep number, not actual time in days). Each time step has a unique symbol. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 ix 5.20 TDepth profiles of the effective plastic strain on the fault at different time steps. Legend shows time in days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.21 Total Effective Plastic Strain for Case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.1 (a) A map view of the Azle region (marked by red star in the inset map) at Ellenburger depth showing 2013-2014 events and the synthetic–antithetic fault pair. Re-drawn with permission from [32]. (b) Stratigraphic section of the Newark East oilfield along X-X’ line in (a), showing Barnett production, Ellenburger injection, and the fault pair. The figure is not to scale. The stars indicate the general locations of hypocenters during induced seismicity. Interaction between contraction of Barnett and expansion of Ellenburger, shown by red circles with inward and outward pointing arrows, and injection-induced pressure diffusion into basement can generate induced stresses on the faults, which will vary laterally across fault blocks and vertically along the faults. . . . . . . . . . . . . . . . 62 6.2 (a) 3D view of the domain showing reservoir (green), aquifer (blue), overburden (yellow), underburden (red), synthetic fault (red surface) and antithetic fault (yellow surface). (b) The corner-point geomechanical mesh and (c) the corner-point flow mesh generated using the proposed meshing algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 Barnett permeability map shown at the Barnett Top surface. It shows the effect of hydraulic fracturing on enhancing the permeability around each gas producer’s horizontal section in the reservoir. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.4 Observed (circles) and simulated (lines) (a) Well API-42367346930000 and (b) Well API- 42497368750000 flowing wellhead pressure calculated from the reported bottom hole pressure [32, 15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.5 Map showing the location of gas producers (highlighted in red) and water injectors (highlighted in blue). The area outside the dashed green rectangle represents a buffer zone where no gas producers are considered in order to minimize the boundary effect. Below, we call the two injectors inside the dashed green box as near-fault injectors and the other three injectors as far-field injectors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.6 Observed (circles) and simulated (lines) (a) gas monthly production rate q g , (b) water monthly production rate q wp , (c) water monthly injection rate q wi , (d) cumulative gas production volumeQ g , (e) cumulative water production volumeQ wp , and (c) cumulative water injection volumeQ wi in the Azle field plotted with time. All volumes are at surface pressure and temperature conditions. (g) Bubble map of cumulative produced and injected volumes per well. The bubble size is linearly proportional to produced or injected fluid’s volume. Produced water circles have been magnified 3.3 times for the ease of visualization. 73 6.7 Evolution of the fault failure magnitude (plastic strain) and the fault pressure in a cell of the synthetic fault for six different values of the vertical-to-horizontal permeability ratio. The cell is located in the basement (I=22, J=12, K=10). . . . . . . . . . . . . . . . . . . . . . 75 x 6.8 Barnett overpressure∆ p at the selected time steps: a, b, c, d, e, and f. Fracturing-induced permeability enhancement around gas producers lead to heterogeneity in the pressure field. Due to smaller injection volumes, the footwall block experiences larger depletion compared to the hanging wall block. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 6.9 Change in the effective mean stress ∆ I ′ 1 (with respect to the initial value) at the selected timesteps: a, b, c, d, e, and f in Barnett. ∆ I ′ 1 , which is positive in compression, shows a monotonic response to∆ p in Figure 6.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.10 Ellenburger overpressure ∆ p fields at different time steps: a, b, c, d, e, and f. The uniform nature of the aquifer permeability causes ∆ p in both fault blocks to evolve monotonically with injection. Low values of the fault horizontal permeabilityk h , relative to the Ellenburger permeability, impede pressure communication between the two fault blocks, however both blocks exhibit relatively balanced overpressure due to water injection. 79 6.11 Changes in Ellenburger’s effective mean stress ∆ I ′ 1 at different time steps: a, b, c, d, e, f. Changes are with respect to the initial value. While the changes mostly follow the behavior of Ellenburger pressure, there are areas with localized higher-order effects in ∆ I ′ 1 that do not directly follow Ellenburger’s pressure and indicates stress propagation from Barnett. These are areas with small magnitudes of ∆ I ′ 1 (closer to zero) where injection-induced tensile stresses are counterbalanced by Barnett’s production-induced compressive stresses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.12 (a) Overpressure ∆ p in the flow domain along a XZ vertical section. (b) Vertical displacement u z in the entire geomechanical domain along the same vertical section. (c) Depth profiles of ∆ p, changes in the effective mean stress ∆ I ′ 1 and von Mises stress ∆ √ J 2 , andu z along a column of cells in the hanging-wall block. (d), (e), and (f) are similar plots in the footwall block. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 6.13 Change in the maximum and minimum effective principal stress angles ( ∆ θ 1 and∆ θ 3 ) from their initial values in Barnett and Ellenburger. (a) shows the Barnett overpressure in Dec 2020. (b) shows∆ θ 1 and∆ θ 3 along section A-A’ of the synthetic fault at Barnett depth at several time steps. (c) shows∆ θ 1 and∆ θ 3 along an east-west section B-B’. (d) shows the Ellenburger overpressure in Dec 2020. (e) shows∆ θ 1 and∆ θ 3 along section C-C’ of the synthetic fault at Ellenburger depth at several time steps. (f) shows∆ θ 1 and ∆ θ 3 along an east-west section D-D’. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 6.14 Rotation of the minimum effective principal stress σ ′ 3 and the intermediate effective principal stressσ ′ 2 on December 2020 in (a) Barnett and (b) Ellenburger. Well rates have been projected to future based on June 2019 rates. Since hydraulic fractures tend to orient perpendicular toσ ′ 3 , rotation in principal stresses have implications for future drilling and fracturing activities in the field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 xi 6.15 (a) Failure magnitudeε p plotted on the synthetic and antithetic faults and overpressure∆ p plotted as contours on a horizontal cross-section of Ellenburger. Bottom-hole locations of the water injectors are shown with triangles. (b) shows a rotated view of (a) without the Ellenburger section. Locations and magnitudes of the failure events on the faults show that deeper section of the synthetic fault experiences larger magnitude events compared to events located at or above the injection depth. . . . . . . . . . . . . . . . . . . . . . . . 85 6.16 Evolution of overpressure ∆ p, change in von Mises stress ∆ √ J 2 , change in effective mean stress∆ I ′ 1 , and instantaneous failure straindε p in synthetic fault cells located at the depths of Ellenburger layers. K=5 through 10 are top six layers in Ellenburger.∆ p,∆ √ J 2 , and∆ I ′ 1 are plotted against the left axis whiledε p is plotted against the right axis. Axis labels, shown only in the bottom left plot, are identical for all other plots. . . . . . . . . . . 87 6.17 Stress trajectories of synthetic fault cells in theI ′ 1 − √ J 2 stress invariant space over the 2000-2020 period. K=5 through 10 are the top six layers in Ellenburger. The K=5 cell is chosen as the reference cell from which other cells are compared. Hence, all cells start at the same initial point, which corresponds to the initial stress state of K=5. The start time of gas production from Barnett (t GP ), the start time of water injection in Ellenburger (t WI ), and the failure onset time (t p i ) are marked. Failure in K=5, which is in Ellenburger, begins before injection begins in Ellenburger. This suggests propagation of production-induced stresses from Barnett into Ellenburger. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.18 Evolution of overpressure ∆ p, change in von Mises stress ∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and failure strain ε p in synthetic fault cells located at the depths of Ellenburger layers for the two simulation runs with different initial stress values (Run-2) and orientations (Run-3). (e) compares the evolution ofε p in synthetic fault cell (22,15,5) for Run-1 (Base Case), Run-2, and Run-3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.19 Top view of model showing the change in effective mean stress distribution at last time step for Barnett and Ellenburger at two different additional simulation runs with initial stress values and orientations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 6.20 Top view of model showing the vertical displacement distribution at last time step at two different additional simulation runs with initial stress values and orientations. . . . . . . . 90 6.21 Spatiotemporal evolution of overpressure ∆ p, change in the effective mean stress ∆ I ′ 1 and change in the von Mises stress∆ √ J 2 , and instantaneous failure straindε p for four cells in the antithetic fault within Ellenburger interval. Top row: locations of the cells in 3D. Middle and bottom rows: ∆ p, ∆ √ J 2 , ∆ I ′ 1 are plotted against the left axis while dε p is plotted against the right axis. Middle row, which corresponds to layer K=5, shows failure events on the dε p curve after water injection begins in the near-fault injectors. The injection start time,t WI , is marked. Bottom row shows that the K=6 layer does not experience failure within the simulated period. . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.22 Stress paths of antithetic fault cells in theI ′ 1 − √ J 2 space. Two cells from Ellenburger’s top layer (K=5) and second layer (K=6) are plotted. The top layer cell shows a more complex stress path (non-monotonicity inI ′ 1 ) due to its proximity to Barnett and hence to production-induced stresses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 xii 6.23 Temporal evolution of the change in effective principal stresses, ∆ σ ′ 1 , ∆ σ ′ 2 and ∆ σ ′ 3 , in selected cells in (a) the antithetic fault and (b) the synthetic fault. Cells that are vertically closer to Barnett exhibits stress evolution after the start of production, while cells that are located deeper do not exhibit the same characteristics. As injection starts near the faults, stress acting on all four cells decrease due to the increase in pressure. The separation between principal stresses increase due to interaction between production- and injection-induced stresses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.24 3D views of the synthetic and antithetic faults showing the induced failure events at the same four timesteps analyzed above: a, b, c, and d. All events occurring at or before a timestep are shown at that timestep. The size of the circle or square is a proxy for the magnitude of the calculated event. The events are more concentrated and are larger in magnitude in the basement interval. Green shaded cells represent fault intersection with Barnett, and blue shaded cells represent fault intersection with Ellenburger. (e) Temporal distribution of simulated and observed events on synthetic and antithetic faults. The Nov 2013–Jan 2014 seismic cluster and a recent Mb 2.8 event in June 2017 are shown with asterisks. (f) Depth and magnitudes of the observed seismic events. . . . . . . . . . . . . . 95 7.1 (a) A map view of the DFW-Irving study area in North East Taxes showing the 2008-2015 seismic events and the major fault lines. County lines and county names are shown in grey color. Two major faults bound the study area from the north and the east, namely Muenster Arch fault and Ouachita Thrust Front. (b) Vertical stratigraphic section (schematic not to scale) of the area along a line AB oriented west-north-west to east-south-east. It shows gas and water production from hydraulically fractured Barnett wells and water injection into Ellenburger wells. Eight out of a total of nine faults are shown. Fault-2 does not become visible in this cross-section. The stars indicate the general hypocenters locations. Volumetric contraction of Barnett is shown by inward pointing arrows on a red circle. Volumetric expansion of Ellenburger is shown by outward pointing arrows. Fluid and pressure diffusion into faults from Ellenburger to basement is shown by wiggly arrows. . 100 7.2 (a) 3D view of the Bend Arch-Fort Worth basin flow and geomechanics domains. The domain is bounded by Muenster arch in the northeast and by Ouachita thrust front in the southeast. (b) 3D view of the DFW-Irving flow and geomechanics domains showing Barnett formation (green) Ellenburger formation (blue) and faults ((red) for portions that are represented by flow and geomechanis mesh and (purple) for portions that are represented by geomechanis mesh only). Triangles (red) represent the surface location of the gas and brine producers and inverted triangles (blue) represent the surface location of water injection wells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.3 Coupling of the two models, (a) Bend Arch-Fort Worth basin and (b) DFW-Irving model. Injection-induced pressure front in the basin model propagates to western and southern boundaries of the Ellenburger formation in the DFW-Irving model. Coupling ensures that the pressure boundary condition for the faults in the DFW-Irving model evolves in sync with the long-term injection activity in the basin.. . . . . . . . . . . . . . . . . . . . . . . . 106 xiii 7.4 (a) Histogram of Barnett porosity, (b) histogram of Barnett log-permeabilities logk x in I direction and logk y in J direction, (c) histogram of Ellenburger porosity, (d) histogram of Ellenburger log-permeability in I and J directions. Vertical-to-horizontal permeability ratio k z /k x is 0.1. (e) Top view of the Barnett permeability map showing the effect of hydraulic fracturing on enhancing the permeability field around each horizontal gas producer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.5 Observed (circles) and simulated (lines) (a) gas monthly production rate q g , (b) water monthly injection rateq wi , (c) cumulative Basin and DFW water injection volumeQ wi in addition to the south and west water influx in the DFW-Irving field plotted with time. All volumes are at surface pressure and temperature conditions. . . . . . . . . . . . . . . . . . 108 7.6 (a) and (b) Basin overpressure ∆ p at the selected time steps. (c) and (d) change in the effective mean stress ∆ I ′ 1 (with respect to the initial value). Ouachita thrust front bounding the basin from southeast and low horizontal permeability across DFW airport fault helps in channeling pressure front towards Irving area. . . . . . . . . . . . . . . . . . 109 7.7 Barnett (top) and Ellenburger (Bottom) overpressure ∆ p at selected time steps: a, b, c, and d. While fracturing-induced permeability enhancement around gas producers lead to heterogeneity in the pressure field in Barnett, ∆ p evolves more uniformly due to Ellenburger uniform permeability. The effect of Basin pressure front is clearly visible as ∆ p evolves more rabidly due to the generally low permeability of faults inx direction and specifically fault 5 and fault 8 channeling pressure front toward fault 6. . . . . . . . . . . . 111 7.8 Barnett (top) and Ellenburger (Bottom) change in the effective mean stress ∆ I ′ 1 (with respect to the initial value) at the selected time steps: a, b, c, and d. ∆ I ′ 1 , which is positive in compression and negative in tensile, shows a generally monotonic response to∆ p in Figure 7.7. Some areas, such as area highlighted with red dashed line, do not directly follow Ellenburger’s pressure and have a smaller magnitudes of∆ I ′ 1 (closer to zero) which indicates stress propagation from Barnett. This shows that injection-induced tensile stresses are counterbalanced by Barnett’s production-induced compressive stresses. . . . . 112 7.9 (a) Overpressure∆ p in Barnett and Ellenburger at the end of simulation time. (b), (c), and (d) are depth profiles of ∆ p, changes in the effective mean stress ∆ I ′ 1 and von Mises stress ∆ √ J 2 , andu z along a column of cells in the selected location shown in (a). . . . . . . . . . 114 7.10 Rotation of the minimum effective principal stress σ ′ 3 and the intermediate effective principal stressσ ′ 2 on December 2020 in (a) Barnett and (b) Ellenburger. . . . . . . . . . . 115 xiv 7.11 (a) Ellenburger∆ p at December 2020 highlighting the location of fault-4 (red) (b) 3D view of fault-4 showing failure magnitudeε p (c) Evolution of overpressure∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure strain dε p in fault cells located at the depths of Ellenburger layers. K=16 through 20 are top five layers in Ellenburger. ∆ p, ∆ √ J 2 , and∆ I ′ 1 are plotted against the left axis whiledε p is plotted against the right axis. Axis labels, shown only in the bottom left plot, are identical for all other plots. (d) Stress trajectories of fault cells in the I ′ 1 − √ J 2 stress invariant space over the 2005-2020 period. The K=16 cell is chosen as the reference cell from which other cells are compared. Hence, all cells start at the same initial point, which corresponds to the initial stress state of K=16. The start time of gas production from Barnett (t GP ), the start time of water injection in Ellenburger (t WI1 ) and (t WI2 ) corresponding to the two water injectors, the failure onset time (t p i ), and the injector shut-in (t WI SI ) are marked. . . 117 7.12 (a) Ellenburger ∆ p at December 2020. (b) Fault 4 failure magnitude ε p . (c) Evolution of overpressure∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure straindε p in different fault cells located at overburden, Barnett, Ellenburger and basement layers. K=15 is located in overburden, K=19 is located in Barnett, K=26 and 28 are located in Ellenburger, and K=36 is located in basement.(d) Stress trajectories of fault cells in theI ′ 1 − √ J 2 stress invariant space over the 2005-2020 period. The failure onset time (t p i ) is marked. (e) fault depth profile of how change in vertical stress inX direction evolves in timeδ x σ v . . . . . . . . . . . . . . . . . . . . . . . 119 7.13 (a) Temporal distribution of simulated and observed failure events showing fault-4 observed failure events (black asterisks), fault-6 observed events (red asterisks), fault-4 simulated failure events (blue) fault-6 simulated failure events (yellow), simulated failure events calculated using Moment Tensor are marked by triangles, and simulated failure events calculated using Potency-Magnitude Scaling are marked by circles. (b) and (c) log of number of eventsN greater than or equal to the magnitudeM versus failure magnitude M for observed (asterisks) and simulated (circles) failure events occurred on fault-4 and fault-6 respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 xv Abstract Propagation of production- and injection-induced stresses from reservoir layers to basement faults is not well understood, especially in reservoirs with hydraulically segregated production and injection units sit- uated across low-permeability structural or stratigraphic boundaries. Having a comprehensive under- standing of how flow-induced stresses propagate and induce fault reactivation has a significant impact in optimizing oil and gas extraction, wastewater injection, and CO2 sequestration strategies, which can promote environmentally safe subsurface operations. Numerical modeling and simulation is the only tool available to forecast stress evolution at field scales and quantify the risk of fault reactivation under different geologic and fluid flow scenarios. Numerical modeling requires the solution of the coupled equations of fluid flow, rock deformation, and fault failure over a computational mesh discretizing the geologic domain. This is challenging in fields with different injection and production intervals, hundreds of hydraulically fractured wells with time-varying flow schedules, and multiple curved and connected faults. We present a framework to model spatiotemporal evolution of production- and injection-induced stress based on a coupled formulation of multiphase flow, poroelasticity, and fault failure over a dual mesh where the flow mesh discretizes the flow region at higher resolution than the resolution of the geomechanics mesh dis- cretizing both flow and no-flow (overburden and basement) regions. We develop a corner point geometry meshing algorithm that can model reservoirs with multiple curved and connected faults. This expands the applicability of many industry-standard simulation software limited to structured meshes and reduces the need for building unstructured meshes at high computational costs. We compare two methods to model xvi seismic events from fault failure events in our simulations, namely, moment tensor and potency-magnitude scaling methods and show that the potency-magnitude scaling method provides a promising alternative to the moment tensor method for event magnitude calculation from quasi-static simulations. We use the proposed framework to study 2008-2015 seismicity around Azle, Dallas, Fort Worth and Irving cities in Texas, USA by modeling the effects of gas and water production from Barnett Shale and wastewater injec- tion into Ellenburger on the stress states of the shale reservoir, the aquifer, and the granitic basement. We analyze stress propagation onto faults to discover new mechanisms of faults reactivation in such complex fields. We also show that flow-induced anisotropy and heterogeneity in the stress state can be analyzed in ultra-low permeability reservoirs undergoing hydraulic fracturing to optimize the trajectories of new and infill wells and their hydraulic fracturing jobs. xvii Chapter1 Introduction With the increasing awareness about stress sensitivity of reservoirs driven by induced seismicity and ground deformation events occurring around oilfields, a comprehensive understanding of induced stress evolution has became essential. Understanding reservoir stresses plays a major role in optimizing reser- voirs depletion strategies and drilling and stimulation programs as they determine rock deformation, con- trol changes in hydraulic and elastic properties, and affect pore pressure. The state of stress in a reservoir is a dynamic quantity that results from different independent pa- rameters such as the regional tectonic stresses and the reservoir pore pressure. Thus, a reservoir that experiences a change in pore pressure due to injection and production operations or experiences earth- quakes will exhibit a change in the principal stress directions and magnitudes. Analysis of borehole images and crossed-dipole sonic logs can provide us with an understanding of the stress field around wellbores. But it is difficult to extrapolate these stresses away from the wellbore into the reservoir because stress is a state variable governed by coupled partial differential equations of momentum and mass balance, not a material property like permeability which can be distributed statistically as random fields. This limits the use of statistical approaches in mapping reservoir stresses especially the change in stress induced by fluid flow. 1 Realizing the importance of having a comprehensive understanding of stress evolution we intended to answer three crucial questions: why stress changes over space and time, how does stress evolution affects faults that are deeper than the reservoir, and what mechanisms are behind it. Stress is a dynamic quantity resulting from different independent parameters, such as regional tectonic stresses and reservoir pore pressure. Thus, injection and production operations can change stress directions and magnitudes in the reservoir. The second question is how can we determine stress evolution in reservoirs? borehole images can provide us with an understanding of the stress field around wellbores. But it is difficult to extrapolate As stress is a dynamic state variable, not a material property like permeability Thus limits the use of statistical approaches, figure 1.1 The objective of this research is to model spatiotemporal evolution in reservoir stress magnitude and directions caused by injection and production activities in order to optimize depletion strategies and to understand the mechanisms behind fault reactivation. We develop such a framework by combining stress inversion and coupled flow and geomechanical simulation. We use a fully coupled poro-elastoplastic sim- ulation approach in CMG STARS to examine production-induced stress changes in a reservoir and plastic failure within a fault. Using the analogy between Mohr-Coulomb plastic failure and Mohr-Coulomb fric- tional failure processes, we compute slip on the fault in terms of plastic strain. Geometric characterisitics of slip is used to derive the focal mechanisms that we use in MSATSI, a MATLAB Package for Stress In- version, to obtain the magnitudes and directions of three principal stresses. The result is a framework that provides a 4D (time-lapse) stress map of a reservoir while honoring both the production/injection and microseismic data. Chapter 2 presents our methodology to develop a sequential workflow to solve the three sub-problems of flow (single-phase and multiphase), mechanical equilibrium, and fracture propagation sequentially. We discuss different types of sequential coupling methods and their advantages and disadvantages. 2 Why stress changes over space and time? How to model stress evolution in reservoirs? Near Wellbore Reservoir Regional Figure 1.1: Why stress changes over space and time? and how to model stress evolution in reservoir? are the two questions that motivated this work 3 Chapter 3 discusses the mathematical modeling and governing equations used in the coupled flow and geomechanics simulator, in addition, discusses the two methods that we utilized to model seismicity. Chapter 4 presents the novel corner point meshing algorithm that we developed to model reservoirs with multiple, curved, and intersecting faults. The developed algorithm provides the 3D coordinates which allow for the representation of faults as 3D objects and to assign physical properties to them. Developing this meshing algorithm was crucial to resolve and model real complex geological setups as traditional structured grids cannot represent curved faults, connected faults, or fault intersections without sacrificing either their geometry or their geomechanical properties. Chapter 5 presents the synthetic model that we utilized to initially examine and validate our workflow upon due to the complexity of real faults and tectonic stress regimes and the errors involved in seismic data interpretation. Chapter 6 showcases the application of our workflow meshing algorithm in a real field. we studied the induced seismic events in Azle that occurred in 2014. the workflow allowed us to examine the causality behind the activation of the two, synthetic and antithetic, faults. Chapter 7 presents how we expanded the utilization of our workflow and our meshing logarithm to develop a dual model with several curved and connected faults to study two seismic clusters that occurred in the DFW-Irving region in 2008-2015. 4 Chapter2 Methodology In order to address our research problem, we proposed a framework, Figure 2.1, that couples fault reacti- vation modeling with coupled flow and geomechanics simulation. The induced seismic events resulting from fault reactivation due to production and injection, allow us to integrate fault reactivation model and flow and geomechanics simulation into a feedback framework. In this chapter, we describe the framework we used to couple fault reactivation modeling with coupled flow and geomechanics simulation and the set of tools and software that we utilized in the process. 2.1 IdentifyingStudyCandidates The framework starts by Identifying and selecting real field candidates to study and to utilize to validate our proposed framework. To select a real field for the application of the framework, three crucial conditions are required to be present. The first condition is the existence of production and injection activities in the vicinity of the field and the existence of production and injection data, this is crucial because our framework is targeting production and injection-induced stresses. the second condition is the existence of faults in the domain of the study, this is important because faults will host the induced seismic events. the third condition is the historical occurrence of observed seismic events. Analyzing observed seismic events will enable us to build and initialize the model. 5 Focal Mechanism Fault Stresses Magnitude and Directions Seismic Event Build/Update Model Reservoir Stresses Magnitude and Directions Production & Injection HASH 1.2 MSATSI FRAMEWORK Meshing Algorithm - Designed for modeling curved & connected faults in high resolution Fault Reactivation Modeling - PROD/INJ Match? - Seismic Events Match? Yes No Calibrate Strategy Optimization Stress Inversion Flow & Geomechanic Simulation Published Literature Reservoir Stress Magnitude & Orientation Modeling Figure 2.1: Stress mapping workflow that combines stress inversion with coupled flow-geomechanical simulation. In our work we started with applying the workflow on a synthetic model, this was required due to the complexities of the geological settings of real fields. This step was critical in paving the road for the application of the framework on real field data as it eliminated the complexities and uncertainties associated with real data. In addition, it provided us with the experience required to utilize the different applications required to achieve our objective. Following that, we introduced the first real field data to our framework, which was the 2014 seismic events that occurred in Azle in the northeast of Texas. This seismic cluster occurred on a pair of synthetic and antithetic curved and connected faults. after that, we elected to expand the application of our workflow to a more geologically complex field. The two seismic clusters that occurred in 2008-2015 and 2015 in DFW-Irving in the northeast of Texas presented a challenge in understanding the triggering mechanism behind these seismic clusters. The two clusters occurred in a heavily faulted domain with nine mapped curved faults. The real fields that we selected to study fulfill all three conditions that are required for the application of the framework. 6 2.2 ProductionandInjectionData Reliable production and injection rates and volumes data are necessary to build a faithful model because stress evolution is derived by pressure changes resulting from fluid extraction and injection. After we identified the field we are interested in studying, we surveyed the literature for detailed production and injection rates and volumes data, while much of the literature was helpful in providing good data to build our model, it was challenging to obtain detailed well by well production and injection data. We obtained the well data including monthly production and injection volumes and well locations from Texas Railroad commission and IHS Markit database. The Production and Injection data consist of well locations, well by well monthly and cumulative production and injection gas and water rates, injectors bottom hole pressures. 2.3 SeismicData Seismic data for the observed events that occurred in the field we are studying were obtained from IRIS database. The seismic data consists of seismic events clusters location, time, the magnitude in addition to the waveform data for the seismic events that occurred in the field for a specific time frame. 2.4 FocalMechanism Deformation associated with seismic events that occur on faults surfaces can be described using focal mechanisms. It describes the dip and strike angles of the slip in addition to the rake, which describes the orientation of the slip on the fault surface. The importance of the focal mechanism is that the direction of slip on a fault plane is governed by the stress field orientation at the time of rupture. Hence, knowing the direction of the slip on a fault plane will provide some insights into the state of the stress at the time of failure that will be used in building and setting up the right model. 7 (a) (b) (c) Strike slip Normal Reverse Figure 2.2: Focal mechanism for different types of faulting regimes (a) strike slip, (b) normal and (c) reverse. The focal mechanism is represented by beach ball representations which are plotted using a lower- hemisphere stereographic projection. The dark part represents the compression, and the white represents tensile. Figure 2.2 shows different representations of focal mechanisms for different faults types. The first step in solving for the focal mechanism is obtaining the seismic waveforms event. The waves generated in the hypocenter by the seismic event propagate in all directions and get received and recorded by the many distributed seismic stations. the seismic station records displacement in all three axes;X,Y , andZ which represent two orthogonal horizontal directions, and one vertical. depending on the location of the seismic station, the first arrived seismic wave polarity can wiggle up or down. The direction of waveform first arrival polarity is different from a station to a station and is dependent on the side of the fault if it experiences compression or tensile. analysis of the first arrival polarity direction for received waveforms can be used to back-calculate the focal mechanism. We used HASH, a FORTRAN program to calculate the focal mechanism based on the polarity of the first arrival waveform [30]. To calculate the focal mechanism, we collected data related to the seismic station location and polarity of each waveform received at such stations. 8 (a) (b) Figure 2.3: (a) Map showing one of Azle 2015 seismic events epicenter (star) and the locations of seismic station (triangles). (b) shows the calculated focal mechanism for the selected event. Figure 2.3 shows the location of a seismic event epicenter and the location of the surrounding seismic stations, in addition, it shows the calculated focal mechanism for the identified seismic event. 2.5 StressInversion The third component is the formal stress inversion, which has proven to be a well-founded methodology to analyze stresses, fault mechanics, and earthquakes. Formal stress inversion utilizes the focal mechanism of a seismic event, the fault plane orientation, and the slip direction on the fault (Figure 5.1 (a)), to calculate the principal stress at the hypocenter of the seismic event. We used MSATSI to handle the formal stress inversion, which is a MATLAB package that enables us to invert the focal mechanisms in a discretized spatial and/or temporal grid and provide the orientation of the three principal stress axes and the relative stress magnitude [45, 46, 44]. The input parameters for MSATSI are the hypocenter coordinates, strike, dip, and rake of each focal mechanism. The results of the formal stress inversion on the synthetic model are presented in Figure 5.5, which shows the orientations of the three principal stresses locally around each of the three hypocenters. 9 (a) (c) (b) Figure 2.4: Stress inversion results plotted using stereonets. Directions of principal stresses can be read from the locations of the red, green, and blue points with respect to the center which represents the vertical direction, and the perimeter which represents the horizontal direction. In addition to stress inversion and solving for focal mechanism at the hypocenter of the seismic events, we surveyed the literature for stress analysis and faults poroelastic properties that would help us to un- derstand the stress magnitude and directions in the reservoir. 2.6 BuildingCoupledFlowandGeomechanicsModel In this section we will describe the steps we followed to build the coupled flow and geomechanics model and to prepare the input file for the simulator. We used CMG STARS coupled flow and geomechanics simulator. We started by selecting the domain that will be modeled. contrary to typical reservoir simulations, the domain in our study should be big enough to cover rock geomechanics evolution. the domain should cover the surface to the basement as faults crossed from surface to basement. Also, the domain should be big enough to minimize the effect of the boundary. It is essential that the model captures the geometry of faults and be able to translate the stresses from the reservoir to the basement. However, currently, no structured meshing techniques can faithfully capture fault geometries. For that, we developed a novel structured corner point meshing algorithm in MATLAB with the objective of modeling complex curved and connected faults. details of the corner point meshing 10 algorithm are presented in Chapter 4. Using this algorithm, we were able to construct a geomechanics mesh and a flow mesh that will be used in simulation. We started by selecting the domain that will be modeled. contrary to typical reservoir simulations, the domain in our study should be big enough to cover rock geomechanics evolution. the domain should cover the surface to the basement as faults crossed from surface to basement. Also, the domain should be big enough to minimize the effect of the boundary. It is essential that the model captures the geometry of faults and be able to translate the stresses from the reservoir to the basement. However, currently, no structured meshing techniques can faithfully capture fault geometries. For that, we developed a novel structured corner point meshing algorithm in MATLAB with the objective of modeling complex curved and connected faults. details of the corner point meshing algorithm are presented in Chapter 4. Using this algorithm, we were able to construct a geomechanics mesh and a flow mech that will be used in simulation. We initialized the model using pressure and stresses gradients based on the stress inversion and the literature for the same or an adjacent field. The pressure gradient dp i /dz was 11.28 kPa/m, the minimum horizontal stress gradientdσ 3,i /dz was 4.5 kPa/m, the maximum horizontal stress gradientdσ 2,i /dz was 4.95 kPa/m, and the vertical stress gradient wasdσ 2,i /dz was 4.95 kPa/m. The stress gradient distribution produces a normal faulting regime. Boundary conditions of the domain are as follows. The normal dis- placement in the geomechanics model is prescribed to be zero on the left, right, front, and back boundaries. The bottom boundary is fixed in all three directions. The top boundary is a traction-free boundary, and it is free to deform as stresses change within the domain. No-flow boundary conditions are prescribed on all six sides of the flow model. We used SGeMS, software for solving problems involving spatially related variables, to model the permeability and porosity fields using a multi-Gaussian random field. In our research, we studied fields with low porosity and low permeability with an average porosityϕ of 0.055-0.06 and average horizontal 11 permeabilityk Z of9.87× 10 − 21 . to allow faults pressure communication with the basement, we assigned different vertical permeability in faults k Z equal to 9.87× 10 − 14 . The initial water saturation S wi for formations in the models were estimated through several iterations of model calibration (see the Model calibration section below) where we fine-tuned the S wi field to honor the observed gas and water rates and volumes. To achieve higher degrees of faithfulness to the real data, each type of formation or geological feature was assigned specific mechanical properties. the mechanical properties we assigned for each formation or geological feature properties are Young’s modulusE, Poisson’s ratioν , Biot coefficient α . internal friction angleψ , and cohesionc. Mechanical properties of the model are based on the published literature. Faults mechanical properties are less certain than those of formations; hence, we used model calibration, where the model predicted pressures and fault failure events are compared against the measured well pressures and recorded seismicity, to estimate these parameters. Due to the ununiform nature of the grid, as each grid cell might have different dimensions than the adjacent cell, we developed a wall placement algorithm that converts the longitude and latitude of cell location toI andJ indices. In order to model the permeability enhancement due to hydraulic fracturing in the gas wells, we altered the permeability field around each gas producer’s horizontal open-hole section. The permeability increases by a factor of10 7 in the well completion block and decreases logarithmically with distance 2.7 FlowandGeomechanicsSimulation We use the CMG-STARS simulator [16], a fully coupled flow and geomechanic simulator quasi-static, to simulate the effect of production and injection on pressure and stresses, see Chapter 3. We started by building the input file which contains the the flow mesh consisting of the corner-point coordinate in a specific required CMG format. our meshing algorithm provides the mesh in this format. The input file 12 should also include a fluid parameter section where we specify the properties of the fluids in the system. Then we enter the relative permeability curves, initial pressure and permeability distribution. We use SGeMS to provide a Gaussian distribution permeability field based on a specific average, standard deviation and correlation length. The geomechanics mesh is provided in the same fashion as the flow mesh. due to scale and complexity of our model, we utilized our meshing algorithm to generate input geomechanics parameters required for the simulation in CMG format. 2.8 ModelCalibration We conduct a history matching analysis with the goal of reproducing observed production and injection volumes, in addition to seismic events magnitudes, location, and onset time undergoing the same condi- tions. This step is crucial to assure the quality of our model In order to ensure model fidelity and confidence in simulation results, we performed a multiscale cal- ibration study. We performed calibration at the well scale, reservoir layer scale, and the entire field scale, which in turn reduced uncertainty in the near wellbore properties, layer-wise properties and field-scale poroelastic properties while honoring the production and injection volumes, the injection wellhead pres- sures, and the magnitude, time, and location of observed induced seismic events. At the well scale, we used the well trajectory and completion information to drive calibration. We started from simple vertical wells with single cell open-hole completion to a more complex representation of horizontal wells with an open-hole completion interval across five cells in the mesh. We found that the latter representation provides a better agreement between predicted and observed gas and water rates and volumes. Moving up to the reservoir layer scale, we used the uncertainty in the layer permeability field, which is altered during hydraulic fracturing or long-term injection, to drive calibration at that scale. The original permeability is altered in a logarithmic-scale fashion mimicking the hydraulic fracturing for each producer. Next, the permeability around the two water injectors is reduced to replicate the effect of scale 13 build-up around each injector. Next, we move up to the field scale where quantities such as the initial water saturation (S wi ) are uncertain. We started with a low initial water saturation of 0.25 and compared the predicted production rates and predicted cumulative volumes to the observed rates and volumes. We increased the initial water saturation gradually until we achieved a better agreement between the predicted and observed quantities. Finally, we used the time and location data of observed seismic events to calibrate the faults vertical permeability and cohesion parameter. We used two fault failure modeling methods, namely the moment tensor method and the potency-magnitude scaling method. Chapter 3 will discuss the two approaches and the difference between them. 14 Chapter3 MathematicalModelingandGoverningEquations 3.1 Coupledmultiphaseflowandgeomechanics We consider quasi-static mechanical equilibrium of the geomechanical domain that encloses the flow do- main saturated with oil, gas and water phases. Production and injection wells are located in the flow domain. The governing equations for coupled multiphase flow and geomechanics in the geomechanical domain are obtained from the mass and momentum conservation principles [17, 41]. We assume that deformation and deformation gradients are infinitesimally small, the geomaterial is isotropic (not homo- geneous), and isothermal conditions prevail. The governing equation for linear momentum balance of the solid/fluid system can be expressed as ∇· σ +ρ b g =0, (3.1) whereσ is the Cauchy total stress tensor,g is the gravity vector,ρ b =ϕ (ρ o S o +ρ w S w +ρ g S g )+(1− ϕ )ρ s is the bulk density,ρ o ,ρ w ,ρ g ,S o ,S w , andS g are the densities and saturations of oil, water and gas phases, ρ s is the density of the solid grains, ϕ is the true or Eulerian porosity calculated as the ratio of the pore 15 volume in the deformed configuration to the bulk volume in the deformed configuration. By definition, S o +S w +S g =1. The mass conservation equation for each phaseα ∈{o,w,g} is 1 J dm α dt +∇· w α =ρ α f α , (3.2) where the accumulation term dm α /dt describes the time variation of Lagrangian fluid mass, defined as m α = JϕS α ρ α in the reference undeformed configuration, w α = ρ α v α is the phase mass flux relative to the motion of the solid skeleton,v α is the phase Darcy velocity given in terms of phase pressure p α , phase density ρ α , phase mobility k rα /µ α , and intrinsic rock permeability k, f α is the phase volumetric source term, andJ is the Jacobian of deformation given by the ratio of the bulk volume in the deformed configuration to that in the reference (undeformed) configuration. For infinitesimal deformation, J ≈ 1+ε v , whereε v is the volumetric strain. Balance equations (3.1) and (3.2) are coupled by Biot’s theory of poroelasticity [10, 23, 17]. On one hand, changes in the pore fluid pressure lead to changes in effective stress, and induce deformation of the porous material and on the other hand, deformation of the porous medium affects fluid mass content and fluid pressure. To model inelastic deformation due to rock failure, we consider the additive decomposition of the total strain into elastic and inelastic parts: ε = ε e +ε p , whereε = ∇ sym u is the infinitesimal total strain tensor given by the symmetric gradient of the displacement vectoru,ε e is the elastic strain tensor, andε p is the inelastic or failure strain tensor. Following [17], constitutive equations of multiphase poro-elastoplasticity for phaseα =o,w,g can be written in an incremental form as dσ =dσ ′ − bdp1, dσ ′ =C dr :dε e (3.3) dm α ρ α = b K dr dσ v + b 2 K dr dp+bdε p v S α + X β =o,w,g N αβ dp β (3.4) 16 where pressure p = P S β p β is the saturation-weighted equivalent pressure, σ v = (1/3)trace(σ ) is the volumetric total stress, b is the Biot coefficient of water-saturated rock, and ε p v is the volumetric strain due to shear failure. Prefix d(·) implies an infinitesimal increment in the quantity (·) which can be approximated during discrete time-stepping as the increment from the previous time step value, for example,dσ ≈ σ n+1 − σ n , where n denotes the old time level and n+1 denotes the new time level. N αβ are the inverse Biot moduli for interacting phasesα andβ and they are functions of the solid grain modulusK s , drained bulk modulusK dr , phase fluid compressibilities c α , and phase fluid saturations S α [41]. The inverse Biot moduli are components of the inverse Biot modulus tensor N = M − 1 where M is the multiphase Biot modulus tensor. The Biot coefficient b = 1− (K dr /K s ) and the inverse Biot modulus tensorN are related to experimentally-measured moduli of the fluid-saturated porous system. For a single-phase flow system, the inverse Biot modulus is N =1/M =ϕc α +(b− ϕ i )/K s , wherec α is the fluid phase compressibility, K s is the solid grain modulus, andϕ i is the initial porosity att=0. Eq. (3.4) relates the increment in fluid phase mass to the change in total volumetric stress and the change in the fluid pressures. Substituting this equation in Eq. (3.2) produces the fluid phase pressure equation. We obtain one equation for each fluid phase. The pressure equation in coupled flow-geomechanics is differ- ent from the pressure equation of flow-only models due to stress-dependent porosity and compressibility functions. Three governing equations can be rewritten in terms of three primary unknowns: gas phase pressure p g , gas saturation S g , and oil saturation S o . Corey-type relations for the gas oil capillary pres- sureP cgo (S g )=p g − p o and the oil-water capillary pressureP cow (S w )=p o − p w , saturation-dependent relative permeability model, k rα (S α ), and a pressure-dependent phase viscosity model, µ α (p α ) are used to provide closure to the mathematical system of equations in terms of the primary unknowns. For an isotropic rock, the rank-4 drained elasticity tensor can be expressed in terms of the drained modulus K dr and the shear modulus G asC dr = (K dr − (2G/3))1⊗ 1+2GI, where1 is the rank-2 identity tensor,I is the rank-4 identity tensor. The drained modulus can be extracted from the elasticity 17 tensor asK dr = (1/9)1 T :C dr :1. For the elastoplastic model used here, the elastoplasticC dr is used. The total stress and strain tensors can be split into volumetric and deviatoric parts: σ =σ v 1+S, ε=ε v 1+e, (3.5) where the volumetric total stress isσ v = (1/3)trace(σ ) and the volumetric strain is ε v = trace(ε), and trace denotes the sum of diagonal components of a tensor. The stress invariants and the principal stresses are scalar measures of a stress tensor [35], which make it easier to interpret results of a deformation experiment or simulation. We define the principal stresses as σ 1 (most compressive), σ 2 (intermediate), andσ 3 (least compressive), and the principal directions (from a suitable reference direction) asθ 1 ,θ 2 , and θ 3 . Stress invariants and principal stresses are defined in terms of each other. In 3D, σ is rank-2 and has three independent invariants. We define the first invariant of σ asI 1 =3σ v =σ 1 +σ 2 +σ 3 . The second and third invariants are denoted byI 2 andI 3 . Similarly, one can define the three invariants of the effective stress tensorσ ′ and the deviatoric stress tensorS. For example,I ′ 1 = 3σ ′ v is the first invariant of σ ′ , and J 2 = (1/2)S :S = 3σ 2 v +I 2 is the second deviatoric stress invariant, which is a commonly used stress measure in shear failure criteria for rocks. Here, we define the von Mises stress as √ J 2 and the effective mean stress asI ′ 1 , although their common definitions [35] include factors of √ 3 and1/3, respectively. A common approach to analyze how different material points in the 3D space evolve towards shear or tensile failure as loading continues is via the stress path or trajectory analysis in either (1) the stress invariant space or (2) the principal stress space. Independence from the choice of coordinate system gives the invariants and principal stresses physical meaning that is useful for two purposes: (1) understanding different deformation styles e.g., an increase in I ′ 1 can indicate production-induced volumetric contrac- tion and the consequent pull on a fault, and an increase in √ J 2 can indicate distortion and shear and an 18 increased likelihood of shear failure on a fault, and (2) developing constitutive equations based on lab ex- periments, e.g., stress invariants are used in the construction of stress-strain laws and failure criteria for the rock matrix and fault zones. In case of anisotropic principal stresses e.g.,σ 1 ̸=σ 3 , the ratio of principal stresses determine the faulting regime, e.g., normal, reverse, or strike-slip faulting on a shear failure plane, and how to orient a new well during drilling in order to create most number of hydraulic fractures or intersect most natural fractures for better productivity. Since here we are interested in both fault slip and the role of hydraulic fracturing on induced seismicity, we use both stress invariants and principal stresses to analyze our simulation results. From Eq. (3.3),dσ v = K dr dε e v − bdp. Decomposing stress into an initial stress tensorσ i and a stress change from initial, Eq. (3.1) can be written as ∇· (C dr :ε e )=b∇(∆ p)− ρ b g−∇· σ i . (3.6) which can be solved for the displacement vectoru as the primary unknown given the strain-displacement relation and a known pressure change field dp. This forms the basis of the sequential solution scheme for the coupled flow-geomechanics problem, which we use here. We use a Mohr-Coulomb shear failure model to approximate the seismic/aseismic failure processes in the rock. The details are given in the sections below. The normal displacements are prescribed to be zero on the west, east, north, and south boundaries of the 3D geomechanical domain The bottom boundary is fixed in all the three directions. The top boundary, which is at the ground surface elevation, is a traction-free boundary, and hence free to deform as stresses change within the domain. No-flow boundary conditions are prescribed on all the six sides of the flow domain, which is enclosed by the overburden, underburden and sideburden parts of the geomechanical domain. 19 We use the standard Galerkin finite element method with trilinear elements and nodal displacement vectors to discretize Eqs. (3.6) and the finite volume method with piecewise constant cell-centered pres- sures and saturations to discretize Eqs. (3.2). Solid skeleton displacementu, volumetric strain due to shear failure ε p v , fluid phase pressures p α , and fluid phase saturations S α are the unknowns of the linearized system of equations of coupled fluid flow, mechanical deformation, and fault failure. The Backward Euler method is used for time integration. We solve the two discrete systems of equations using the fixed stress sequential solution scheme. At a given time step, at a given iteration within the sequential loop, the flow problem is solved first while keeping the total volumetric stress fixed at its previous sequential iteration value. Next, the geomechanical problem is solved keeping the pressure and saturation fixed at their cur- rent iteration values. The new displacements and any incremental failure strain are then used to update the volumetric stress and porosity in the flow problem. This procedure is repeated until convergence is achieved within the sequential loop at which point the simulator marches ahead to the next time step. Each discretized sub-problem–flow and geomechanical–are solved using their respective Newton solver. We use the CMG- STARS simulator [16] to build a 3D coupled flow and geomechanics simulation model based on the numerical model presented above. 3.2 Shearfailuremodel We use the Mohr-Coulomb yield criterion to determine the onset of shear failure and a flow rule to de- termine the failure strain. We use the Mohr-Coulomb model for the yield criterion because it is based on frictional sliding of grains at the microscopic scale, which is relevant for shear failure during fault slip, and it can model the sensitivity of failure response to the volumetric stress, also known as pressure-sensitivity, 20 which cannot be modeled by criteria such as the Tresca criterion [35]. In terms of the principal stresses and the yield strength of the rock, the Mohr-Coulomb yield function is given as f =τ f − σ f sinψ f − ccosψ f , (3.7) whereτ f =(σ ′ 1 − σ ′ 3 )/2 is the maximum shear stress andσ f =(σ ′ 1 +σ ′ 3 )/2 is the mean effective stress in the material. Both are expressed in terms of the maximum effective principal stress σ ′ 1 and the minimum effective principal stress σ ′ 3 (positive in compression). c is cohesion andψ f is the internal friction angle. Whenf < 0, the material deforms elastically with a zero increment in the failure strain, i.e. ˙ ε p = 0, and when f = 0, the material deforms due to shear failure, i.e. ˙ ε p > 0. Here, the dot symbol is defined as ˙ ε=(ε(t+δt )− ε(t))/δt and hence indicates the material time derivative following the motion of the solid skeleton. The consistency condition of failure-induced deformation states that ˙ f = (∂f/∂σ ′ ) T : ˙ σ ′ = 0, when ˙ ε p >0. To determine the increment in failure strain during a loading step, we use the flow rule, ˙ ε p = ˙ λ ∂g ∂σ ′ (3.8) where the potential functiong is defined as, g =τ f − σ f sinψ g − ccosψ g . (3.9) Here,λ is the plastic multiplier andψ g is the dilatancy angle. Computing∂g/∂σ ′ gives us the expression for the incremental failure strain tensor in terms of its volumetric and deviatoric parts: dε p =dλ 2sinψ g √ 3(3− sinψ g ) 1+ 1 2 √ J 2 S ! (3.10) 21 where the first part multiplying the identity tensor is the volumetric strain due to shear failure, and the second part is the deviatoric failure strain tensor. We need the plastic multiplier,λ , to evaluate both. We obtainλ by substituting the flow rule in the consistency condition: ∂f ∂σ ′ T : C dr : ˙ ε− ˙ λ C dr : ∂g ∂σ ′ =0 (3.11) This gives ˙ λ = ∂f ∂σ ′ T :C dr : ˙ ε ∂f ∂σ ′ T :C dr : ∂g ∂σ ′ (3.12) Next, the volumetric strain due to failure can be computed by applying the trace operator to Eq. (3.8). It requires differentiating Eq. (3.9) with respect to the effective stress to obtain the plastic potential function derivative. 3.3 ModelingFaultFailureEvents 3.3.1 MomentTensorMethod Following [4, 49], we define the moment density tensor m as the difference between the linearized isotropic elastic model stress tensor and the actual physical stress tensor, also known as the ‘stress glut’. This allows the moment density tensor to be written as the product of rank-4 elasticity tensor and the incremental in- elastic strain tensor, which, in our model presented in Appendix A, are given byC dr anddε p , respectively. Then the seismic moment tensor can be modeled as M = Z V mdV = Z V C dr :dε p dV (3.13) where we define V to be the dynamically changing volume of the fault region experiencing shear failure. This model has been used in the literature [49] to obtain elastic and inelastic strain fields from geodetic 22 measurements, e.g. GPS array data, as well as to obtain the ground displacement field due to fault slip along a predefined fault surface. Assuming that the fault slip can be represented by a distribution of moment density tensorsm f localized on the fault surface and given by coordinatesx f [4], wherem = m f δ (x− x f ) andδ (·) is the Dirac delta function, we can use this model to compute the seismic moment tensor from the work done during failure. To apply this seismic model to a quasi-static geomechanics model discretized in space and time, we define three additional quantities: the failure state variable θ f , the start time of e-th shear failure event t e i , and the end time of e-th failure event t e f , which are defined for each cell. Referring to Eq. (3.7),θ f =0 whenf <0 andθ f =1 whenf ≥ 0. t e i is defined as the timestep at whichθ f jumps from 0 to 1 andt e f is defined as the timestep at which θ f jumps from 1 to 0 such that θ f =1 fort e i <t<t e f . In the discrete form, for fault cellk and for evente,M is approximated to beM e k calculated as the product ofC dr,k ,dε p k and cell volumedV k . M e k = t e f X t=t e i C dr,k (t)dε p k (t)dV k (3.14) where the time dependence of the elasticity tensor and the incremental failure strain tensor has been expressed explicitly. Each cell can undergo multiple events due to the production- and injection-induced stresses, and each event is considered as a proxy for one seismic event located at that cell and occurring at timet e f . The moment magnitude for an event can be computed asM 0 =∥M k ∥ and the seismic magnitude M w is calculated using the Kanamori equation: M w = [(logM 0 − 16.1)/1.5] + 4.667. Reducing the timestep size and choosing other methods to compute seismic magnitudes, e.g., the potency-magnitude scaling method [8], can improve the accuracy of magnitude calculation. 23 3.3.2 Potency-MagnitudeScalingMethod The Potency-magnitude scaling method was introduced by [7, 8, 6] to overcome the dependency of the Moment tensor method on the assumptions of material properties at the source, which are subject to spatiotemporal variations [7]. The relation was derived in the form logP 0 =γM 2 +cM +d (3.15) WhereP 0 is the potency,M is the magnitude, andγ ,c andd are fitting parameters (Table 3.1). P 0 can be defined by the integral of fault slip displacement l over the fault rupture areaA and is given by P 0 = Z A ldA (3.16) Assuming fault slip displacementl can be approximated by calculating the resultant difference between the displacement of the cellsu r on both sides of the fault.We calculatedu r as follows u r (I)= q (u y (I +1)− u y (I− 1)) 2 +(u z (I +1)− u z (I− 1)) 2 (3.17) A is calculated by dividing fault cells volume by the fault width. Eq. (3.16) can be rewritten for each fault cell as P 0 =u r A (3.18) The magnitude for each fault cell failure event was calculated by solving for M from Eq. (3.15) and using Table 3.1 for the fitting parameters. We used the quadratic set of parameters. 24 Table 3.1: Potency-magnitude scaling correlation parameters Linear Quadratic c 1.35± 0.03 0.988± 0.213 d -5.33± 0.06 -4.87± 0.25 γ 0.0612± 0.0385 25 Chapter4 CornerPointMeshingAlgorithm 4.1 Introduction With the increasing development in reservoir description tools and the mounting volumes of data and knowledge constantly acquired, the overall understanding of the complex structure of reservoirs is ex- panding. However, the usefulness of such understanding is often limited by our capability in honoring the structural complexity during dynamic modeling. Reservoir stratigraphy and fault geometry are examples of such structural information. Often the accuracy of dynamic simulation results depend on the structural fidelity of the model. Structured grids have been used in traditional reservoir flow simulation models, are widely accepted in all commercial simulation software, and lead to fast simulation runtimes because of the ease in building and using the grid data structures. The pillar gridding method [26, 27] which is an example of a structured grid- ding method available in commercial simulators such as Eclipse (Schlumberger) and STARS (CMG), creates rectangular cells in 2D stacked vertically along straight line pillars. This leads to zig-zag or stair-case fault geometry with sudden and geologically implausible changes in fault dip and aperture. Unstructured grid- ding techniques have been proposed as the best candidate to represent structurally complex reservoirs because they can generate tetrahedral and hexahedral meshes of 3D domains with complex internal struc- tures and geometries. However, they have not been adopted widely in the reservoir simulation community 26 due to their high cost, steep learning curve, tedious meshing steps, and incompatibility with many com- mercial simulation software. As a result, we must find a workflow that balances the modeling accuracy with the computing cost. Usually, the balance has been found by sacrificing properties such as the fault geometry in order to represent the reservoir in a structured pillar meshing format. STARS offers the op- tion of building a corner point grid with trapezoidal cells in order to better represent more complex and faulted reservoirs. However, this option requires the knowledge of the eight corner points of each cell in the grid. Pillar gridding may be acceptable for structures which have all its internal and boundary surfaces aligned with one of the coordinate directions e.g. reservoirs containing no faults or pinch-outs or vertical faults with uniform dip. However, the increasing need for coupled flow and geomechanical simulations has made the following three points very clear: (1) all faults, including the ones that are geologically less certain and hydraulically less important, must be represented in the geomechanical model to accurately model the transmission of induced stresses, (2) fault curvature and intersections must be honored because the fault stresses depend on the fault geometry by definition, and (3) faults must be extended above and be- low the production or injection interval i.e. in the overburden and basement to allow modeling of induced seismicity in basement and mechanical integrity failure in caprock. We developed a corner point structured meshing algorithm that can model a geologically complex reservoir with multiple, curved and intersecting faults for high-resolution coupled multiphase flow-geomechanical simulations using commercial simulators. Using the routine, we build a dual flow-geomechanical grid in CMG STARS and demonstrate its accuracy in modeling flow-induced stresses and fault stability in a cou- pled simulation. Our meshing routine incorporates the reservoir geometry and the 3D shape of curved, connected faults to generate corner point grid coordinates in the input format accepted by STARS. 27 Figure 4.1: A schematic (not to scale) showing a horst and graben region. 4.2 PhysicalModel In order to examine our meshing algorithm, we elected to simulate a case with multiple levels of com- plexities, such as multiple curved and intersecting faults. We consider a horizontally-layered reservoir in a Horst and Graben region, as shown in Figure 4.1 [24], confined between overburden and underburden and bounded by synthetic-antithetic faults. The reservoir is located in a normal faulting stress regime, where the downthrown side of the reservoir is located in the hanging wall block of the two faults and the upthrown sides are located in the footwall blocks (Figure 4.2) . Figure 4.3 shows the exact proportions of the physical model. A well, located in the center of the reservoir, produces oil and water at a constant flow rate subject to a minimum bottomhole pressure constraint. We assume that the two faults are seal- ing (non-conductive), and the reservoir is hydraulically isolated from the adjacent formations due to low permeability of the fault and the overburden-underburden layers. Boundary conditions on the geomechanical domain are as follows. The normal displacements are prescribed to be zero on the left, right, front, and back boundaries. The bottom boundary is fixed in all the three directions. The top boundary, which is at the ground surface elevation, is a traction-free boundary, and hence free to deform as stresses change within the domain. No-flow boundary conditions 28 Figure 4.2: A schematic of the physical model showing the reservoir in yellow, synthetic fault in red and antithetic fault in green and the overall geomechanical domain Figure 4.3: Cross sectional views of the domain in theX− Z andY − Z planes showing the proportional dimensions of the reservoir and the faults compared to the geomechanical domain. It also shows the curvature of the faults and the semi trapezoidal geometry of the reservoir. 29 Table 4.1: Meshing Model Initial Conditions Property Value DatumPressure(psi) 4000 DatumDepth(ft) 5741.5 VerticalTotalStressGradient(psi/ft) 1 MaximumHorizontalTotalStressGradient(psi/ft) 0.8 MinimumHorizontalTotalStressGradient(psi/ft) 0.7 Figure 4.4: Cross sectional views of the domain in the X-Z and Y -Z planes showing the mechanical boundary conditions. The reservoir is bounded by no-flow boundaries from all six sides. are prescribed on all the six sides of the reservoir (Figure 4.4). Table 4.1 lists the initial conditions for the model. A Mohr-Columb elastoplastic failure model was used for both the matrix and the fault regions (Ta- ble 4.2). The fault in the simulation model is modeled as a finite thickness region with one cell across its aperture. The model is assumed to be saturated with oil and water phases. Initial pressure is distributed using the hydrostatic law and a datum pressure at a datum depth. Table 4.2: Meshing Model Poroelastic Properties Property Overburden Reservoir Underburned Faults Porosity 0.25 0.25 0.25 0.25 Permeability(md) .01 10 .01 .01 Youngmodulus(psi) 3.0E+4 3.0E+4 3.0E+4 3.0E+4 Poissonratio 0.3 0.45 0.35 0.2 Cohesion(psi) 50000 50000 50000 500 Frictionangle(O) 45 45 45 10 30 4.3 Approach In this section, we explain the workflow (Figure 4.5) that we followed to generate a dual geomechanical and flow mesh. The first step is to define the spatial parameters of the geomechanical domain mesh. This includes theX,Y , andZ dimensions, the desired number of cells in the three dimensions, and∆ Z (thickness) and∆ Y dimensions of the cells.∆ X of elements are calculated internally. Next, we define the number of faults, fault starting and ending points in the domain, and fault geometry including its curvature, aperture and the dip angle. We define fault positions on the front boundary ( − Y ), back boundary (+Y ), and interior planes where the strike changes. Our framework allows us to specify different fault properties on front and back faces which allows the fault curvature to change along theY -axis. We also allow flexibility in specifying starting and ending points of the fault in the domain. These allow a greater freedom in modeling geologically complex domains. To honor the continuity requirement along the logical meshing directions in a structured mesh, we introduce an imaginary fault that starts at the bottom of the minor (antithetic) fault, which is at a short horizontal distance from the major (synthetic) fault it intersects, and extends parallel to the major fault and connects with the bottom boundary of the domain. Physical properties of the imaginary fault are set equal to those of the host rock around it to honor the mechanical integrity of the model. The logical meshing directions, I, J and K in our structured mesh are quasi-parallel to the Cartesian directions, X, Y , and Z. The origin of the mesh is the top left front corner of the 3D domain. Next, we divide the domain horizontally into a series of fault block compartments and faults. Each compartment is further sub-divided uniformly into mesh elements with element∆ X calculated from the compartment dimensions. The X, Y and Z coordinates of the layers in-between are calculated through linear interpolation based on the Y and Z dimensions. At the end of this step, we have the (X,Y,Z) coordinates of each corner point in the 3D mesh. 31 Figure 4.5: Dual Geomechanics and Flow Coordinates Meshing Routine Workflow. 32 Figure 4.6: Model Coordinate system and the point of origin. This cube shows all six sides, Front (F), Back (B), Right (R), Left (L), Up (U) and Down (D). The next step is to identify the dimensions of the flow domain and its position inside the geomechanical domain. This is achieved by utilizing the corner point coordinates. The flow domain is uniformly meshed except for some locally refined regions near wells, bed boundaries, and faults to obtain better resolution in pressure and saturation results in those regions. The proposed meshing method generates a simulation-ready, 3D structured mesh that honors the curvature and thickness of multiple intersecting faults and the stratigraphy of the overburden-reservoir- underburden system shown in Figure 4.8. This is significant because of three reasons. First, the fault curvature determines the orientation of the fault normal vector, which is used to calculate the fault trac- tions from the stress tensor in the host rock. An accurate calculation of tractions is required to determine the onset time, location and magnitude of failure events on the fault. Second, the fault thickness determines the hydraulic properties and pressure propagation along and across the fault, which is the driving force behind induced stresses and mechanical failure. Third, by robustly generating 3D meshes with multiple in- tersecting and curved surfaces, the proposed method addresses a well-known challenge in computational 33 Figure 4.7: Corner point ordering algorithm (STARS User Guide, 2016). mechanics. Figure 4.9 shows a comparison between the proposed corner point geometry (CPG) meshing method and the two most common structured meshing methods, classical and pillar. 4.4 ResultandDiscussion In order to demonstrate the robustness and accuracy of our corner point approach, we generated another model using the typical pillar meshing method. Both models have the same flow and petrophysical prop- erties. The fault in the pillar mesh is represented as a “stair-case” object with a one-to-two cell thick 34 Figure 4.8: Corner point representation resulting from our proposed meshing routine. The reservoir is represented by the flow mesh and the non-reservoir domain (overburden, underburden, and sideburden) is represented by the geomechanical mesh surrounding the reservoir. Coupling between flow and geome- chanics is honored within the reservoir mesh. (a) (c) (b) Y X Z Figure 4.9: Curved fault representations by different meshing techniques showing the accuracy in de- scribing curved faults and the principal stress acting on the surfaces of individual fault mesh elements. (a) Classical structured meshing technique. (b) Pillar meshing technique. (c) Corner-point meshing technique. 35 Figure 4.10: Poisson’s ratio distribution in (a) the Pillar Mesh and (b) the Corner Point Mesh. Young’s Modulus in (c) the Pillar Mesh and (d) the Corner Point Mesh. aperture in the X-direction, thus resulting in an approximate representation that may be unrealistic of typical faults. In the corner point mesh, the fault is represented as a one cell thick object that perfectly honors the aperture and curvature of the fault. Figure 4.10 shows the Poisson ratio and Young modulus distribution for both meshes. Figure 4.11 shows the permeability distribution for both meshes. For the pillar mesh, the permeability is defined for the whole grid. The permeability for the corner point mesh is defined over the reservoir flow domain only. This is because no flow calculations are carried in the non-reservoir region of the dual mesh. Figure 4.12 shows the pressure profile for the two meshes at different time steps. The first time step when there is a substantial pressure difference in the pillar mesh is at t = 1200 day. For the corner point 36 Figure 4.11: Permeability (md) distribution in (A) the Pillar mesh and (B) the Corner point mesh. mesh, it is at t = 2000 day. This can be attributed to the slight volume difference between the two rep- resentations of the reservoir in the two meshes and the resulting spatial distortion in the pillar mesh. The pressure depletion follows a uniform pattern in both meshes as it follows the uniform permeability distribution. The vertical displacement figure 4.13 shows that the reservoir undergoes a compaction process during the pressure depletion, thus resulting in a subsidence on the top boundary which is free to deform. We also notice that the reservoir exerts pulling motion on the underburden which is arrested by the fixed bottom boundary. Figure 4.14 shows the evolution of the Von Mises stress and the effective mean stress in a fault cell on the synthetic fault adjacent to the production well perforation. The initial values of the two stresses are different (smaller for the pillar mesh cell) for the two meshes because stress calculation depends on the cell geometry (fault dip) and size (fault aperture). Mesh-induced difference in pressure (Figure 11) also causes a difference in the stresses. To evaluate the effect of meshing on the dynamic response (fault stability), we shift all the stress values of the pillar mesh to match with the corner mesh values at t = 0. The two different stress paths of the similarly-positioned fault cell in the two meshes suggest a difference in the 37 psi Figure 4.12: Pressure (psi) profile for the two meshes at different time steps. (a) Pressure profile for pillar mesh at t = 1200 days (b) at t = 10000 days. (c) Pressure profile in the flow grid only for corner point mesh at t = 2000 day and (d) Pressure at t = 10000 day 38 ft ft Figure 4.13: Vertical Displacement (ft.) for the two meshes at different time steps. (a) Vertical displacement for the pillar mesh at t = 1200 days (b) at t = 10000 days. (c) Vertical displacement for the corner point mesh at t = 2000 days and (d) at t = 10000 days. 39 Figure 4.14: Stress path in the failure space of Von Mises stress vs. the effective mean stress. The failure envelope is shown as a thick dash line. Stress paths of the fault cells for the two mesh types are shown for different time steps. Green triangles represent the stress path for the pillar mesh and the purple triangles represent the extrapolated failure point. Red circles represent the stress path for the corner point mesh and the blue circle represent the extrapolated failure point. Approximating the curved fault geometry (corner point) with a stair-stepped profile (pillar) leads to inaccurate stresses and fault stability predictions in case of the pillar mesh. To allow comparison between the two stress trajectories, the pillar mesh stress values are shifted to match the corner point stress values at t = 0. induced failure times for the synthetic fault. This is shown by the extrapolated failure points corresponding to the two meshes. This demonstrates the effect of meshing accuracy on the model accuracy. 40 Chapter5 SyntheticModel 5.1 PhysicalModel We consider a horizontally-layered reservoir confined between overburden and underburden and inter- sected by a planar fault. We consider a normal faulting stress regime, where the downthrown side of the reservoir is located in the hanging wall block of the fault, and the upthrown side is located in the footwall block of the fault. A well, located in either the footwall reservoir or the hanging wall reservoir, produces oil at a constant flow rate subject to a minimum bottomhole pressure constraint. We assume that the two sides of reservoir are not in hydraulic communication due to low permeability of the fault and overburden-underburden layers. Figure 5.1 (a)shows the seismological model of the fault for stress inver- sion, and Figure 5.1 (b)shows the poro-elastoplastic model for coupled simulation. The objective of the seismological model is to reproduce the focal mechanisms of the seismic events caused by the fault move- ment. The objective of the poro-elastoplastic model is to compute the changes in pressure and stresses in the reservoir and track the evolution of plastic failure within the fault. The poro-elastoplastic simulation model is defined in a quasi-3D domain (Table 5.1). It is a vertical XZ section with one cell in the Y direction. Figure 5.1 (b) shows the initial and boundary conditions of the domain which are as follows. The normal displacement is prescribed to be zero on the left, right, front, and 41 Z X a) b) Figure 5.1: (a) Seismological model representing a normal fault with specified strike, rake, and dip an- gles. (b) The poroelastic model shows the reservoir in the footwall block (Reservoir (A)) and the hanging wall block (Reservoir (B)). The overburden and underburden layers, shown with a hashed region, have a permeability value lower than that of the reservoir. The mechanical boundary conditions are shown. All boundaries are no-flow boundaries. A well is located in either the footwall reservoir (Case 1) or the hanging wall reservoir (Case 2). Table 5.1: Synthetic Seismological Model Fault Parameters FaultProperty Value Strike 0 ◦ Rake -90 ◦ Dip 10 ◦ 42 Table 5.2: Synthetic Model Initial Conditions Parameter Value Initial effective vertical stress 9000 Initial effective horizontal stresses in X and Y directions 6000 Initial datum pressure at the top 3000 Well production rate, minimum BHP 1000, 100 Table 5.3: Synthetic Model Matrix Poroelastic Properties Matrix Property Value Reservoir Permeability 10 Overburden-underburden Permeability 0.01 Porosity 0.25 Model Size 5000 x 500 x 500 Number of cells in X,Y,Z 51 x 1 x 50 Young Modulus 3.0E4 Poisson Ratio 0.45 Cohesion in elastoplastic model 50000 back boundaries. The bottom boundary is fixed in all the three directions. The top boundaryis a traction- free boundary, and it is free to deform as stresses change within the domain. No-flow boundary conditions are prescribed on all the six sides of the model. Initial conditions are listed in Table 5.2. The horizontal tectonic stresses in the X and Y directions are equal and the vertical effective stress is 1.5 times the hor- izontal effective stress. This results into a normal faulting stress regime with maximum principal stress oriented vertically. The stress ratio of 1.5 comes from the stress inversion results. Initial displacements in the model are zero and satisfy the mechanical equilibrium. We use a Mohr-Coulomb elastoplastic failure model for both the matrix and the fault regions (Table 5.3 and Table 5.4). The fault is modeled as a finite thickness region with one cell across its aperture. The model is assumed to be saturated completely with oil. Pressure in the model is distributed using hydrostatics and a datum pressure at a datum depth.s In order to initiate a fault movement, the hydrostatic and mechanical equilibrium need to be disturbed. That will be demonstrated by two cases: the first one is through pressure depletion by the oil producer located in Reservoir-A and the second one is through a producer located in Reservoir-B (Figure 5.1 (b)). 43 Table 5.4: Synthetic Model Fault Poroelastic Properties Fault Property Value Aperture 2 Along- and across-fault permeability 0.01 Porosity 0.25 Young Modulus 3.0E4 Poisson Ratio 0.2 Cohesion in the elastoplastic model 500 Friction Angle in the elastoplastic model 10 Figure 5.2: The seismological model showing three seismic events with the respective focal mechanisms on the fault plane. 5.2 Approach We consider a sequence of seismic events occurring on the fault below the reservoir. Due to the complexity of real faults and tectonic stress regimes, and the errors involved in seismic data interpretation, we generate synthetic seismic events. The locations of these seismic events are shown in Figure 5.2. In order to have a representative sample of the real data, we generate 1500 focal mechanisms, 500 at each of the 3 locations, based on Gaussian distributions assumed for each of the three parameters (strike, dip, rake). Figure 5.3 (a), (b), and (c) shows the distributions of the three parameters. FigureFigure 5.3 (d) shows a sample of the 44 -40 -20 0 20 40 -10 10 30 -130 -110 -90 -70 -50 a) b) c) d) Figure 5.3: (a) Gausian distributions of strike, dip, and rake angles. (b) sample of synthetic focal mecha- nisms data used for stress inversion. focal mechanisms generated through the synthetization process. All 1500 data points are utilized during the stress inversion to produce stress field orientations and relative magnitudes. We use these results from stress inversion as inputs to our poro-elastoplastic simulation model. We examine two simulation cases: Case 1 with the producer located in Reservoir (A) (Figure 5.4 (A)), and Case 2 with the producer located in Reservoir (B) (Figure 5.4 (B)). The two cases show the effect of well location on fault failure and subsequent seismicity. The simulation end time is 10,000 days. 5.3 ResultandDiscussion 5.3.1 FormalStressInversion The results of the formal stress inversion are presented in Figure 5.5, which shows the orientations of the three principal stresses locally around each of the three hypocenters. σ 1 represents the maximum principal stress and σ 3 represents the minimum principal stress. The steronet plots for the three hypocenters show thatσ 1 is closer to the center of the steronet, which denotes that the maximum stress is approximately vertical ( 30 degrees from vertical). The angle between sigma 1 and sigma 3 is 90 degrees. This distribution of principal directions agrees with a normal faulting stress 45 a) b) Figure 5.4: Permeability distribution and well location in the poro-elastoplastic model. (a) Case 1 and (b) Case 2 correspond to two different simulation models with two different well locations. (a) (c) (b) Figure 5.5: Stress inversion results plotted using stereonets. Directions of principal stresses can be read from the locations of the red, green, and blue points with respect to the center which represents the ver- tical direction and the perimeter which represents the horizontal direction. The difference between the maximum and minimum principal stress directions for the 1500 events is 90 degree 46 3202 159 1681 a) b) Figure 5.6: Pressure at t=1428.6 and t=10000 day for Case 1. regime. In addition to the stress orientations, stress inversion also provides us with relative stress magni- tudes of the three principal stresses at the hypocenter of the seismic event. In order to honor the seismicity before well activities, we use these stress inversion results to initialize stresses in our coupled flow and geomechanical simulator. 5.3.2 CoupledFlowandGeomechanicsSimulation 5.3.2.1 Case1 Figure 5.6 shows the pressure map of the reservoir before the fault failure begins and at the end of simu- lation after pressure depletion in the reservoir. Since we have a uniform permeability distribution in the reservoir, the pressure depletion is also uniform. Figure 5.7 present the effective mean stress. It shows that the effective mean stress is increasing in the reservoir (compressive stresses are positive, tensile stresses are negative) and follows the distribution of 47 a) b) 9606 6927 8267 Figure 5.7: Effective Mean Stress at t=1428.6 and t=10000 day for Case 1. Compression is positive. pressure depletion. The effective mean stress is related to the total mean stress and the pore pressure as follows: σ ′ v =σ v − αp (5.1) whereσ v is the total mean stress,σ v prime is the effective mean stress, p is the pore pressure, andα is the Biot coefficient. At t = 0, the effective mean stress is given by the average of the initial effective stresses i.e. (9000+6000+6000)/3 = 7000 psi. Figure 5.8 presents the effective plastic stain before fault failure and at the end of simulation. The fault fails plastically under high von Mises stress due to production-induced reservoir contraction, which causes changes in the total stress and an increase in the deviatoric stresses. Although the increase in effective mean compression tends to stabilize the fault, the increase in von Mises stress is large enough to induce plastic failure and cause positive plastic strains in the fault. Moreover, using results at multiple 48 a) b) 0.056 0 0.028 Figure 5.8: Effective Plastic Strain at t=1428.6 and t=10000 for Case 1. time steps, we conclude that the failure propagates both updip and downdip along the fault. The failure span is longer above the reservoir because of the traction-free top boundary which allows larger displace- ments above the reservoir than below the reservoir. However, the effective plastic strain, a non-negative value, does not indicate the direction of the slip, which is given by the relative displacement across the two sides of the fault. Therefore, we plot the vertical displacement (Figure 5.9) to understand the fault slip direction. It shows that the reservoir undergoes a compaction process during pressure depletion. Since the top boundary is free to deform, and the bottom boundary is fixed, the vertical displacement is highest near the top boundary. A larger downward displacement in the footwall block and a smaller downward displacement in the hanging wall block indicates a reverse motion on the fault. Figure 5.10 shows the von Mises stress before fault failure and at the end of simulation. The von Mises stress quantifies the deviatoric stresses and is used to predict plastic yielding of materials under loading. The von Mises stress becomes high around the perforation and becomes very high within the fault. This 49 a) b) 0.5 -5.4 -2.5 Figure 5.9: Vertical Displacement at t=1428.6 and t=10000 for Case 1. Unit is feet. 2863 4288 (a) (b) Figure 5.10: von Mises stress at t=1428.6 and t=10000 for Case 1. 50 is a result of reservoir contraction, which applies downdip shear on the fault above the well perforation depth and updip shear on the fault below the well perforation depth. Reservoir contraction also applies a tensile pull on the fault, however the net change in fault normal stresses is compressive because the depletion-induced increase in the effective mean stress (Eq. (5.1)) overwhelms the contraction-induced pull. Figure 5.11 shows the evolution of the von Mises stress and the effective mean stress for four fault cells located at the top of the reservoir, middle of the reservoir, bottom of the reservoir, and adjacent to the well perforation. Prior to fault failure, all points are below the failure envelope which is a straight line given by the cohesion and friction angle of the fault failure model. As the cell pressure decrease with time, the von Mises stress and the effective mean stress increase until the point hits the failure envelope at which time plastic failure is initiated and a positive plastic strain is computed. After failure is initiated, the points continue to move along the failure envelope indicating that the rate of depletion is large enough to prevent the von Mises stress from relaxing. Figure 5.12 shows the effective plastic strain in the fault. It provides us with an understanding of how failure nucleates and propagates on the fault. The failure nucleates at the point adjacent to the well perfora- tion where pressure change is the highest. It propagates both updip and downdip with an asymmetry that is driven by the assymetry in the vertical displacement field caused by asymmetric boundary conditions on the top and bottom boundaries. The moment magnitude,Mw, of an equivalent seismic event (equivalent in terms of the total mechan- ical energy released) can be computed from the failure area and the cumulative plastic strain [29, 36]: G is the shear modulus in the fault. d is the cumulative slip magnitude calculated from the plastic strain and the fault cell size.Figure 5.13 presents the time evolution of failure, which shows a trigger time att = 1428.6 day and a final Mw of 2.3. 51 Figure 5.11: Stress path in the failure space of von Mises stress vs. effective mean stress. The failure envelope is shown as a dash line. Stress paths of four fault cells are shown by plotting the stresses at five timesteps (value of t indicates the timestep number, not actual time in days). Each time step has a unique symbol. At the first timestep (t=1), all four cells have the same von Mises stress (3000 psi) and the effective mean stress (7000 psi) because the initial stresses are uniform in the model (see Table 2). At t=2, the cell next to well perforation moves, and other cells stay at the initial point 52 Figure 5.12: Depth profiles of effective plastic strain at different time steps (time in days shown in the legend) during the simulation. The nucleation point (hypocenter) is at the well perforation depth and failure propagates both updip and downdip, with a larger updip span because of the traction-free top boundary. Figure 5.13: Total Effective Plastic Strain for Case 1 53 (a) (b) 3,202 195 psi Figure 5.14: Pressure at t=1666 and t=10000 for Case 2. Unit is psi 5.3.2.2 Case2 Compared to Case 1, the oil producer is located in the hanging wall block of the fault. We show similar plots as Case 1 to facilitate comparison and analyze the effect of well location. Figure 5.14 shows pressure maps of the reservoir before fault failure and at the end of simulation. Figure 5.15 presents the effective mean stress, which increases due to pressure depletion. Figure 5.16 presents the effective plastic strain before and after fault failure. Failure nucleates around the depth of 2,360 ft. adjacent to the well perforation and propagates updip until it hits the top boundary of the reservoir on the hanging wall side. The downdip propagation of failure is arrested by the fixed bottom boundary. Figure 5.17 shows the vertical displacement in the domain as a result of pressure depletion. It indicates a normal slip motion on the fault because the hanging wall block moves downward by a larger magnitude than the footwall block. Figure 5.18 presents the von Mises stress before fault failure and at the end of simulation. 54 (a) (b) 9,658 6,931 psi Figure 5.15: Effective Mean Stress at t=1666 andt=10000 day for Case 2. Unit is psi (a) (b) 0.060 0 Figure 5.16: Effective Plastic Strain at t=1666 and t=10000 day for Case 2 55 (a) (b) 0 -6.1 Figure 5.17: TVertical Displacement at t=1666 day and t=10000 day for Case 2. Figure 5.20 shows the effective plastic strain in the fault and Figure 5.21 presents the time evolution of failure, which shows a trigger time at t = 3333.3 day and a final Mw of 2.3 that is similar to Mw of Case 1. 56 (a) (b) 4.297 2.854 Figure 5.18: Von Mises Stress at t=1666 and t=10000 for Case 2. Figure 5.19: Stress path in the failure space of von Mises stress vs. effective mean stress. The failure envelope is shown as a dash line. Stress paths of four fault cells are shown by plotting the stresses at five timesteps (value of t indicates the timestep number, not actual time in days). Each time step has a unique symbol. 57 Figure 5.20: TDepth profiles of the effective plastic strain on the fault at different time steps. Legend shows time in days. Figure 5.21: Total Effective Plastic Strain for Case 2. 58 Chapter6 AzleInducedSeismicEvents 6.1 Introduction Induced seismicity and aseismic deformation events around subsurface reservoirs are often associated exclusively with either injection or production of fluids. We know that production-induced stresses pro- mote reverse faulting and injection-induced stresses promote normal faulting on basement faults [57, 58], and fluid diffusion into faults lead to their unclamping due to a drop in effective compression [34, 39, 62]. However, in case of simultaneous injection and production in the same reservoir, e.g., in joint CO 2 injection–brine production operations [9, 14] and in oilfields with produced water disposal wells, it can be difficult to extract the dominant causative mechanism for such deformation events and the individual roles played by production and injection in modulating that mechanism. Understanding and predicting the stress propagation from injection and production locations onto basement faults is essential to forecast and control the induced seismicity hazard at such extraction–disposal sites. This requires understanding how production- and injection-induced stresses interact to destabilize a fault intersecting or bounding the production and injection regions. This is further complicated by production and injection intervals that are stratigraphically and hy- draulically segregated, i.e. they occur in different layers or fault blocks of a single reservoir, or in two 59 different reservoirs within the same basin. Such is the case of Barnett Shale production–Ellenburger in- jection in the Newark East field located in the Fort Worth basin of Texas, which was the epicenter of 2013-2014 Azle-Reno earthquake sequence with a maximum magnitude of Mw 3.8 [32] (Figure 6.1). The seismic events were located on a pair of synthetic–antithetic normal faults crossing through Barnett and Ellenburger and extending into the basement. Gas and brine producing Barnett wells operate on both sides of the synthetic and antithetic faults. Two water disposal wells in Ellenburger operated in close vicinity of the hypocenters in the footwall block of the antithetic fault, although there are more injectors located farther away from the faults. Studies have been conducted to explore the causality of these seismic events in terms of Ellenburger injection [32, 33] and Barnett production [15] using a single-phase pressure dif- fusion model and a coupled multiphase flow-geomechanics model, respectively. These studies have been crucial in establishing the anthropogenic nature of the seismic events. However, they also raised impor- tant questions that need further investigation. In particular, there are three questions related to (1) the role of Barnett production, both gas and water, and Ellenburger injection in inducing the events, (2) the spatiotemporal propagation of production- and injection-induced stresses from wells to the hypocentral locations on the faults, and (3) the role of model’s resolution in space and time on the prediction accuracy. Below we discuss these issues in more detail. 6.1.1 Roleofwastewaterinjectionandwaterproduction Two of the previous studies [32, 33] concluded that Ellenburger injection-induced pressurization of the faults played a critical role in inducing the Azle events. A third study [15] concluded that the injection pressure did not propagate into the fault, and the events were caused indirectly by the difference in poroe- lastic stresses on the two sides of the fault, with gas production playing the dominant role. Hence, further investigation is warranted. Also, the role of water production is not clear in prior studies. 60 Knowing whether the water produced from Barnett wells is the Barnett pore water, or the flowback water that was injected during hydraulic fracturing of Barnett wells, or the Ellenburger water flowing upwards via vertical fractures connecting the two formations, or a combination of these, will be helpful in modeling induced stresses and subsequent fault reactivation. Since water is less compressible than gas, for the same volume produced, water can cause a larger change in pressure and stress than gas. Previous studies assumed that all the water produced from Barnett wells came from Ellenburger through fractures. However, this assumption has not been fully supported by any fracture height data or geochemical analysis of the produced water that can tie the source of water to Ellenburger depths. Pressure-rate analysis of the Barnett well data, required to rule out the possibility of water production from Barnett, is also missing from public domain. Therefore, assuming a widespread presence of vertical fractures providing a conduit for the Ellenburger water to flow out of Barnett wells may require further investigation [51, 33]. This is relevant for modeling induced seismicity in Azle because making that assumption of water pro- duction exclusively from Ellenburger implies the net overpressure in Ellenburger due to water disposal and Barnett’s production may not be sufficient, or solely responsible, for Azle’s induced seismicity. This ob- servation is supported by earlier studies, e.g. a maximum excess pressure of 0.008 MPa in Ellenburger [32] seems insufficient to induce a Mw 3.8 event. This is also supported by the measured overpressure values in Ellenburger, which are much higher and in the range of 1.7 to 4.5 MPa [33]. Another argument against the assumption of Barnett’s produced water coming from Ellenburger is that, before the seismicity, pres- sures in many Barnett wells were on a declining trend for a long time, which indicates a lack of aquifer pressure support from Ellenburger to Barnett, which implies poor hydraulic connectivity between the two. Here, we hypothesize that a fraction of the produced water is the flowback liquid that was injected during hydraulic fracturing of Barnet wells. We test the hypothesis in this paper. 61 Marble Falls Basement Caprock Antithetic fault Synthetic fault Water injection Gas & water production Gas & water Production Hyd. Fractures a X X’ Ellenburger Barnett b Newark East Fault Zone Injection well X X’ 32.92 33.0 Texas Oklahoma New Mexico Azle * –97.60 –97.52 Antithetic fault Synthetic fault Injection well 2 km diusion Natural Fractures Figure 6.1: (a) A map view of the Azle region (marked by red star in the inset map) at Ellenburger depth showing 2013-2014 events and the synthetic–antithetic fault pair. Re-drawn with permission from [32]. (b) Stratigraphic section of the Newark East oilfield along X-X’ line in (a), showing Barnett production, Ellenburger injection, and the fault pair. The figure is not to scale. The stars indicate the general locations of hypocenters during induced seismicity. Interaction between contraction of Barnett and expansion of Ellenburger, shown by red circles with inward and outward pointing arrows, and injection-induced pres- sure diffusion into basement can generate induced stresses on the faults, which will vary laterally across fault blocks and vertically along the faults. 6.1.2 Interactionandpropagationofinducedstresses For the general case of production and injection from the same reservoir layer, the interaction between production- and injection-induced stresses and their propagation from production and injection depths to basement faults has been discussed. This is not the case in Barnett, where production and injection reservoirs are stratigraphically segregated and, due to a large difference in permeability (microdarcy in Barnett vs. 30-50 millidarcy in Ellenburger), also hydraulically segregated over time scales of few years. We hypothesize that the fault stability in the segregated case cannot be determined using pressure diffusion models uncoupled to geomechanics, and poroelastic stresses must play a dominant role. We also hypothesize that coupled flow–geomechanics–fault failure models, which capture the geologic structure and well flow at higher spatiotemporal resolutions than afforded by older models, are necessary to compute the induced stresses and the timing and location of failure events accurately. There are three observations that support this hypothesis. (1) An accurate determination of poroelastic stresses in injection 62 and production layers requires a model that includes both the layers as well as the caprock, basement, and sideburden, because the domain of geomechanical equilibrium includes all these regions [58, 2, 43, 19, 37, 36]. In other words, it is not accurate to isolate a hydraulic unit to study induced stresses because the mechanical boundary conditions on the isolated unit is unknown and most likely changing in time due to the mechanical action of adjacent hydraulic units. (2) The interaction between production- and injection-induced stress tensor fields is different from interactions between two scalar fields of pressure [35]. This is the reason one cannot estimate the change in maximum principal stress direction and magnitude by looking at the changes in pressure in production and injection layers, even when the driving force behind the change in stress state is the change in pressure field. (3) The 3D mesh of the coupled model must resolve the fault geometry, including changes in fault dip and strike and intersections with stratigraphic horizons, to accurately project the induced stress tensor as shear and normal traction vectors on the fault surface. This is true because the tractions are required to calculate the Coulomb stress and the seismicity rate [58]. Here, we present novel results on the spatiotemporal distribution of principal stresses in production/injection layers and on the two faults shed new light on the role of Barnett production on fault stability in the Azle region and understanding the causality of induced earthquakes in general. 6.1.3 Spatiotemporalresolutioninmultiphysicsmodels Recent observations of production- and injection-induced seismicity around aquifers [25, 32, 2], gas storage reservoirs [22], and petroleum reservoirs [55, 37] has motivated the development of advanced multiphysics models that can honor the geologic structure of a faulted caprock–reservoir–basement complex and the multiphysics processes activated during production and injection in such a complex. In case of Barnett, almost all of its hundreds of producers are horizontal and hydraulically fractured. Therefore, to account for the effect of hydraulic fracturing on the spatial distribution of pressure, the permeability field in the 63 multiphysics model should resolve the effect of fracturing. Number of wells used in the model also affects how resolved the induced stresses are in space and time. For example, prior studies [15] considered 70 Barnett producers for a 12 year production period in a 10 km× 10 km area centered around the recorded hypocenters. However, Barnett has hundreds of producing wells spread over a larger area that are likely affecting its pressure and stress fields. Here, we consider all ∼ 350 wells within a 25 km× 15 km area around the faults, including three additional far-away injectors, to understand the effect of high-resolution well data on the spatiotemporal resolution of model predicted stress fields. The material constitutive model used to capture rock’s deformation response, especially in presence of time-variable fluid withdrawal and injection rates in wells affects the spatial distribution and temporal evolution of induced stresses, and, hence, model’s resolution. Barnett shale shows inelastic deformation behavior [59]. Therefore, to account for this effect on model resolution, we consider plastic failure in the proposed model. In summary, previous studies on induced seismicity modeling in Azle call for more research on the topic by including more data, higher resolution, and coupled physics in the model. Here, we answer that call by conducting a modeling study that relaxes some of the assumptions made in prior studies, includes data from more wells in the region, and analyzes the coupled processes at resolutions higher than reported before. We hypothesize that increasing the space and time resolution of geologic objects and well flow data that drive such advanced models can help us understand the propagation of induced stresses from the location of production/injection, e.g. a well in a reservoir layer, to the location of an induced event, e.g. an earthquake on a basement fault, which may be hydraulically isolated from the reservoir. Here, we test this hypothesis by building a high-resolution, coupled multiphase flow-geomechanics-rock failure model of the Azle region that uses daily production/injection rates from∼ 350 wells (sources: Texas Railroad Commission, DrillingInfo, and IHS Energy) and seismic data from USGS and SMU stations [32, 33]. 64 6.2 Geologicstructure,fluidvolumesandwellconfiguration The Newark East field is located in the north Texas about 50 Miles west of Dallas. Geological structure of our selected area of study is given by a∼ 100 m-thick Barnett shale overlying the∼ 1 km-thick Ellenburger dolomite formation, which overlies a thick granite basement. The Barnett shale is a Mississipian-aged black, organic-rich shale at depths of about 1980 m to 2590 m. It is an unconventional petroleum resource because the source, seal, and reservoir rocks are all the same and the rock has ultra-low matrix permeabil- ity in the range of approximately 0.01 to 1 microdarcy. Hydraulic fracturing is required for the wells to flow at commercially viable rates. Barnett’s overburden is composed of Marble Falls limestone and shaly clastic/carbonate layers above that. The region is dominated by normal faulting [31] with the maximum principal stress being vertical. The Azle-Reno earthquake sequence is hosted by a synthetic-antithetic nor- mal fault pair, analogous to a horst-graben structure. The synthetic fault, Newark East Fault Zone (NEFZ) dipping 60-70 ◦ and towards northwest, extends through Barnett and Ellenburger, and the antithetic fault, which dips 70-80 ◦ and towards southeast, is smaller in extent and terminates at the synthetic fault (Fig- ure 6.2a) at approximately 3500 m depth. NEFZ consists of two fault surfaces of different strike angles in the general northeast-southwest direction. The original gas-in-place (OGIP) in Barnett has been estimated to vary between 120 and 200 billion standard cubic feet per square mile [20] depending on the reservoir thickness, organic content, porosity, and gas compressibility which varies with the hydrocarbon composition (fractions of methane, ethane, propane, and other higher hydrocarbons). One of the estimates place the recoverable gas volume at ap- proximately 43 trillion cubic feet [61]. Since 2003, more than 300 horizontal, hydraulically fractured wells, have been drilled in Barnett that produce gas and brine from both sides of the synthetic and antithetic faults. Since June 2009, two water disposal wells have been operating in the footwall block of the anti- thetic fault to inject water into Ellenburger. Additionally, there are two large-volume injectors further east of the faults and one large-volume injector further west of the faults. 65 (b) (a) (c) Figure 6.2: (a) 3D view of the domain showing reservoir (green), aquifer (blue), overburden (yellow), un- derburden (red), synthetic fault (red surface) and antithetic fault (yellow surface). (b) The corner-point geomechanical mesh and (c) the corner-point flow mesh generated using the proposed meshing algorithm. 66 This complexity in spatiotemporal distribution of production and injection wells provides motivation for building a high-resolution model to study induced seismicity in the area (Figure 6.5). 6.3 Coupledmultiphaseflow–geomechanics–faultfailuremodel We provide the details of our mathematical and numerical model in appendix. We build a 3D model since the use of 2D models in induced seismicity studies limits the fault slip propagation to be 1D i.e. only along the dip or the strike direction. This may not represent real-world 3D cases, where production- and injection-induced contraction and expansion of the rock are three-dimensional, and the slip direction can evolve in time from along-strike to along-dip into overburden or underburden. A fault surface can also be non-planar or curved and can intersect another fault, as in the case of Azle’s synthetic-antithetic fault system. Such intersections activate site-specific mechanisms of pressure and stress transfer and should be honored during modeling. However, creating a 3D mesh that can honor intersection between nonplanar surfaces is non-trivial because the geometrical complexity can lead to (1) failure of the meshing algorithm, (2) large aspect ratio mesh cells causing poor mesh quality, (3) ill-conditioned finite element matrices, and/or (4) too many mesh elements causing long simulation times. Creating a geologically faithful and computationally efficient 3D mesh of a subsurface volume with multiple curved and/or intersecting faults is a challenge. Therefore, many geomechanical models either discard curved faults or simplify them as planes [38, 32, 37, 40, 63] We develop a modified corner point geometry (CPG) meshing algorithm that can incorporate faults with changing dip and strike. Details are provided in the Appendix. The dimensions of our model domain are 25 km× 15 km× 10 km in X, Y , and Z directions, respectively. It contains the two major faults mentioned above: the northeast-southwest striking Newark East Fault Zone (NEFZ) synthetic fault and a minor antithetic fault. The synthetic fault dips 81-86 ◦ to northwest and the antithetic fault dips 86 ◦ to southeast (Figure 6.2(a)). The synthetic fault surface changes its strike and bends eastward. 67 Table 6.1: Azle Model parameters h(km) k (m 2 ) ϕ S wi E (kPa) ψ ( ◦ ) ν b c (kPa) Overburden (OB) 1.8 - - - 1.44× 10 7 30 0.2 0.7 2× 10 4 Marble Falls (MF) 0.1 9.87× 10 − 18 0.20 0.5 6× 10 7 30 0.2 0.7 2× 10 4 Barnett (BT) 0.1 9.87× 10 − 21 0.06 0.5 4× 10 7 30 0.23 0.7 2× 10 4 Ellenburger (EB) 1 2.96× 10 − 14 0.055 1 6× 10 7 30 0.2 0.7 2× 10 4 Basement-A (BA) 2 9.87× 10 − 21 0.05 1 4.3× 10 7 30 0.27 0.7 2× 10 4 Basement-B (BB) 5 - - - 4.3× 10 7 30 0.27 0.7 2× 10 4 Synthetic Fault 10 k h =9.87× 10 − 19 k v =9.87× 10 − 14 0.5 0.5(MF) 0.5(BT) 1(EB) 4× 10 7 30 (OB-MF- BT-BB) 31(EB-BA) 0.25 0.7 6894.76 (OB-MF- BT-BB) 0 (EB-BA) Antithetic Fault 3 k h =9.87× 10 − 19 k v =9.87× 10 − 14 0.5 0.5(MF) 0.5(BT) 1(EB) 4× 10 7 30 (OB-MF- BT) 31(EB) 0.25 0.7 6894.76 (OB-MF- BT) 0 (EB) dp i /dz =11.28 kPa/m dσ 3,i /dz =4.5 kPa/m dσ 2,i /dz =4.95 kPa/m dσ 1,i /dz =14.7 kPa/m Petrophysical properties of the model are based on the published literature [32] and are shown in Table 6.1. Mechanical and hydraulic properties of the two faults e.g., permeability, are less certain than those of the formations. We use model calibration, where the model predicted pressures and fault failure events are compared against the measured well pressures and recorded seismicity, to estimate these pa- rameters. In order to model pressure communication between the aquifer and the basement, we increase the along-fault permeability, which models the presence of fractures located along the fault [5, 48, 1]. We also perform a sensitivity analysis, where we analyze the effect of variation in the vertical-to-horizontal fault permeability ratio on the fault pressure and failure event magnitude; see supplementary information. In order to model the permeability enhancement due to hydraulic fracturing in the gas wells, we altered the permeability field around each gas producer’s horizontal open-hole section. The permeability increases by a factor of10 7 in the well completion block and decreases logarithmically with distance from the wellbore in both horizontal and vertical directions as shown in Figure 6.3. The permeability enhancement is limited to the top four layers of Barnett and do not extend to Ellenburger or the two faults, which is expected from the average effective length of hydraulic fractures reported in the literature, which is in the range of 50-75 m approximately [47]. Based on the data available in the Azle region, we initialize pressures, saturations and stresses in the domain such that hydrostatic and mechanical equilibria are achieved at initial time, t = 0. Table 6.1 68 -20 -19 -18 -17 -16 -15 -14 10 10 10 10 10 10 10 N32.90 W97.71 2 4km Figure 6.3: Barnett permeability map shown at the Barnett Top surface. It shows the effect of hydraulic fracturing on enhancing the permeability around each gas producer’s horizontal section in the reservoir. shows the initial pressure gradient with depthdp i /dz, initial maximum horizontal stress gradientdσ 2,i /dz, initial minimum horizontal stress gradientdσ 3,i /dz, and the initial vertical stress gradientdσ 1,i /dz, where σ 1 > σ 2 > σ 3 are the principal stresses (positive under compression). The initial stress field is based on the normal faulting regime prevalent in the Azle region. Hence, the initial principal stress directions,θ 1 , θ 2 , andθ 3 , point vertically downward, eastward, and northward, respectively. The initial water saturation S wi for Barnett and Ellenburger formations and for the two faults are also shown. S wi in Barnett and Ellenburger are estimated through several iterations of model calibration (see the Model calibration section below) where we fine-tuned the S wi field to honor the observed gas and water rates and volumes shown in Figure 6.6. 6.3.1 Productionandinjectionhistory Barnett production started in year 2003 and the earthquakes occurred in 2013, 2014 and 2017. We consider the 21-year period from January 1, 2000 to December 31, 2020 to conduct our analysis. Dynamic simu- lation of the model is driven by gas and water production history from 313 horizontal and hydraulically 69 fractured wells in Barnett Shale and by water injection in five wells in Ellenburger. Figure 6.5 shows the well locations. To address the effect of model’s no-flow boundaries on well performance, we designate the 2-km wide zone along the periphery of the model as a buffer zone with no gas producers and only consider the producers located within the inner zone of the model. To model the effect of three far-away injectors, which are located at the boundaries of the model and expected to influence the Ellenburger pressure, we assigned a prorated injection rate based on their actual injection rate and the ratio of the effective injection area to the total injection area for each injector. This is based on the concept of well’s radius-of-influence in the theory of pressure transient analysis [3] and is necessary to account for the boundary effect re- sulting from the three injectors. All the well data including monthly production and injection volumes and well locations are obtained from the IHS Markit database and provided in the manuscript and the supplementary material. Figure 6.6(a), (b), (d) and (e) show the monthly and cumulative produced gas and water. The volume of the total gas produced is more than 14× 10 9 m 3 and the volume of total water produced is more than 5× 10 6 m 3 . Figure 6.6(e) also shows the cumulative volume of hydraulic fracturing fluid (primarily water) used to stimulate the producers. The noticeable correlation between the cumulative injected volume of fracturing fluid and the produced volume of total water suggests that the flowback of the injected fracturing fluid in Barnett wells is an important source of produced water. As mentioned, understanding the origin of the produced water is important to increase the accuracy of model predicted pressure and stress fields. Below we uset GP to denote the start time of gas production from Barnett. Figure 6.6(c) and (f) show the monthly and cumulative water injected in two regions: the study area and the two injectors located near the faults, which are the only injectors included in prior studies [32, 15]. The difference in the injected volume of the two regions is significant because Figure 6.6(f) shows that the total injected water volume in the study area (13.9× 10 6 m 3 ) is more than twice the volume of water injected near the faults (6.1× 10 6 m 3 ). Figure 6.6(g) shows the volume of gas and water produced and injected for 70 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 (a) (b) Figure 6.4: Observed (circles) and simulated (lines) (a) Well API-42367346930000 and (b) Well API- 42497368750000 flowing wellhead pressure calculated from the reported bottom hole pressure [32, 15]. each well relative to the volume produced and injected by the largest producer or injector. This figure shows that while some production is allocated from the hanging-wall block, the majority of production is from the footwall block. Figure 6.6(g) also shows that both sides of the synthetic fault experience water injection. We uset WI to denote the start time of injection near the faults. Figure 6.4 shows the flowing well head pressure for the two injectors located near the faults. The flowing well head pressures are calculated from the reported bottom hole pressures for these wells [32, 15]. A permeability multiplier of 0.33 was imposed on the injector well permeability to capture the effect of scale build-up (e.g., salt precipitation) around the injectors [3]. 6.4 Modelcalibration In order to ensure model fidelity and confidence in simulation results, we performed a multiscale cali- bration study. We performed calibration at the well scale, reservoir layer scale, and the entire field scale, which in turn reduced uncertainty in the near wellbore properties, layer-wise properties and field-scale 71 2km 2 4km N32.90 W97.71 Figure 6.5: Map showing the location of gas producers (highlighted in red) and water injectors (highlighted in blue). The area outside the dashed green rectangle represents a buffer zone where no gas producers are considered in order to minimize the boundary effect. Below, we call the two injectors inside the dashed green box as near-fault injectors and the other three injectors as far-field injectors. poroelastic properties while honoring the production and injection volumes and the injection wellhead pressures. At the well scale, we used the well trajectory and completion information to drive calibration. We started from simple vertical wells with single cell open-hole completion to a more complex representation of horizontal wells with an open-hole completion interval across five cells in the mesh. We found that the latter representation provides a better agreement between predicted and observed gas and water rates and volumes. Moving up to the Barnett/Ellenburger layer scale, we used the uncertainty in the layer permeability field, which is altered during hydraulic fracturing or long-term injection, to drive calibration at that scale. The original permeability is altered in a logarithmic-scale fashion mimicking the hydraulic fracturing for each producer as shown in Figure 6.5. Next, the permeability around the two water injectors is reduced to replicate the effect of scale build-up around each injector. Next, we move up to the field scale where quantities such as the initial water saturation (S wi ) in Barnett are uncertain. We started with a 72 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 (a) (c) (b) (d) (f) (e) 0 2 4 6 8 10 12 14 16 10 4 0 1 2 3 4 5 6 7 8 9 10 10 4 0 2 4 6 8 10 12 14 10 7 0 2 4 6 8 10 12 14 16 10 9 0 1 2 3 4 5 6 7 10 6 0 2 4 6 8 10 12 14 16 10 6 (g) N32.90 W97.71 2 4km Figure 6.6: Observed (circles) and simulated (lines) (a) gas monthly production rateq g , (b) water monthly production rate q wp , (c) water monthly injection rate q wi , (d) cumulative gas production volume Q g , (e) cumulative water production volume Q wp , and (c) cumulative water injection volume Q wi in the Azle field plotted with time. All volumes are at surface pressure and temperature conditions. (g) Bubble map of cumulative produced and injected volumes per well. The bubble size is linearly proportional to produced or injected fluid’s volume. Produced water circles have been magnified 3.3 times for the ease of visualization. 73 low initial water saturation of 0.25 and compared the predicted production rates and predicted cumulative volumes to the observed rates and volumes. We increased the initial water saturation gradually until we achieved a better agreement between the predicted and observed quantities. While calibratingS wi , we investigated three possible sources of the produced water from Barnett wells: (1) water seeping from the Ellenburger formation through natural fractures that may exist and hydrauli- cally connect Barnett Shale and Ellenburger dolomite, (2) hydraulic fracturing flow-back water that was injected during fracturing jobs and got produced after the fractured wells were put on production, and (3) formation water of Barnett. We increased the vertical permeability in the bottom layer of Barnett Shale to allow hydraulic communication between Ellenburger and Barnett and found that this results in a poor agreement between observed and predicted fault failure event times because it allows injection-induced overpressure to diffuse upward into the shallower Barnett formation. This finding, along with the strong correlation between the injected fracturing fluid volume and the produced water volume in Figure 6.6, supports the hypothesis that the flowback water, not water from Ellenburger, is an important contributor to the produced water from Barnett wells. This reduced uncertainty onS wi in the model. Finally, we used the time and location data of observed seismic events to calibrate the faults vertical permeability and co- hesion parameter. Figure 6.7 shows the effect of variation in the fault vertical-to-horizontal permeability ratio on the timing and magnitude of fault failure events and the fault pressure. We used this sensitivity analysis to use a permeability ratio value in the model that can yield pressure, stress and failure results in alignment with the observed data. 74 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 0 1 2 3 4 5 6 7 8 9 10 -6 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 7420 7440 7460 7480 7500 7520 7540 7560 7580 (a) (b) Figure 6.7: Evolution of the fault failure magnitude (plastic strain) and the fault pressure in a cell of the synthetic fault for six different values of the vertical-to-horizontal permeability ratio. The cell is located in the basement (I=22, J=12, K=10). 6.5 Results 6.5.1 Pressureandstressdistributions Figure 6.8 shows reservoir overpressure∆ p, with respect to the initial pressure, at six selected timesteps starting approximately two years aftert GP until the end of year 2020. Hydraulic fracturing-induced perme- ability enhancement around each gas producer affects depletion and pore pressure distribution in Barnett. Even though the gas producers are distributed in both fault blocks of the synthetic fault, the footwall block exhibits a higher pressure depletion than the hanging-wall block, which agrees with the observed produc- tion volumes. Ellenburger injection has only a minor effect on Barnett’s pressure because Barnett’s low permeability slows down the propagation of pressure diffusion front. Figure 6.9 shows that the change in Barnett’s effective mean stress, ∆ I ′ 1 , which measures the production-induced volumetric contraction of Barnett, evolves in sync with Barnett’s∆ p shown in Figure 6.8, i.e. more depletion corresponds to more compressive stresses and higher∆ I ′ 1 values. 75 (a) 01-May-2005 (b) 01-May-2009 (c) 01-Jan-2011 (d) 01-Dec-2013 (e) 01-Aug-2017 (f) 31-Dec-2020 X Y -2000 -1500 -1000 -500 0 2 4km N32.90 W97.71 Figure 6.8: Barnett overpressure ∆ p at the selected time steps: a, b, c, d, e, and f. Fracturing-induced permeability enhancement around gas producers lead to heterogeneity in the pressure field. Due to smaller injection volumes, the footwall block experiences larger depletion compared to the hanging wall block. 76 (a) 01-May-2005 (b) 01-May-2009 (c) 01-Jan-2011 (d) 01-Dec-2013 (e) 01-Aug-2017 (f) 31-Dec-2020 -100 0 100 200 300 400 500 X Y 2 4km N32.90 W97.71 Figure 6.9: Change in the effective mean stress ∆ I ′ 1 (with respect to the initial value) at the selected timesteps: a, b, c, d, e, and f in Barnett. ∆ I ′ 1 , which is positive in compression, shows a monotonic re- sponse to∆ p in Figure 6.8. 77 Figure 6.10 shows Ellenburger’s∆ p at the same time steps. Ellenburger’s∆ p evolves more uniformly compared to Barnett due to a more homogeneous permeability field. While low horizontal permeabil- ity k h impedes pressure communication between the two fault blocks, both sides of the synthetic fault exhibit pressure rise due to injection on both sides of the fault. The aquifer area between the synthetic and antithetic faults communicates with the hanging-wall block and experiences a slight overpressure. Decrease in ∆ I ′ 1 in Figure 6.11 indicates tensile or expansion response in regions with high ∆ p. An in- teresting observation is that production-induced compressive stresses in Barnett have propagated deeper into Ellenburger and can counterbalance the injection-induced tensile stresses in Ellenburger. This creates distinct regions of induced stress at a length scale that is larger than well’s radius-of-influence and smaller than the fault block size. Boundaries between these regions play an important role in rotating the principal stress directionsθ 1 andθ 3 , which is discussed below. 6.5.2 Interactionandpropagationofstressesalongdepthandwithtime Figure 6.12 shows two vertical cross-sections along an XZ plane in the model that combines ∆ p in (a) and (d) and the vertical displacementu z in (b) and (e) at the end of Year 2020. In (c) and (f), we see depth profiles of the quantities along a highlighted column of cells in the vertical direction. In the map view, the profiles are taken at locations that are experiencing highest amount of pressure depletion in Barnett. To understand time evolution, the profiles are shown at four selected time steps. We observe that Barnett’s pressure depletion in the footwall block is larger than that in the hanging wall block, while Ellenburger’s pressure rise in the two fault blocks are relatively similar in magnitude. The change in pressure drives the change in the state of stress. This is observed in the change of the effective mean stress ∆ I ′ 1 and the von Mises stress ∆ √ J 2 . The latter measures the shear or distortion response in the rock; see the discussion in 3 Eq. (3.5). This is also observed in theu z plot, where we see that the footwall block experiences larger subsidence than hanging-wall block, which is aligned with the location of the wells. 78 0 200 400 600 800 1000 1200 (a) 01-May-2005 (b) 01-May-2009 (c) 01-Jan-2011 (d) 01-Dec-2013 (e) 01-Aug-2017 (f) 31-Dec-2020 X Y 2 4km N32.90 W97.71 Figure 6.10: Ellenburger overpressure∆ p fields at different time steps: a, b, c, d, e, and f. The uniform nature of the aquifer permeability causes ∆ p in both fault blocks to evolve monotonically with injection. Low values of the fault horizontal permeability k h , relative to the Ellenburger permeability, impede pressure communication between the two fault blocks, however both blocks exhibit relatively balanced overpressure due to water injection. 79 -350 -300 -250 -200 -150 -100 -50 0 (a) 01-May-2005 (b) 01-May-2009 (c) 01-Jan-2011 (d) 01-Dec-2013 (e) 01-Aug-2017 (f) 31-Dec-2020 X Y 2 4km N32.90 W97.71 Figure 6.11: Changes in Ellenburger’s effective mean stress ∆ I ′ 1 at different time steps: a, b, c, d, e, f. Changes are with respect to the initial value. While the changes mostly follow the behavior of Ellen- burger pressure, there are areas with localized higher-order effects in ∆ I ′ 1 that do not directly follow Ellenburger’s pressure and indicates stress propagation from Barnett. These are areas with small mag- nitudes of∆ I ′ 1 (closer to zero) where injection-induced tensile stresses are counterbalanced by Barnett’s production-induced compressive stresses. 80 500 1000 1500 2000 2500 3000 3500 4000 4500 1000 500 0 -500 -1000 -1500 -2000 -0.03 -0.01 0.01 0.03 0.05 0.07 -2000 0 2000 -500 0 500 -500 0 500 -0.02 0 0.02 -2000 0 2000 -500 0 500 -500 0 500 -0.02 0 0.02 500 1000 1500 2000 2500 3000 3500 4000 4500 (a) (b) (d) (e) (f) 0 m 1800 m 3000 m 10000 m May 1, 2009 January 1, 2011 December 1, 2017 December 31, 2020 (c) Synthetic Fault Hanging-Wall Block Synthetic Fault Foot-Wall Block Figure 6.12: (a) Overpressure∆ p in the flow domain along a XZ vertical section. (b) Vertical displacement u z in the entire geomechanical domain along the same vertical section. (c) Depth profiles of ∆ p, changes in the effective mean stress ∆ I ′ 1 and von Mises stress∆ √ J 2 , andu z along a column of cells in the hanging- wall block. (d), (e), and (f) are similar plots in the footwall block. 81 6.5.3 Stressrotation Figure 6.13(a-1) and (a-2) show overpressure∆ p in Barnett and Ellenburger formations on December 31, 2020. Figure 6.13(a-2) shows the change in maximum and minimum effective principal stress angles ( ∆ θ 1 and ∆ θ 3 ) along the trace A-A’ of the synthetic fault in Barnett formation. Figure 6.13(b-2) shows this in Ellenburger. Figure 6.13(a-3) shows the changes in stress orientations along an east-west line B-B’ in Bar- nett. Figure 6.13(b-3) shows this in Ellenburger. We observe that∆ θ 3 is significantly larger in magnitude than∆ θ 1 . This agrees with our expectation because the maximum principal stressσ ′ 1 in the Azle region acts vertically downwards (normal faulting regime). Also, we find that the principal directions rotate more in Barnett than in Ellenburger. This can be attributed to low permeability areas, which surround hydrauli- cally stimulated high permeability areas around Barnett wells. Pressure diffusion is impeded across the boundary of Stimulated Rock Volume (SRV), which creates a contrast in principal stresses and directions between inside and outside the SRV. Ellenburger, on the other hand, has a more uniform and higher per- meability and pressure values and hence a smaller rotation of the principal directions. Next observation is regarding the effect of fault (position and geometry) on stress rotation. Figure 6.13(a-2) and (b-2) show that stresses reorient along the synthetic fault with a larger rotation near the center of the field closer to producers and injectors and fault’s bending point. Figure 6.13(a-3) and (b-3) show that the rotation along the east-west direction is concentrated near the synthetic fault. This agrees with the fact that the fault hinders pressure communication between the two fault blocks. Figure 6.14(a) and (b) show the Dec 2020 orientations of the minimum effective principal stress σ ′ 3 and the intermediate effective principal stress σ ′ 2 for all cells that exhibited a rotation greater than 1 degree in Barnett and Ellenburger, respectively. There are two regions where the rotation of principal directions is severe. The first region is along the two faults, especially along the synthetic fault. This is due to the fault acting as a barrier against pressure commu- nication between the two fault blocks. The second region is around the Barnett producers and between Ellenburger injectors. While production from Barnett causes pressure depletion and compressive stresses 82 -2000 -1500 -1000 -500 0 1000 1100 1200 9 8 7 6 5 4 3 2 1 1 9 8 7 6 5 4 3 2 1 1 9 8 7 6 5 4 3 2 1 1 9 8 7 6 5 4 3 2 1 1 (a) (d) (b) (e) (c) (f) 01-May-2009 01-Jan-2011 01-Dec-2013 31-Dec-2020 2 4km 2 4km N32.90 W97.71 N32.90 W97.71 Figure 6.13: Change in the maximum and minimum effective principal stress angles ( ∆ θ 1 and∆ θ 3 ) from their initial values in Barnett and Ellenburger. (a) shows the Barnett overpressure in Dec 2020. (b) shows ∆ θ 1 and∆ θ 3 along section A-A’ of the synthetic fault at Barnett depth at several time steps. (c) shows∆ θ 1 and∆ θ 3 along an east-west section B-B’. (d) shows the Ellenburger overpressure in Dec 2020. (e) shows ∆ θ 1 and∆ θ 3 along section C-C’ of the synthetic fault at Ellenburger depth at several time steps. (f) shows ∆ θ 1 and∆ θ 3 along an east-west section D-D’. 83 a) b) -2000 -1500 -1000 -500 0 1000 1100 1200 2 4km 2 4km N32.90 W97.71 N32.90 W97.71 Figure 6.14: Rotation of the minimum effective principal stress σ ′ 3 and the intermediate effective principal stress σ ′ 2 on December 2020 in (a) Barnett and (b) Ellenburger. Well rates have been projected to future based on June 2019 rates. Since hydraulic fractures tend to orient perpendicular toσ ′ 3 , rotation in principal stresses have implications for future drilling and fracturing activities in the field. within the hydraulically fractured area, this effect is localized and does not extend to the un-fractured rock due to the ultra-low permeability of Barnett rock matrix. This causes pressure depletion to localize around a cluster of closely spaced producers and causes tensile stresses in the surrounding area away from the cluster. Lower depletion also translates into a relatively higher pore pressure in the un-fractured, low- permeability area. This phenomenon of induced anisotropy and heterogeneity in the stress state is most clearly observed in ultra-low permeability reservoirs undergoing hydraulic fracturing as it creates a clear contrast in permeability between fractured and un-fractured rocks. This information is needed to opti- mize the trajectories of new and infill horizontal wells and their hydraulic fracturing jobs. Consequently, the figure shows the importance of modeling spatiotemporal changes in principal stress magnitudes and directions. 6.5.4 Faultstressesandinducedfailuremechanisms Figure 6.15 shows a 3D view of the model interior: contours of overpressure ∆ p in Ellenburger on De- cember 31 2020, fault failure locations identified by ε p plotted on synthetic and antithetic faults, and the 84 Ellenburger Basement Overburden (a) (b) 0 0.2 0.4 0.6 0.8 1 1.2 10 -5 X 0 3 10 2 1.9 1.8 X Y Figure 6.15: (a) Failure magnitude ε p plotted on the synthetic and antithetic faults and overpressure ∆ p plotted as contours on a horizontal cross-section of Ellenburger. Bottom-hole locations of the water injec- tors are shown with triangles. (b) shows a rotated view of (a) without the Ellenburger section. Locations and magnitudes of the failure events on the faults show that deeper section of the synthetic fault experi- ences larger magnitude events compared to events located at or above the injection depth. locations of the five injectors (the 5th injector near the northeast corner of Figure 6.5 is hidden). The verti- cal location of failure events on the synthetic fault and their magnitudes in Figure 6.15(b) show that deeper sections of the fault host larger magnitude events compared to events at or above the well depths in Bar- nett and Ellenburger. This supports our hypothesis of induced stresses propagating vertically downward from Barnett into the basement, where a velocity-weakening frictional rheology of the fault is expected to support earthquake nucleation [56]. 6.5.4.1 SyntheticFault In order to understand how the stresses propagate to the basement section of the faults, Figure 6.16 shows the spatiotemporal evolution of overpressure ∆ p, change in von Mises stress ∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure strain dε p for two columns of cells located on the synthetic fault. In this and subsequent plots, we use the mesh axes (I, J, K), which are quasi-parallel to Cartesian axes (X,Y,Z) (see below appendix on meshing), to indicate the location of mesh cells being analyzed in a plot. For example, the two columns of cells in Figure 6.16 are at J=14 and 15, which we select because 85 they are locations along the strike of the synthetic fault between the two near-fault injectors. From the line plots, we see that∆ p response to injection is similar across depth for all the cells. Note that pressure starts increasing beforet WI , marked with an inverted triangle, because the far-field injectors start injecting before the near-fault injectors. However, the failure onset time, i.e. the time whendε p becomes positive for the first time, is different across depth: the failure onset times for K=6 through 10 are later than that of K=5. Also, failure at K=5 occurs beforet WI , i.e. before injection begins in the two near-fault injectors. We hypothesize that this behavior results from propagation of production-induced stresses from Barnett into Ellenburger. We test the hypothesis in Figure 6.17, which shows the stress paths of the same cells at J=14. The paths are shown in the stress invariant space of I ′ 1 vs. √ J 2 , and the direction of time is implied by the curved arrow shown in the plot. The stress path of the J=14, K=5 cell is unique from other cells, that are in deeper layers, as shown by the spacing between thet GP andt WI markers corresponding to the start of gas production and water injection, respectively. Production-induced stresses dominate stress evolution at K=5 causing a longer spacing betweent GP andt WI and larger deviatoric stresses ( √ J 2 ) compared to the other cells. The onset of failure is marked by the asterisk symbol, which is prior to injection in case of K=5 layer. To account for the uncertainty in the initial state of stress (magnitude and direction), we performed two additional simulations, Run-2 and Run-3. In Run-2, we increased initial stresses by 30% and in Run-3 we rotated the initial stress directions such that the maximum horizontal principal stress is along NE32 ◦ . Figure 6.18 (a), (b), (c) and (d) shows the evolution of pressure, stress and failure quantities (for comparison to base case Run-1, Figure 6.16) for cells located on the synthetic fault in layers K=5 and 6. From Figure 6.18 (e), we can observe that the fault failure time of Run-2 is earlier than that of Run-1, which is expected from the higher initial stresses in that run. 86 2000 2005 2010 2015 2020 0 500 1000 0 1 2 3 10 -7 J=15 K=5 J=15 K=6 J=15 K=7 J=15 K=8 J=15 K=9 J=15 K=10 J=14 K=5 J=14 K=6 J=14 K=7 J=14 K=8 J=14 K=9 J=14 K=10 Figure 6.16: Evolution of overpressure ∆ p, change in von Mises stress ∆ √ J 2 , change in effective mean stress∆ I ′ 1 , and instantaneous failure straindε p in synthetic fault cells located at the depths of Ellenburger layers. K=5 through 10 are top six layers in Ellenburger. ∆ p,∆ √ J 2 , and∆ I ′ 1 are plotted against the left axis whiledε p is plotted against the right axis. Axis labels, shown only in the bottom left plot, are identical for all other plots. 87 K=5 K=8 K=6 K=7 K=9 K=10 1.78 1.785 1.79 1.795 1.8 1.805 1.81 1.815 10 4 1.28 1.282 1.284 1.286 1.288 1.29 1.292 1.294 1.296 1.298 10 4 1.3 time Jan 1, 2000 Figure 6.17: Stress trajectories of synthetic fault cells in theI ′ 1 − √ J 2 stress invariant space over the 2000- 2020 period. K=5 through 10 are the top six layers in Ellenburger. The K=5 cell is chosen as the reference cell from which other cells are compared. Hence, all cells start at the same initial point, which corresponds to the initial stress state of K=5. The start time of gas production from Barnett (t GP ), the start time of water injection in Ellenburger (t WI ), and the failure onset time (t p i ) are marked. Failure in K=5, which is in Ellenburger, begins before injection begins in Ellenburger. This suggests propagation of production- induced stresses from Barnett into Ellenburger. 88 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 -200 -150 -100 -50 0 50 100 150 200 -4 -3 -2 -1 0 1 2 3 4 10 -6 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 -200 -150 -100 -50 0 50 100 150 200 -4 -3 -2 -1 0 1 2 3 4 10 -6 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 -200 -150 -100 -50 0 50 100 150 200 -4 -3 -2 -1 0 1 2 3 4 10 -6 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 2022 -200 -150 -100 -50 0 50 100 150 200 -8E-08 -6E-08 -4E-08 -2E-08 0 (a) (22,15,5) (c) (22,15,5) (b) (22,15,6) (d) (22,15,6) 2E-08 4E-08 6E-08 8E-08 1E-07 (e) Figure 6.18: Evolution of overpressure ∆ p, change in von Mises stress ∆ √ J 2 , change in effective mean stress∆ I ′ 1 , and failure strainε p in synthetic fault cells located at the depths of Ellenburger layers for the two simulation runs with different initial stress values (Run-2) and orientations (Run-3). (e) compares the evolution ofε p in synthetic fault cell (22,15,5) for Run-1 (Base Case), Run-2, and Run-3. 89 Figure 6.19: Top view of model showing the change in effective mean stress distribution at last time step for Barnett and Ellenburger at two different additional simulation runs with initial stress values and ori- entations. Figure 6.20: Top view of model showing the vertical displacement distribution at last time step at two different additional simulation runs with initial stress values and orientations. 90 2005 2010 2015 2020 0 500 1000 0 1 2 3 10 -7 2000 Figure 6.21: Spatiotemporal evolution of overpressure ∆ p, change in the effective mean stress ∆ I ′ 1 and change in the von Mises stress∆ √ J 2 , and instantaneous failure straindε p for four cells in the antithetic fault within Ellenburger interval. Top row: locations of the cells in 3D. Middle and bottom rows: ∆ p, ∆ √ J 2 ,∆ I ′ 1 are plotted against the left axis whiledε p is plotted against the right axis. Middle row, which corresponds to layer K=5, shows failure events on thedε p curve after water injection begins in the near- fault injectors. The injection start time, t WI , is marked. Bottom row shows that the K=6 layer does not experience failure within the simulated period. 6.5.4.2 AntitheticFault Next, we analyze the evolution of state variables on the antithetic fault, which extends from surface to the bottom of the aquifer. Figure 6.21 shows the evolution of overpressure, changes in the effective mean stress and the von Mises stress, and the failure event magnitudes. The overpressure evolution does not show any effect of far-field injectors on the antithetic fault because the synthetic fault, which has a low across-fault permeability, impedes pressure communication between the antithetic fault and the far-field injectors. This is in contrast to the overpressure evolution on the synthetic fault shown in Figure 6.16. Figure 6.22 shows a close-up view of the stress paths of two cells in the antithetic fault in Ellenburger layers K=5 and 6. As pressure depletes in Barnett and production-induced stresses propagate into Ellenburger, √ J 2 increases 91 1.28 1.284 1.288 1.292 1.296 10 4 1.785 1.79 1.795 1.8 1.805 1.81 1.815 x10 4 Figure 6.22: Stress paths of antithetic fault cells in theI ′ 1 − √ J 2 space. Two cells from Ellenburger’s top layer (K=5) and second layer (K=6) are plotted. The top layer cell shows a more complex stress path (non- monotonicity inI ′ 1 ) due to its proximity to Barnett and hence to production-induced stresses. in Ellenburger layer K=5, but not in layer K=6 because the former is closer to Barnett. The curved stress path of layer K=5 also shows that production-induced compression causes positive changes in I ′ 1 while injection-induced expansion causes negative changes inI ′ 1 . 6.5.4.3 Evolutionofprincipalstressesduetointeractionbetweeninducedstresses In Figure 6.23, we analyze the time evolution of the three effective principal stresses, σ ′ 1 ,σ ′ 2 ,σ ′ 3 , in two cells in the synthetic and antithetic faults. We can see that the cell K=5 on the antithetic fault, which lies just below Barnett, exhibits a larger change in stress during the period 2006-2009, which is prior to near fault injection, compared to the deeper cell K=6 closer to the basement. This effect is caused by a gas producer that is located in the vicinity of cell K=5. This gas producer came online in mid 2006 and caused the pressure in cell K=5 to deplete which in turn caused∆ σ ′ 1 ,∆ σ ′ 2 and∆ σ ′ 3 to increase. We can see that this effect begins to diminish after 2008 when injection starts from the far-away injector on the western boundary, causing cell K=5 pressure to increase. This event is followed by a near-fault injector coming 92 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 -350 -300 -250 -200 -150 -100 -50 0 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 -250 -200 -150 -100 -50 0 J=10 K=5 J=14 K=6 J=10 K=6 J=14 K=10 (b) (a) Figure 6.23: Temporal evolution of the change in effective principal stresses, ∆ σ ′ 1 ,∆ σ ′ 2 and∆ σ ′ 3 , in selected cells in (a) the antithetic fault and (b) the synthetic fault. Cells that are vertically closer to Barnett exhibits stress evolution after the start of production, while cells that are located deeper do not exhibit the same characteristics. As injection starts near the faults, stress acting on all four cells decrease due to the increase in pressure. The separation between principal stresses increase due to interaction between production- and injection-induced stresses. 93 online in mid 2009 (indicated by t WI marker in the plot), causing further pressure increase in cell K=5 which in its turn causes∆ σ ′ 1 ,∆ σ ′ 2 and∆ σ ′ 3 to decrease or become less compressive. We know from the Mohr Circle analysis [35] that as the principal stresses decrease they get closer to the Coulomb failure line. Also, as the difference between the principal stresses increase, which is obvious from the plots after t WI , the circle becomes larger and more likely to touch the failure line. Both mechanisms induce shear failure. Our main observation is that, in case of Azle, production-induced stress, injection-induced stress, and fault pressurization interacted with each other to drive failure in time, and this interaction varies with depth. Figure 6.24 shows the location and magnitudes of the failure events cumulatively at different times on both synthetic and antithetic faults. While only the first layer of the aquifer exhibits failure events due to production, the majority of events occurred in the basement by the end of 2013, which is aligned with the recorded seismicity data. 94 (a) 01-May-2009 (b) 01-Jan-2011 (c) 01-Dec-2013 (d) 01-Aug-2017 Synthetic fault events Antithetic fault Events 0 3 10 2 2.33 0.0641 1.8305 0.3664 2000 2002 2004 2006 2008 2010 2012 2014 2016 2018 2020 0 1 2 3 4 Observed Synthetic fault events Antithetic fault events 2013 2014 2015 2016 2017 2018 0 1 2 3 4 5 6 7 8 9 10 2.8 3.1 2.2 3.3 3.3 3.1 2.8 2.7 3.4 2.1-3.6 3.6 (e) (f) Figure 6.24: 3D views of the synthetic and antithetic faults showing the induced failure events at the same four timesteps analyzed above: a, b, c, and d. All events occurring at or before a timestep are shown at that timestep. The size of the circle or square is a proxy for the magnitude of the calculated event. The events are more concentrated and are larger in magnitude in the basement interval. Green shaded cells represent fault intersection with Barnett, and blue shaded cells represent fault intersection with Ellenburger. (e) Temporal distribution of simulated and observed events on synthetic and antithetic faults. The Nov 2013– Jan 2014 seismic cluster and a recent Mb 2.8 event in June 2017 are shown with asterisks. (f) Depth and magnitudes of the observed seismic events. 95 Chapter7 DFW-IrvingIncucedSeismicEvents 7.1 Introduction The Irving-Dallas earthquake sequence in North East Texas during 2008-2015 has been hypothesized as induced seismicity [33, 42]; faults destabilized by wastewater injection into the Ellenburger dolomite for- mation, which lies below the gas-producing Barnett Shale formation. The first cluster of seismic events occurred in 2008 near DFW and were located on a set of normal faults crossing through Barnett and El- lenburger and extending into the basement [21]. Gas and brine producing Barnett wells operate adjacent to the faults in addition to two water disposal wells in Ellenburger that operate in close vicinity of the hypocenters of the events. The second cluster of seismic events occurred in 2015 near Irving, Texas and were located on a normal fault crossing through Barnett and Ellenburger and extending into the basement. These events share similar characteristics with other injection-induced seismic events reported in the lit- erature, e.g. the 2013-2014 Azle-Reno events. As a result, it was hypothesized that Ellenburger injection and the confinement created by the low permeability Ouachita fold and thrust belt on the east led to an increase in the basin pressure, which in turn led to the onsetting of seismic events [33]. However, the activation mechanism of Irving-Dallas events is not well understood as it is shrouded in ambiguity due to many event locations being far away from any extraction–disposal wells [18]. The nearest injector is located 11 miles from the Irving event epicenters. This is further complicated by gas and water production 96 from horizontal, hydraulically fractured Barnett Shale wells and how the resulting production-induced stress interacted with the Ellenburger injection-induced stress. The presence of extensive faulting [31], stratigraphic offset across faults, and a dipping stratigraphy (south-eastward dip) make data interpretation and causative analysis complex. There is no mechanistic model of Irving-Dallas seismicity that accounts for the production and injection-induced stress changes and their effect on fault stability while honoring the complex regional geology. To test the injection-induced seismicity hypothesis, a call for detailed mod- eling studies was made [33]. Moreover, for the model to correctly identify the physical mechanisms and test the induced seismicity hypothesis, it must be based on the coupled physics of fluid flow (gas and wa- ter phases), poroelasticity, and fault failure. Here, we address this challenge by developing a mechanistic modeling approach to study Irving-Dallas sequence and use the approach to test new hypotheses. Understanding and predicting the pressure and stress propagation from injection and production lo- cations onto nearby basement faults and their propagation onto far away faults is essential to forecast and control the induced seismicity hazard at such extraction–disposal sites. This requires understanding how production- and injection-induced stresses interact to destabilize a fault intersecting or bounding the production and injection regions. In this regard, some of the prior studies [50] motivated our work to ex- pand the poroelastic modeling to include DFW and Irving events and to see how they are correlated with increasing pore pressure. We hypothesize that the key to discovering the causative mechanisms behind DFW and Irving seismic sequences is to understand the spatiotemporal evolution of pressure and stress fields at both the basin scale and the fault block scale. This requires a modeling approach that can quantify spatiotemporal propagation of production- and injection-induced stresses from wells to the faults while resolving fault geometry and stratigraphy, and well activity. However, constructing one such detailed model for the entire basin is computationally prohibitive due to millions of grid cells needed to discretize the basin at that resolution. Based on our analysis of the data on well activity and fault position, we propose a novel two-model approach that exploits the disparity in scales 97 between the basin-scale injection analysis and the well-scale fault reactivation analysis. We demonstrate the effectiveness of the approach by constructing a coarse-scale model of Ellenburger injection in the basin, which provides time-dependent pressure boundary conditions to a fine-scale model of the DFW- Irving region containing the faults, which hosted the seismicity, and the production/injection wells in that region. Analysis of the model simulation results provides evidence for interaction between Barnett’s production- and Ellenburger’s injection-induced stresses as well as pressure diffusion from Ellenburger into basement along the through-going faults. It allows us to successfully test the hypothesis of injection- induced reactivation as the causative mechanism for the Irving events. 7.2 Geologicstructure,fluidvolumesandwellconfiguration The location of the study area is in the north of Texas, in the vicinity of DFW Airport and Irving city (Figure 7.1a). Primary features of the geological structure of the selected area are a 100-150 m-thick Bar- nett shale layer overlying 700-1000 m-thick Ellenburger Group dolomite formation, which overlies a thick granite basement. Barnett and Ellenburger formations are present throughout the Fort Worth basin (ap- proximately 55000 square km). See the stratigraphic section in Figure 7.1b. The region is dominated by normal faulting [31] with the maximum principal stress being vertical [52]. The DFW-Irving earthquake sequences were hosted by multiple such normal faults. The study area encloses eight normal faults that generally strike in the northeast-southwest direction and extends in depth through Barnett and Ellen- burger. Figure 7.1b shows seven faults (in red color), out of the eight, that appear in this cross-section. It also shows a potential fault, in yellow color, that dips antithetic to the existing faults. This fault, which is unmapped in the literature, can be understood as the nodal plane passing through the locations of seismic events in this depth interval as shown in Figure 2a of [31]. The original gas-in-place (OGIP) in Barnett has been estimated to vary between 120 and 200 billion standard cubic feet per square mile [20] depending on the reservoir thickness, organic content, porosity, 98 and gas compressibility which varies with the hydrocarbon composition (fractions of methane, ethane, propane, and other heavier hydrocarbons). One of the estimates places the recoverable gas volume at approximately 43 trillion standard cubic feet [60]. More than 12000 wells, mostly horizontal hydraulically fractured, have been drilled into Barnett shale during 2000-2010. Since 2007, more than 170 hydraulically fractured wells have been drilled in Barnett that produce gas and brine from the vicinity and west of the DFW Airport. The injection of disposed saltwater into Ellenburger has been going on since 2006. In addition to the two water disposal wells in the vicinity of DFW Airport, more injection wells are located to the west and south of our study area. This complexity in the spatiotemporal distribution of production and injection wells augments the complexity imposed by the geologic structure, in particular faulting, on estimating stress changes and fault destabilization. A high-resolution multiphysics model is required to incorporate such complexities to study induced seismicity in the area. 7.3 Coupledmultiphaseflow–geomechanics–faultfailuremodel Production- and injection-induced contraction and expansion of the rock are three-dimensional in geom- etry. Similarly, fault failure can propagate along both strike and dip directions as stresses propagate from reservoir intervals into overburden or basement [36, 63]. A fault surface can be non-planar (changing dip and strike) or intersect/connect to another fault, e.g., in the case of synthetic-antithetic fault pairs. An advanced 3D meshing algorithm is required to honor such geometric complexities present in our Irving- Dallas model. We developed and improved our previous meshing algorithm [11] to address these chal- lenges. Using this algorithm, we build two models: (1) A coarse-scale model named the Bend Arch–Fort Worth Basin model and (2) the fine-scale model named the DFW-Irving model. 99 Wise Denton Tarrant Ellis Parker Hood Johnson Somervell Dallas Fort Worth Ouachita Thrust Front Muenster Arch 0 40 30 20 10 km Injection Well Seismic Event Fault N A B ELLENBURGER BARNETT BASEMENT CAPROCK Gas & water Production Water injection Hyd. Frac. Unmapped Antithetic Fault? DFW Events Irving Events Fracs into basement A B (a) (b) 0 48.73 0 8 1 4 3 5 6 8 9 7 km Figure 7.1: (a) A map view of the DFW-Irving study area in North East Taxes showing the 2008-2015 seismic events and the major fault lines. County lines and county names are shown in grey color. Two major faults bound the study area from the north and the east, namely Muenster Arch fault and Ouachita Thrust Front. (b) Vertical stratigraphic section (schematic not to scale) of the area along a line AB oriented west-north- west to east-south-east. It shows gas and water production from hydraulically fractured Barnett wells and water injection into Ellenburger wells. Eight out of a total of nine faults are shown. Fault-2 does not become visible in this cross-section. The stars indicate the general hypocenters locations. Volumetric contraction of Barnett is shown by inward pointing arrows on a red circle. Volumetric expansion of Ellenburger is shown by outward pointing arrows. Fluid and pressure diffusion into faults from Ellenburger to basement is shown by wiggly arrows. 100 7.3.1 BendArch–FortWorthBasinmodel (Figure 7.2a) are 150 km× 150 km× 1 km in X, Y , and Z directions, respectively. It contains three major faults: Muenster Arch fault with a northwest-southeast strike, Ouachita Thrust Front fault with a northeast-southwest strike, and the DFW Airport fault with a quasi-parallel northeast-southwest strike. The DFW Airport fault is the longest fault west of Ouachita crossing Tarrant, Johnson and Somervell counties in Figure 7.2a. 7.3.2 DFW–Irvingmodel The DFW–Irving model domain dimensions are 49 km× 38 km× 8 km inX,Y , andZ directions, respec- tively. It contains eight normal faults with a general northeast-southwest strike in addition to a potential antithetic fault (Figure 7.2c). Eight major synthetic faults that dips west to east in addition to one minor potential antithetic fault that dips in the opposite direction are studied in this model (Figure 7.2). All fault lines curve along the North-South direction, and some faults curve more than others. Next, we discuss the distribution of petrophysical properties in the model. To address the limitation in prior studies associated with assuming a homogeneous permeability, we model Barnett and Ellenburger log-permeability fields as multigaussian random fields characterized by an average permeability and standard deviation of the per- meability. We generate these fields using SGeMS, a geostatistical tool [53]. Since porosity and permeability are often related to each other via rock’s porous texture, we generate multigaussian fields for distributing porosity as well. The average permeability in Barnett is9.87× 10 − 21 m 2 and the standard deviation of per- meability is6.76× 10 − 21 m 2 . The average permeability in Ellenburger is2.96× 10 − 14 m 2 and the standard deviation of permeability is4.74× 10 − 15 m 2 (Figure 7.4). Some of the values are based on the published literature [32]. We assume transverse isotropy of the depositional environment, i.e.k x =k y ̸=k z with the anisotropy ratiok z /k x =0.1, which is a commonly assumed value in subsurface flow simulation models. Overburden, Marble Falls and Basement layers are not part of the flow model. Mechanical and hydraulic 101 properties of the faults are less certain than those of the formation layers. We use model calibration, where the model predicted pressures and fault failure events are compared against the measured well pressures and recorded seismicity, to estimate such parameters. In order to model pressure communication between the Ellenburger aquifer and the basement, we increase the along-fault permeability, which models the presence of fractures located in the fault damage zone and oriented along the fault [5, 48, 1]. To model the permeability enhancement due to hydraulic fracturing in the Barnett shale gas wells, we altered the permeability field around each gas well’s horizontal open-hole section. Values for the permeability en- hancement factor in the hydraulically stimulated volume, which is 3D and includes the vertical direction, are determined during model calibration discussed below. The main observation from Figure 7.4e is that the enhanced permeability decreases with distance from the wellbore. This is because the effect of fractur- ing ceases as we move away from the well. As a result, the bottom layer within Barnett formation remains at the original Barnett permeability value. The permeability enhancement does not extend to the two faults, which is expected from the average effective length of hydraulic fractures reported in the literature, which is in the range of 50-75 m approximately [47]. Table 7.1 lists poromechanical properties. Generally, no-flux boundary condition is prescribed on the flow domain boundaries except for portions of south and west boundaries of Ellenburger layers where in order to investigate the sensitivity of DFW- Irving model to Basin model pressure evolution, a variable pressure boundary condition was prescribed. The pressure in the variable pressure boundary condition was prescribed by the evolving pressure of Basin model (Figure 7.2). This was achieved by defining two imaginary injection wells, one in each variable pressure boundary, with a defined evolving bottom pressures defined by the pressure of the same location in the Basin model. in order to mimic the effect of pressure front, we increased the horizontal permeability the row of cells where we imposed south variable pressure boundary and in the column of cells where we imposed west variable pressure boundary. Thus, allowing pressure to diffuse in a line front fashion. The vertical boundaries on east, west, north and south of the geomechanical domain are roller boundaries with 102 zero normal displacement. The bottom boundary is fixed with zero displacement in the three directions, and the top boundary is a traction-free surface. We used an initial pressure gradientdp i /dz of 11.28 kPa/m, minimum horizontal stress gradientdσ hi /dz of 4.5 kPa/m, maximum horizontal stress gradientdσ Hi /dz of 4.95 kPa/m, and vertical stress gradient ofdσ vi /dz of 14.7 kPa/m. 7.3.3 CouplingBendArch-FortWorthBasinandDFW-Irvingmodels Figure 7.3 shows how we managed to couple the two models by passing pressure from Bend Arch-Fort Worth Basin model to the DFW-Irving model. The pressure in the Basin model is driven by injection in Basin while the pressure in the DFW-Irving model is driven by production and injection in Barnett and Ellenburger formation in addition to the boundary pressure resolved from the Basin model. In order to feed the Basin Model pressure output into DFW-Irving Model, we defined two high permeability sets of cells located at the south and west edges of the DFW-Irving model with two pseudo water injectors located at the center of each set. The function of these two sets of cells is to emulate the effect of Basin model pressure acting on the DFW-Irving model. The bottomhole pressure in the high permeability edges in the DFW-Irving model is prescribed through the pressure performance of the cell that shares the same location in the Basin model as shown in Figure 7.3. 7.3.4 ModelCalibration To calibrate the model, we start with the well representation. We modeled the horizontal wells with an open-hole completion interval across five cells in the mesh. It provides a better agreement between pre- dicted and observed gas and water rates and volumes. To mimic the effect of hydraulic fracturing, as wells come online during simulation, we alter the permeability field around each gas producer wellbore in a logarithmic-with-distance fashion. Starting with a near-wellbore permeability of 9.86× 10 − 14 m 2 , the 103 Table 7.1: DFW-Irving Model parameters including the initial stress and pressure gradient values.h is the vertical thickness of a layer or a fault. E is Young’s modulus. ν is Poisson’s ratio. b w is Biot coefficient. ψ is internal friction angle. c is cohesion. h (km) E (kPa) ν b w ψ ( ◦ ) c (kPa) Overburden 1.8 1.44× 10 7 0.2 0.7 30 2× 10 4 Marble Falls 0.1 6× 10 7 0.2 Barnett 0.1 4× 10 7 0.23 Ellenburger 1 6× 10 7 0.2 Basement 2 4.3× 10 7 0.27 Fault 1 3.99 4× 10 7 0.25 0 Fault 2 5.64 Fault 3 5.64 Fault 4 6 Fault 5 6.58 Fault 6 6.8 Fault 7 1.78 Fault 8 6.4 Fault 9 5.53 permeability of the adjacent cell as we move away from the wellbore decreases by a factor of 10 (Fig- ure 7.4e). Next, we move up to the field scale where quantities such as the initial water saturation ( S wi ) in Barnett are uncertain due to the limited information available in the literature. Based on geographic prox- imity and geologic continuity, we estimate that formation properties in the DFW-Irving area are analogous to those in the Azle area, which has been studied more extensively in the literature. We defined Barnett’s initial water saturation S wi as 0.5 and assume Ellenburger to be completely saturated with water. Next, we investigated the sensitivity of the DFW-Irving model to the pressure boundary condition values due to coupling to the basin model. The final result of model calibration is shown in Figure 7.5 by comparing the observed (circles) and simulated (lines) monthly gas production rateq g and monthly water injection rate q wi . The excellent agreement between simulated and observed values supports calibration and suggests confidence in the model. 104 Fault 1 Fault 2 Fault 3 Fault 4 Fault 5 Fault 6 Fault 8 Fault 9 Fault 7 8 6 4 2 0 35 30 25 45 40 20 35 30 15 25 20 10 15 5 10 5 0 0 Figure 7.2: (a) 3D view of the Bend Arch-Fort Worth basin flow and geomechanics domains. The domain is bounded by Muenster arch in the northeast and by Ouachita thrust front in the southeast. (b) 3D view of the DFW-Irving flow and geomechanics domains showing Barnett formation (green) Ellenburger formation (blue) and faults ((red) for portions that are represented by flow and geomechanis mesh and (purple) for portions that are represented by geomechanis mesh only). Triangles (red) represent the surface location of the gas and brine producers and inverted triangles (blue) represent the surface location of water injection wells. 105 0 200 400 600 800 1000 1200 1400 1600 0 200 400 600 800 1000 2005 2007 2009 2011 2013 2015 2017 2019 0 200 400 600 800 1000 Figure 7.3: Coupling of the two models, (a) Bend Arch-Fort Worth basin and (b) DFW-Irving model. Injection-induced pressure front in the basin model propagates to western and southern boundaries of the Ellenburger formation in the DFW-Irving model. Coupling ensures that the pressure boundary con- dition for the faults in the DFW-Irving model evolves in sync with the long-term injection activity in the basin.. 7.3.5 ProductionandInjectionHistory Barnett production and Ellenburger Injection in the vicinity of DFW-Irving area started in the year 2007 and the earthquakes occurred in the years 2008-2015 in DFW and in the year 2015 in Irving. We consider the 21-year period from January 1, 2000, to December 31, 2020, to conduct our analysis. Dynamics of the Bend Arch–Fort Worth Basin model is driven by the reported water injection volumes for the seven coun- ties that are located west of Muenster Arch fault and Ouachita Thrust Front fault as shown in Figure 7.1. Dynamics of the DFW–Irving model is driven by gas and water production history from 170 horizon- tal and hydraulically fractured wells in Barnett Shale and two Ellenburger water injection wells shown in Figure 7.3c in addition to boundary pressure prescribed by Basin model. All the well data including monthly production and injection volumes and well locations are obtained from the IHS Markit database and provided in the manuscript and the supplementary material. Figure 7.5a and b show DFW-Irving model’s monthly produced gas and injected water. The volume of the total gas produced is more than 7× 10 9 m 3 and the volume of total water injected is more than 3× 10 6 m 3 as shown in Figure 7.5c. Figure 7.5c also shows the cumulative volume of Basin water injected which 106 10 15 20 25 30 35 40 45 0 5 10 15 20 25 30 35 5 0 -20 -19 -18 -17 -16 -15 -14 10 10 10 10 10 10 10 0.045 0.05 0.055 0.06 0.065 0 0.01 0.02 0.03 0.04 0.05 0.06 10 -19 0 0.05 0.1 0.15 0.2 0.25 0.987 1.974 2.961 3.948 4.935 5.922 6.908 7.895 0 10 -14 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 1.5 2 2.5 3 3.5 4 4.5 5 0.05 0.055 0.06 0.065 0.07 0 0.01 0.02 0.03 0.04 0.05 0.06 Figure 7.4: (a) Histogram of Barnett porosity, (b) histogram of Barnett log-permeabilities log k x in I di- rection and logk y in J direction, (c) histogram of Ellenburger porosity, (d) histogram of Ellenburger log- permeability in I and J directions. Vertical-to-horizontal permeability ratiok z /k x is 0.1. (e) Top view of the Barnett permeability map showing the effect of hydraulic fracturing on enhancing the permeability field around each horizontal gas producer. 107 2005 2007 2009 2011 2013 2015 2017 2019 0 0.5 1 1.5 2 2.5 3 3.5 10 6 0 1 2 3 4 5 6 7 8 10 4 (a) (b) 10 4 10 5 10 6 10 7 10 8 10 9 (c) Figure 7.5: Observed (circles) and simulated (lines) (a) gas monthly production rateq g , (b) water monthly injection rate q wi , (c) cumulative Basin and DFW water injection volume Q wi in addition to the south and west water influx in the DFW-Irving field plotted with time. All volumes are at surface pressure and temperature conditions. is close to 2× 10 8 m 3 in addition to the DFW-Irving model cumulative south water influx which is close to 7× 10 6 m 3 and cumulative west water influx which is close to 3 × 10 6 m 3 . 7.4 Results 7.4.1 BendArch-FortWorthBasinmodelpressureandstressdistributions Figure 7.6 shows basin overpressure∆ p with respect to initial pressure at time close to the time of (a) DFW seismic events, that occurred towards the end of 2008, and (b) Irving seismic events time, that occurred during 2015. This figure shows that the water disposal activities in basin lead to an overall increase in basin pressure, especially in the region south of Irving between DFW fault and Ouachita thrust front. The relatively higher disposed water volume between DFW fault an Ouchita thrust front in addition extremely low fault permeability across thess faults channeled pressure front towards Irving Area. Figure 7.6 (c) and (d) shows Basin change in effective mean stress, ∆ I ′ 1 , at the same timesteps. ∆ I ′ 1 measures the injection-induced volumetric expansion of Basin which evolves in sync with Basin∆ p shown in Figure 7.6 (a) and (b). 108 0 200 400 600 800 1000 1200 1400 1600 -2500 -2000 -1500 -1000 -500 0 (a) 01-Jan-2009 (c) 01-Jan-2009 (b) 01-Jan-2016 (d) 01-Jan-2016 Figure 7.6: (a) and (b) Basin overpressure∆ p at the selected time steps. (c) and (d) change in the effective mean stress∆ I ′ 1 (with respect to the initial value). Ouachita thrust front bounding the basin from southeast and low horizontal permeability across DFW airport fault helps in channeling pressure front towards Irving area. 109 7.4.2 DFW-Irvingmodelpressureandstressdistributions Figure 7.7 shows the Barnett and Ellenburger∆ p with respect to initial pressure at four selected timesteps starting around the time of DFW seismic events until January 1 st , 2021. Barnett depletion and pore pres- sure distributions is affected by hydraulic fracturing-induced permeability enhancement around each gas producer horizontal wellbore. Gas production from Barnett in this model is localized mainly to the west side and in DFW airport area (between fault 1 and 3). On the other hand, Ellenburger∆ p evolves more uniformly due to a more homogenous permeability field. looking at ∆ p evolution in Ellenburger formation in figure 7.7 (a), we can see the localized effect of water injection between fault 1 and 3 and fault 3 and 4 on pore pressure. In addition to the two water injectors near DFW airport, south and west pressure fronts caused by water injection in the basin causes pressure to increase. Red arrows shown in Figure 7.7 (d) in Ellenburger indicates the direction at which Basin pressure front evolves due to water injection across the basin. We can see that the effect of south pressure front is more prominent due to the higher volume of water injected south of Irving. Ellenburger injection has only a minor effect on Barnett’s pressure because Barnett’s low permeability slows down the propagation of vertical pressure diffusion front. Figure 7.8 shows the effective mean stress, ∆ I ′ 1 , which measures the production and injection-induced volumetric contraction and expansion. Looking at Barnett ∆ I ′ 1 , we can see how it evolves in sync with ∆ p shown in figure 7.7, i.e. more depletion corresponds to more compressive stresses and higher ∆ I ′ 1 values. Ellenburger∆ I ′ 1 field generally evolves in sync with ∆ p. i.e more recharge corresponds to more tensile stresses and lower ∆ I ′ 1 . Looking at figure 7.8 (d), we can see some areas in Barnett ,such as the one circled with dashed red line, that have higher ∆ I ′ 1 .This behavior can be explained as production- induced compressive stresses in Barnett have propagated deeper into Ellenburger and can counterbalance the injection-induced tensile stresses in Ellenburger. 110 0 5 10 15 20 25 30 35 40 15 20 25 30 (a) 01-Jan-2009 (b) 01-Apr-2013 (c) 01-Jan-2016 (d) 01-Jan-2021 -2000 -1000 0 0 400 800 Figure 7.7: Barnett (top) and Ellenburger (Bottom) overpressure∆ p at selected time steps: a, b, c, and d. While fracturing-induced permeability enhancement around gas producers lead to heterogeneity in the pressure field in Barnett, ∆ p evolves more uniformly due to Ellenburger uniform permeability. The effect of Basin pressure front is clearly visible as∆ p evolves more rabidly due to the generally low permeability of faults inx direction and specifically fault 5 and fault 8 channeling pressure front toward fault 6. 111 0 200 400 600 800 0 5 10 15 20 25 30 35 40 45 0 10 20 30 -400 -300 -200 -100 0 (a) 01-Jan-2009 (b) 01-Apr-2013 (c) 01-Jan-2016 (d) 01-Jan-2021 Figure 7.8: Barnett (top) and Ellenburger (Bottom) change in the effective mean stress ∆ I ′ 1 (with respect to the initial value) at the selected time steps: a, b, c, and d.∆ I ′ 1 , which is positive in compression and negative in tensile, shows a generally monotonic response to∆ p in Figure 7.7. Some areas, such as area highlighted with red dashed line, do not directly follow Ellenburger’s pressure and have a smaller magnitudes of∆ I ′ 1 (closer to zero) which indicates stress propagation from Barnett. This shows that injection-induced tensile stresses are counterbalanced by Barnett’s production-induced compressive stresses. 112 7.4.3 Interactionandpropagationofstressesalongdepthandwithtime Figure 7.9 (b), (c), and (d) shows three depth profiles of ∆ p, ∆ I ′ 1 , ∆ √ J 2 , andu z at locations highlighted in figure 7.9 (a). The profiles are taken at locations that experienced the highest amount of ∆ I ′ 1 evolution in Barnet and Ellenburger. The profiles are shown at four time steps to understand time evolution. we can see that pressure is depleting in Barnett due to production in figure 7.9 (b), (c), while Ellenburger is experiencing pressure increase due to water injection near DFW airport as shown in (c) and due to boundary Basin pressure as shown in (b) and (d). by comparing (b) and (c) we can see that the magnitude and effect of the south pressure boundary is more sever than the west pressure boundary. This translates into localized compression in Barnett and overall tensile in Ellenburger, which is observed in the change of the effective mean stress ∆ I ′ 1 and the von Mises stress∆ √ J 2 . The latter measures the shear or distortion response in the rock; Eq. (7.1). comparingu z for the three depth profiles ,we can see that both (b) and (c) are experiencing net subsidence at the surface.Figure 7.9 also shows how the meshing algorithm was able to honor changing formation thickness and elevation at different locations by comparing the elevations and thickness for Barnet and Ellenburger at (b), (c), and (d). σ =σ v 1+S, ε=ε v 1+e, (7.1) 7.4.4 Stressrotation Figure 7.10 shows the Dec 2020 orientations of the minimum effective principal stress σ ′ 3 and the interme- diate effective principal stress σ ′ 2 for all cells that exhibited a rotation greater than 0.5 degree in (a) Barnett and (b) Ellenburger, respectively. Studying the spatial distribution of the principal stress direction rotation, we can notice that there are two distinct regions where the rotation is sever. The first region is around the hydraulically fractured area in Barnett as production causes pressure depletion and compressive stresses 113 -2000 0 1000 2000 3000 4000 5000 6000 7000 8000 0 400 800 0 500 1000 0 10 -3 5 -10 10 -3 1000 2000 3000 4000 5000 6000 7000 8000 -2000 0 0 400 800 0 500 1000 0 10 -3 5 -10 10 -3 -2000 0 0 400 800 0 500 1000 0 10 -3 5 -10 10 -3 Figure 7.9: (a) Overpressure∆ p in Barnett and Ellenburger at the end of simulation time. (b), (c), and (d) are depth profiles of ∆ p, changes in the effective mean stress ∆ I ′ 1 and von Mises stress ∆ √ J 2 , and u z along a column of cells in the selected location shown in (a). 114 10 15 20 25 30 35 40 0 5 10 15 20 25 30 0 200 400 600 800 1000 -2500 -2000 -1500 -1000 -500 0 5 0 Figure 7.10: Rotation of the minimum effective principal stress σ ′ 3 and the intermediate effective principal stressσ ′ 2 on December 2020 in (a) Barnett and (b) Ellenburger. within the area. Due to the surrounding unfractured ultra low matrix impeding pressure communication and localizing pressure depletion, tensile stresses develops in the unfractured surrounding area. The sec- ond region is along the faults, this is most clearly around faults 5 through 8, where south pressure front is acting. This is due to the faults acting as a barrier impeeding pressure communication between the re- gions directly in communication with pressure boundary and surrounding regions. This phenomenon of induced anisotropy and heterogeneity in the stress state is most clearly observed in ultra-low permeability reservoirs undergoing hydraulic fracturing as it creates a clear contrast in permeability between fractured and un-fractured rocks [13, 12]. This information is needed to optimize the trajectories of new and infill horizontal wells and their hydraulic fracturing jobs. Consequently, the figure shows the importance of modeling spatiotemporal changes in principal stress magnitudes and directions. 7.4.5 Faultstressesandinducedfailuremechanisms 7.4.5.1 DFW Figure 7.11 (a) shows the Ellenburger overpressure∆ p at December 2020 with the location of fault-4 high- lighted in red which underlies DFW area. A 3D view of fault-4 showing failure locations identified by ε p 115 is shown in Figure 7.11 (b). This shows that failure in fault-4 is limited only to Ellenburger formation and did not propagate downward in the basement. Figure 7.11 (c) shows the evolution of overpressure ∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure strain dε p for the column of cells highlighted in figure 7.11 (b). In this and subsequent plots, we use the mesh axes (I, J, K), which are quasi-parallel to Cartesian axes (X,Y,Z), to indicate the location of mesh cells being analyzed in a plot. The location of the fault cells column, K=16 through 20, is closest to the adjacent injector. failuredε p occurred right after pore pressure increase indicated by∆ p caused by the two water injectors coming online. Moredε p is noticeable in shallower layers (K=16:18) in comparison to (K=19:20). as water injection is ceased from the injector nearest to the cells column, pressure declines hindering any further failure, however, looking at pressure evolution on later times, we can see that it increases due to overall basin pressure increase. Figure 7.11 (d) shows the stress path of the same cells. The paths are shown in the stress invariant space ofI ′ 1 vs. √ J 2 , and the direction of time is implied by the curved arrow shown in the plot. Stress path of cells K=16 to K=20 generally follows the same pattern. The stress path of all cells moves towards failure envelop staring at the time the far injector comes onlinet WI1 . As the close injector comes online t WI2 the stress path continue moving toward failure with much higher rate. The onset of failure is marked by the asterisk symbol. As water injection ceasest WI SI , stress paths for all fault move toward more stable region. 7.4.5.2 Irving Figure 7.12 (a) shows the Ellenburger overpressure∆ p at December 2020 with the location of fault-6 high- lighted in red which underlies Irving area. A 3D view of fault-6 showing failure locations identified by ε p is shown in (b). This shows that, on the contrary to the failure in fault-4, failure in fault-6 propagates downward to the basement. Figure 7.12 (c) shows the evolution of overpressure∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure straindε p for key cells in the 116 0 400 800 0 1 2 3 4 x10 -6 K=16:20 2005 2007 2009 2011 2013 2015 2017 2019 2021 -100 0 100 200 300 0 2 4 6 x10 -7 K=16 K=17 K=18 K=19 K=20 2.63 2.632 2.634 2.636 2.638 2.64 2.642 2.644 10 4 1.887 1.888 1.889 1.89 1.891 1.892 1.893 1.894 10 4 time K=16 K=17 K=18 K=19 K=20 Figure 7.11: (a) Ellenburger∆ p at December 2020 highlighting the location of fault-4 (red) (b) 3D view of fault-4 showing failure magnitudeε p (c) Evolution of overpressure∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure straindε p in fault cells located at the depths of Ellenburger layers. K=16 through 20 are top five layers in Ellenburger. ∆ p,∆ √ J 2 , and∆ I ′ 1 are plotted against the left axis whiledε p is plotted against the right axis. Axis labels, shown only in the bottom left plot, are identical for all other plots. (d) Stress trajectories of fault cells in theI ′ 1 − √ J 2 stress invariant space over the 2005-2020 period. The K=16 cell is chosen as the reference cell from which other cells are compared. Hence, all cells start at the same initial point, which corresponds to the initial stress state of K=16. The start time of gas production from Barnett (t GP ), the start time of water injection in Ellenburger (t WI1 ) and (t WI2 ) corresponding to the two water injectors, the failure onset time (t p i ), and the injector shut-in (t WI SI ) are marked. 117 highlighted column in figure 7.12 (b) that are located in the overburden (K=15), Barnett (K=19), Ellenburger (K=26,28), and basement (K=36). we can see that failuredε p is only noticeable in fault cells located in El- lenburger and basement and occurred after pore pressure increased due to overall basin pressure increase. Figure 7.12 (d) shows the stress path of the same key fault cells highlighted in figure 7.12 (c). We can see that the stress paths for the cells can be categorized based on their general behavior which is found to be related to their locations. the first category is overburden cells, where we can see no failure was noticeable. the second category is the Ellenburger, where both∆ I ′ 1 and∆ √ J 2 are decreasing causing the stress path to move along failure envelope until failure point. The third category is Basment cells, where ∆ I ′ 1 was decreasing and∆ √ J 2 was decreasing causing the stress path to move in a collision path toward the failure envelop until failure point then continue with the slope of failure envelop. To more understand the three stress path behavior categories, we examined the change in vertical stress inX direction that evolves in timeδ X σ V (Eq. (7.2)). This sheds some light on how fault cells responds to changing in pore pressure in comparison to the adjacent matrix cells. we can see that Ellenburger formation experience more tensile than the fault attributed to fault’s low horizontal permeability impeding pressure from entering faults cells. on the other hand we can see that the fault cells experienced more tensile in comparison to basement cells due to absence of basement permeability isolating pressure to the fault cells only. δ X σ V =σ V (I+1)− σ V (I) (7.2) 7.4.5.3 Potency-MagnitudeScalingandMomentTensorcomparison Figure 7.13(a) shows the observed and calculated failure events resulting from employing the two above mentioned methods. we can see that Potency-Magnitude Scaling method provided higher degree of agree- ment with the observed date for both fault-4 and fault-6. 118 time 2005 2007 2009 2011 2013 2015 2017 2019 2021 0 400 800 0 0.5 1 1.5 2 10 -5 1000 2000 3000 4000 5000 6000 7000 8000 0 100 -100 -200 0 200 400 600 0 1 2 3 10 -7 K=15 K=26 K=28 K=36 K=19 2.615 2.62 2.625 2.63 2.635 2.64 2.645 10 4 1.882 1.884 1.886 1.888 1.89 1.892 1.894 1.896 10 4 15 19 26 28 36 Figure 7.12: (a) Ellenburger∆ p at December 2020. (b) Fault 4 failure magnitudeε p . (c) Evolution of over- pressure∆ p, change in von Mises stress∆ √ J 2 , change in effective mean stress ∆ I ′ 1 , and instantaneous failure straindε p in different fault cells located at overburden, Barnett, Ellenburger and basement layers. K=15 is located in overburden, K=19 is located in Barnett, K=26 and 28 are located in Ellenburger, and K=36 is located in basement.(d) Stress trajectories of fault cells in theI ′ 1 − √ J 2 stress invariant space over the 2005-2020 period. The failure onset time (t p i ) is marked. (e) fault depth profile of how change in vertical stress inX direction evolves in timeδ x σ v 119 One measure that is used to differentiate induced seismic events from naturally occurring events is Gutenberg–Richter (GR) [28] earthquake frequency magnitude relationship B-value which correlate the log of number of events N greater than or equal to the magnitude M versus failure magnitude M. Fig- ure 7.13 (b) and (c) shows the b value for the observed and calculated events based on Potency-Magnitude Scaling method. 120 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0.5 1 1.5 2 2.5 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 1.6 3.6 Figure 7.13: (a) Temporal distribution of simulated and observed failure events showing fault-4 observed failure events (black asterisks), fault-6 observed events (red asterisks), fault-4 simulated failure events (blue) fault-6 simulated failure events (yellow), simulated failure events calculated using Moment Tensor are marked by triangles, and simulated failure events calculated using Potency-Magnitude Scaling are marked by circles. (b) and (c) log of number of eventsN greater than or equal to the magnitudeM versus failure magnitudeM for observed (asterisks) and simulated (circles) failure events occurred on fault-4 and fault-6 respectively. 121 Chapter8 Conclusions Determining the causality of induced seismicity at sites with segregated production and injection units is more challenging than determining it at sites where production and injection are from one hydraulically connected unit. The reason is the difference in the style and magnitude of the injection- and production- induced stresses and how they interact within a given reservoir-basement complex. This interaction affects the stability of a fault intersecting or bounding these regions and must be understood to improve our forecasts of the seismicity hazard, especially in regions with low natural seismicity or strain rate, e.g., central US. We studied the problem of modeling production- and injection-induced stress evolution in low porosity low permeability faulted reservoirs and investigated the role of induced stresses in activating basement faults. A coupled multiphase flow—geomechanics–fault reactivation modeling framework has been proposed, successfully implemented, and presented. The framework has been validated and applied on a number of cases with increasing complexities starting with a simple synthetic case to a real field with one seismic cluster and two connected faults to a real field with two seismic clusters with nine curved faults. A novel MATLAB structured meshing routine that is able to capture the geometric complexity of reser- voirs with multiple, curved, and intersecting faults is successfully demonstrated. The meshing routine combines the model building and runtime advantages of a structured mesh with the geometry-conforming 122 capability of unstructured meshing routines. This can increase the geologic fidelity in subsurface models and decrease the simulation runtime of the coupled flow-geomechanics models, which will benefit the petroleum engineering and computational geoscience communities. One major achievement in our research is the use of seismic event magnitude, location, and time to calibrate coupled flow and geomechanics models, which are often based on the quasi-static equilibrium assumption. This assumption, which is suitable for pre and post seismic periods dominated by elastic loading processes and leads to computational gains compared to full elastodynamics simulation, does not allow seismic velocity rupture and radiation. Therefore, we require a model to compute seismic events from plastic failure events in our coupled model. We used and compared two methods: (1) the moment tensor method and (2) the potency-magnitude scaling method. The moment tensor method was originally developed for seismic events of magnitude 5 and above, which might lead to some inaccuracy as most induced seismic events are of magnitude 4 and below. Also, this method requires the knowledge of material properties at the hypocenter of the event, which is difficult to estimate. On the other hand, the potency- magnitude scaling method was developed for a much wider range of event magnitudes and does not depend on the assumption of material properties at the source. In our studies, we show that the potency-magnitude scaling method performs better than the moment tensor method in terms of calibrating the model to the observed data. Optimally placing horizontal producer wells and planning a hydraulic fracturing program to create most connected fracture networks are of utmost importance in developing low permeability reservoirs such as shales. Magnitudes and directions of principal stresses are important parameters for that purpose. However, their values change with time because production from fractured wells leads to redistribution and rotation of the stress state. We presented a method to map spatiotemporal changes in principal stress magnitudes and orientations due to nearby production or injection. This can proactively provide stress information before drilling a new well. 123 8.1 Recommendations 1. Couplingwithelastodynamics Incorporating elastodynamics or approximating its effect via the quasi-dynamics framework [54] during fault failure, instead of assuming quasi-static equilibrium, will better capture the physics of seismic fail- ure process. The additional computational cost of this enhancement is justified by the improvement in accuracy of stress solution and the event magnitude calculation. This will also allow for a better model calibration with observed seismic events. 2. Expandingthemeshingalgorithm Expanding the meshing algorithm’s spatial degrees of freedom to capture more complex stratigraphy such as domes and pinch-out will allow the proposed framework to model a wider range of geologic sites. 3. Additionaltestingofthepotency-magnitudescalingmethod The potency-magnitude scaling method, which has been derived based on small-to-moderate sized earth- quakes, proved to be more capable in modeling induced failure events in the DFW-Irving case study. The method should be tested on more induced seismicity models to evaluate its suitability and superiority over the moment tensor method. 124 Bibliography [1] Fabrizio Agosta, M. Prasad, and A. Aydin. “Physical properties of carbonate fault rocks, fucino basin (Central Italy): Implications for fault seal in platform carbonates”. In: Geofluids 7.1 (2007). issn: 14688115.doi: 10.1111/j.1468-8123.2006.00158.x. [2] Matteo Albano, Salvatore Barba, Gabriele Tarabusi, Michele Saroli, and Salvatore Stramondo. “Discriminating between natural and anthropogenic earthquakes: Insights from the Emilia Romagna (Italy) 2012 seismic sequence”. In: Scientific Reports 7.1 (2017).issn: 20452322.doi: 10.1038/s41598-017-00379-2. [3] Khalid Aziz and Antonín Settari. Petroleum Reservoir Simulation. Elsevier, 1979.isbn: 9780853347873. [4] George Backus and Marjorie Mulcahy. “Moment tensors and other phenomenological descriptions of seismic sources—II. Discontinuous displacements”. In: Geophysical Journal of the Royal Astronomical Society 47.2 (1976).issn: 1365246X.doi: 10.1111/j.1365-246X.1976.tb01275.x. [5] C. A. Barton, M. D. Zoback, and D. Moos. “Fluid flow along potentially active faults in crystalline rock”. In: Geology 23.8 (1995).issn: 00917613.doi: 10.1130/0091-7613(1995)023<0683:FFAPAF>2.3.CO;2. [6] Yehuda Ben-Zion. “Appendix 2 Key formulas in earthquake seismology”. In: International Geophysics 81.PART B (2003).issn: 00746142.doi: 10.1016/S0074-6142(03)80304-2. [7] Yehuda Ben-Zion and Vladimir Lyakhovsky. “Representation of seismic sources sustaining changes of elastic moduli”. In: Geophysical Journal International 217.1 (2019).issn: 1365246X.doi: 10.1093/gji/ggz018. [8] Yehuda Ben-Zion and Lupei Zhu. “Potency-magnitude scaling relations for Southern California earthquakes with 1.0 < ML < 7.0”. In: Geophysical Journal International 148.3 (2002).issn: 0956540X.doi: 10.1046/j.1365-246X.2002.01637.x. [9] Per Eirik S. Bergmo, Alv Arne Grimstad, and Erik Lindeberg. “Simultaneous CO2 injection and water production to optimise aquifer storage capacity”. In: International Journal of Greenhouse Gas Control 5.3 (2011).issn: 17505836.doi: 10.1016/j.ijggc.2010.09.002. 125 [10] Maurice A. Biot. “General theory of three-dimensional consolidation”. In: Journal of Applied Physics 12.2 (1941).issn: 00218979.doi: 10.1063/1.1712886. [11] Abdulrahman Bubshait and Birendra Jha. “A novel meshing routine for modeling reservoirs with multiple curved and intersecting faults in coupled flow-geomechanical dual 3D grids”. In: SPE Western Regional Meeting Proceedings. Vol. 2019. 2019.doi: 10.2118/195256-ms. [12] Abdulrahman Bubshait and Birendra Jha. “Revisiting 2013-2014 Azle seismicity to understand the role of Barnett production on stress propagation and fault stability”. In: Geophysics (2021). [13] Ahmed Bubshait and Birendra Jha. “Coupled poromechanics-damage mechanics modeling of fracturing during injection in brittle rocks”. In: International Journal for Numerical Methods in Engineering 121.2 (2020).issn: 10970207.doi: 10.1002/nme.6208. [14] Thomas A. Buscheck, Yunwei Sun, Mingjie Chen, Yue Hao, Thomas J. Wolery, William L. Bourcier, Benjamin Court, Michael A. Celia, S. Julio Friedmann, and Roger D. Aines. “Active CO 2 reservoir management for carbon storage: Analysis of operational strategies to relieve pressure buildup and improve injectivity”. In: International Journal of Greenhouse Gas Control 6 (2012).issn: 17505836. doi: 10.1016/j.ijggc.2011.11.007. [15] Rongqiang Chen, Xu Xue, Jaeyoung Park, Akhil Datta-Gupta, and Michael J. King. “New insights into the mechanisms of seismicity in the Azle area, North Texas”. In: Geophysics 85.1 (2020).issn: 19422156.doi: 10.1190/geo2018-0357.1. [16] CMG. CMG User Manual. 2019. [17] Olivier Coussy. Mechanics of Porous Continua. John Wiley and Sons, 1995. [18] Steve Everley. Irving Earthquake Study Relies on Contradictions and Dubious Assumptions. July 2016.url: https://www.energyindepth.org/irving-earthquake-study-contradictions-dubious- assumptions/?154 (visited on 12/29/2021). [19] Zhiqiang Fan, Peter Eichhubl, and Julia F.W. Gale. “Geomechanical analysis of fluid injection and seismic fault slip for the Mw4.8 Timpson, Texas, earthquake sequence”. In: Journal of Geophysical Research: Solid Earth 121.4 (2016).issn: 21699356.doi: 10.1002/2016JB012821. [20] J. H. Frantz, J. R. Williamson, W. K. Sawyer, D. Johnston, G. Waters, L. P. Moore, R. J. MacDonald, M. Pearcy, S. V. Ganpule, and K. S. March. “Evaluating barnett shale production performance using an integrated approach”. In: SPE Eastern Regional Meeting. 2005.doi: 10.2523/96917-ms. [21] Cliff Frohlich, Chris Hayward, Brian Stump, and Eric Potter. “The Dallas-Fort Worth earthquake sequence: October 2008 through May 2009”. In: Bulletin of the Seismological Society of America 101.1 (2011).issn: 00371106.doi: 10.1785/0120100131. [22] Beatriz Gaite, Arantza Ugalde, Antonio Villaseñor, and Estefania Blanch. “Improving the location of induced earthquakes associated with an underground gas storage in the Gulf of Valencia (Spain)”. In: Physics of the Earth and Planetary Interiors 254 (2016).issn: 00319201.doi: 10.1016/j.pepi.2016.03.006. 126 [23] J. Geertsma. “The Effect of Fluid Pressure Decline on Volumetric Changes of Porous Rocks”. In: Transactions of the AIME 210.01 (1957).issn: 0081-1696.doi: 10.2118/728-g. [24] Peter O. Gold, Eric Cowgill, Oliver Kreylos, and Ryan D. Gold. “A terrestrial lidar-based workflow for determining three-dimensional slip vectors and associated uncertainties”. In: Geosphere 8.2 (2012).issn: 1553-040.doi: 10.1130/GES00714.1. [25] Pablo J. González, Kristy F. Tiampo, Mimmo Palano, Flavio Cannavó, and José Fernández. “The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading”. In: Nature Geoscience 5.11 (2012).issn: 17520894.doi: 10.1038/ngeo1610. [26] E. Gringarten, B. Arpat, A. Haouesse, A. Dutranois, L. Deny, S. Jayr, A. L. Tertois, J. L. Mallet, A. Bernal, and L. Nghiem. “New grids for robust reservoir modeling”. In: Proceedings - SPE Annual Technical Conference and Exhibition. Vol. 6. 2008.doi: 10.2118/116649-ms. [27] E. Gringarten, A. Haouesse, B. Arpat, and L. Nghiem. “Advantages of using vertical stair step faults in reservoir grids for flow simulation”. In: SPE Reservoir Simulation Symposium Proceedings. Vol. 2. 2009.doi: 10.2118/119188-ms. [28] B. Gutenberg and C. F. Richter. “Frequency of earthquakes in California”. In: Bulletin of the Seismological Society of America 34.4 (1944).issn: 0037-1106.doi: 10.1785/bssa0340040185. [29] Thomas C. Hanks and Hiroo Kanamori. “A moment magnitude scale”. In: Journal of Geophysical Research B: Solid Earth. Vol. 84. B5. 1979.doi: 10.1029/JB084iB05p02348. [30] Jeanne L. Hardebeck and Peter M. Shearer. “A new method for determining first-motion focal mechanisms”. In: Bulletin of the Seismological Society of America 92.6 (2002).issn: 00371106.doi: 10.1785/0120010200. [31] Peter H. Hennings, Jens Erik Lund Snee, Johnathon L. Osmond, Heather R. Deshon, Robin Dommisse, Elizabeth Horne, Casee Lemons, and Mark D. Zoback. “Injection-induced seismicity and fault-slip potential in the fort worth basin, Texas”. In: Bulletin of the Seismological Society of America 109.5 (2019).issn: 19433573.doi: 10.1785/0120190017. [32] Matthew J. Hornbach, Heather R. Deshon, William L. Ellsworth, Brian W. Stump, Chris Hayward, Cliff Frohlich, Harrison R. Oldham, Jon E. Olson, M. Beatrice Magnani, Casey Brokaw, and James H. Luetgert. “Causal factors for seismicity near Azle, Texas”. In: Nature Communications 6 (2015).issn: 20411723.doi: 10.1038/ncomms7728. [33] Matthew J. Hornbach, Madeline Jones, Monique Scales, Heather R. DeShon, M. Beatrice Magnani, Cliff Frohlich, Brian Stump, Chris Hayward, and Mary Layton. “Ellenburger wastewater injection and seismicity in North Texas”. In: Physics of the Earth and Planetary Interiors 261 (2016).issn: 00319201.doi: 10.1016/j.pepi.2016.06.012. [34] S. Horton. “Disposal of hydrofracking waste fluid by injection into subsurface aquifers triggers earthquake swarm in central arkansas with potential for damaging earthquake”. In: Seismological Research Letters 83.2 (2012).issn: 08950695.doi: 10.1785/gssrl.83.2.250. [35] John Conrad Jaeger and Neville G. W. Cook. No Title. London: Chapman and Hall, 1979. 127 [36] Birendra Jha and Ruben Juanes. “Coupled multiphase flow and poromechanics: A computational model of pore pressure effects on fault slip and earthquake triggering”. In: Water Resources Research 50.5 (2014).issn: 19447973.doi: 10.1002/2013WR015175. [37] R. Juanes, B. Jha, B. H. Hager, J. H. Shaw, A. Plesch, L. Astiz, J. H. Dieterich, and C. Frohlich. “Were the May 2012 Emilia-Romagna earthquakes induced? A coupled flow-geomechanics modeling assessment”. In: Geophysical Research Letters 43.13 (2016).issn: 19448007.doi: 10.1002/2016GL069284. [38] K. M. Keranen, M. Weingarten, G. A. Abers, B. A. Bekins, and S. Ge. “Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection”. In: Science 345.6195 (2014).issn: 10959203.doi: 10.1126/science.1255802. [39] Katie M. Keranen, Heather M. Savage, Geoffrey A. Abers, and Elizabeth S. Cochran. “Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence”. In: Geology 41.6 (2013).issn: 00917613.doi: 10.1130/G34045.1. [40] Xinglin Lei, Dongjian Huang, Jinrong Su, Guomao Jiang, Xiaolong Wang, Hui Wang, Xin Guo, and Hong Fu. “Fault reactivation and earthquakes with magnitudes of up to Mw4.7 induced by shale-gas hydraulic fracturing in Sichuan Basin, China”. In: Scientific Reports 7.1 (2017).issn: 20452322.doi: 10.1038/s41598-017-08557-y. [41] Roland W. Lewis and B. A. Schrefler. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. Second. Chichester: Wiley, 1998. [42] Maria Beatrice Magnani, Michael L. Blanpied, Heather R. DeShon, and Matthew J. Hornbach. “Discriminating between natural versus induced seismicity from long-term deformation history of intraplate faults”. In: Science Advances 3.11 (2017).issn: 23752548.doi: 10.1126/sciadv.1701593. [43] B. Maillot, S. Nielsen, and I. Main. “Numerical simulation of seismicity due to fluid injection in a brittle poroelastic medium”. In: Geophysical Journal International 139.2 (1999).issn: 0956540X.doi: 10.1046/j.1365-246X.1999.00933.x. [44] Patricia Martínez-Garzón, Yehuda Ben-Zion, Niloufar Abolfathian, Grzegorz Kwiatek, and Marco Bohnhoff. “A refined methodology for stress inversions of earthquake focal mechanisms”. In: Journal of Geophysical Research: Solid Earth 121.12 (2016).issn: 21699356.doi: 10.1002/2016JB013493. [45] Patricia Martínez-Garzón, Marco Bohnhoff, Grzegorz Kwiatek, and Georg Dresen. “Stress tensor changes related to fluid injection at the Geysers geothermal field, California”. In: Geophysical Research Letters 40.11 (2013).issn: 00948276.doi: 10.1002/grl.50438. [46] Patricia Martínez-Garzón, Grzegorz Kwiatek, Michèle Ickrath, and Marco Bohnhoff. “MSATSI: A MATLAB package for stress inversion combining solid classic methodology, a new simplified user-handling, and a visualization tool”. In: Seismological Research Letters 85.4 (2014).issn: 19382057.doi: 10.1785/0220130189. 128 [47] M. J. Mayerhofer, E. P. Lolon, J. E. Youngblood, and J. R. Heinze. “Integration of microseismic fracture mapping results with numerical fracture network production modeling in the Barnett shale”. In: Proceedings - SPE Annual Technical Conference and Exhibition. Vol. 2. 2006.doi: 10.2118/102103-ms. [48] Aderson F. do Nascimento, Rebecca J. Lunn, and Patience A. Cowie. “Modeling the heterogeneous hydraulic properties of faults using constraints from reservoir-induced seismicity”. In: Journal of Geophysical Research: Solid Earth 110.9 (2005).issn: 21699356.doi: 10.1029/2004JB003398. [49] Akemi Noda and Mitsuhiro Matsu’ura. “Physics-based GPS data inversion to estimate three-dimensional elastic and inelastic strain fields”. In: Geophysical Journal International 182.2 (2010).issn: 0956540X.doi: 10.1111/j.1365-246X.2010.04611.x. [50] Paul O. Ogwari, Heather R. DeShon, and Matthew J. Hornbach. “The Dallas-Fort Worth Airport Earthquake Sequence: Seismicity Beyond Injection Period”. In: Journal of Geophysical Research: Solid Earth 123.1 (2018).issn: 21699356.doi: 10.1002/2017JB015003. [51] Richard M. Pollastro, Daniel M. Jarvie, Ronald J. Hill, and Craig W. Adams. “Geologic framework of the Mississippian Barnett Shale, Barnett-Paleozoic total petroleum system, Bend arch-Fort Worth Basin, Texas”. In: American Association of Petroleum Geologists Bulletin 91.4 (2007).issn: 01491423. doi: 10.1306/10300606008. [52] Louis Andrew Quinones, Heather R. DeShon, Maria B. Magnani, and Cliff Frohlich. “Stress orientations in the fort worth basin, texas, determined from earthquake focal mechanisms”. In: BulletinoftheSeismologicalSocietyofAmerica 108.3 (2018).issn: 19433573.doi: 10.1785/0120170337. [53] Nicolas Remy, Alexandre Boucher, and Jianbing Wu. Applied Geostatistics with SGeMS. Cambridge: Cambridge University Press, 2009.isbn: 9781139150019.doi: 10.1017/CBO9781139150019. [54] J. R. Rice. “Spatio-temporal complexity of slip on a fault”. In: 98 (1993), pp. 9885–9907. [55] Pablo F. Sanz, Suvrat P. Lele, Kevin H. Searles, Sheng Yuan Hsu, Jorge L. Garzon, Jason A. Burdette, William E. Kline, Bruce A. Dale, and Paul D. Hector. “Geomechanical analysis to evaluate production-induced fault reactivation at Groningen Gas Field”. In: Proceedings - SPE Annual Technical Conference and Exhibition. Vol. 2015-Janua. 2015.doi: 10.2118/174942-ms. [56] C. Scholz. “Mechanics Of Faulting”. In: Annual Review of Earth and Planetary Sciences 17.1 (1989). issn: 00846597.doi: 10.1146/annurev.earth.17.1.309. [57] P. Segall. “Earthquakes triggered by fluid extraction”. In: Geology 17.10 (1989).issn: 00917613.doi: 10.1130/0091-7613(1989)017<0942:ETBFE>2.3.CO;2. [58] Paul Segall. Earthquake and volcano deformation. 2010.doi: 10.5860/choice.48-0287. [59] Hiroki Sone and Mark D. Zoback. “Mechanical properties of shale-gas reservoir rocks - Part 1: Static and dynamic elastic properties and anisotropy”. In: Geophysics 78.5 (2013).issn: 00168033. doi: 10.1190/GEO2013-0050.1. 129 [60] J. Speight. Shale Oil and Gas Production Processes. Cambridge, MA USA: Elsevier, 2020.doi: 10.1016/C2015-0-02189-X. [61] James Speight. Shale oil and gas production processes. 2019.doi: 10.1016/C2015-0-02189-X. [62] Yipeng Zhang, Mark Person, John Rupp, Kevin Ellett, Michael A. Celia, Carl W. Gable, Brenda Bowen, James Evans, Karl Bandilla, Peter Mozley, Thomas Dewers, and Thomas Elliot. “Hydrogeologic controls on induced seismicity in crystalline basement rocks due to fluid injection into basal reservoirs”. In: Groundwater 51.4 (2013).issn: 0017467X.doi: 10.1111/gwat.12071. [63] Xiaoxi Zhao and Birendra Jha. “Role of Well Operations and Multiphase Geomechanics in Controlling Fault Stability During CO2 Storage and Enhanced Oil Recovery”. In: Journal of Geophysical Research: Solid Earth 124.7 (2019).issn: 21699356.doi: 10.1029/2019JB017298. 130 Appendix ** Irving0830201 ** ============== INPUT/OUTPUT CONTROL ====================== ** 2018-01-29, 9:49:37 PM, abubshait RESULTS SIMULATOR STARS 201313 *interrupt *stop *title1 ’STARS Geomechanics Template #8’ *title2 ’Corner Point Grid Deformation with Displacement Vectors’ *INUNIT *FIELD *OUTUNIT *FIELD 131 **WRST 40 **SRFORMAT *SR3 *WPRN *GRID 1 *OUTPRN *GRID *VPOROS *PRES *SO *SG *SW *WPRN *SECTOR 1 *PRNTORIEN 2 0 *WPRN *ITER 1 *OUTPRN *ITER *NEWTON *WSRF *GRIDDEFORM 2 ** keyword for output grid deformation *WSRF *GRID 1 *OUTSRF GRID POISSON PRES STRESI STRESJ STRESK STRESSHIJ STRESSHIK STRESSHJK STRESEFF STRESINVF STRESINVS STRESSM STRNEPL VDISPL VMSTRESS VPOROS YOUNG STRAINVOL YLDSTATE BULKVOL SW SG *WSRF *SECTOR 1 132 ** ========= DEFINITION OF FUNDAMENTAL CARTESIAN GRID ========== *GRID *CORNER 248 120 129 *KDIR *DOWN *CORNERS *INCLUDE FlowMesh.txt *NULL *IJK *INCLUDE Null.txt *SECTOR ’S1’1:61 1:120 1:129 *SECTOR ’S2’62:92 1:120 1:129 *SECTOR ’S3’93:108 1:120 1:129 *SECTOR ’S4’109:139 1:120 1:129 *SECTOR ’S5’140:155 1:120 1:129 *SECTOR ’S6’156:171 1:120 1:129 *SECTOR ’S7’172:187 1:120 1:129 *SECTOR ’S8’188:218 1:120 1:129 *SECTOR ’S9’219:248 1:120 1:129 133 *POR *ALL *INCLUDE Porosity.txt *PERMI *ALL *INCLUDE Permeability.txt *MOD *INCLUDE FaultKi.txt *INCLUDE HPAE.txt PERMJ EQUALSI * 1 *MOD *INCLUDE FaultKj.txt *INCLUDE HPAE.txt PERMK EQUALSI * 0.1 *MOD *INCLUDE FaultKk.txt **INCLUDE HPAE.txt *END-GRID 134 **============= FLUID DEFINITIONS ============= *model 2 2 1 1 *compname ’Water’’CH4’ ** --------------- *cmm 18 16.043 *pcrit 0 667.1 *tcrit 0 -116.59 *massden 0 *cp 0 *VISCTABLE 62.4 0 100 0 *kv1 0 135 *kv4 0 *kv5 0 ** Reference conditions *prsr 14.65 *temr 62.4 *psurf 14.65 *tsurf 62.4 *xnacl 0 ** ============== ROCK-FLUID PROPERTIES ====================== *rockfluid *swt ** Water-oil relative permeabilities ** Sw Krw Krow Pcow ** ----- ----- ------ ------ 0.0500000000000000 0 1 0.287500000000000 0.250000000000000 0.750000000000000 0.525000000000000 0.500000000000000 0.500000000000000 136 0.762500000000000 0.750000000000000 0.250000000000000 1 1 0 *slt ** Liquid-gas relative permeabilities ** Sl Krg Krog Pcog ** ----- ------ ------ ------ 0.0500000000000000 0.800000000000000 0 0.287500000000000 0.600000000000000 0.250000000000000 0.525000000000000 0.400000000000000 0.500000000000000 0.762500000000000 0.200000000000000 0.750000000000000 1 0 1 ** ============== INITIAL CONDITIONS ====================== *INITIAL *VERTICAL *DEPTH_AVE *INITREGION 1 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 137 *DGOC 9.1864e+03 *DWOC 9.1864e+03 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 2 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 9.6238e+03 *DWOC 9.6238e+03 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 3 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 9.6238e+03 *DWOC 9.6238e+03 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL 138 *INITREGION 4 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 1.0061e+04 *DWOC 1.0061e+04 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 5 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 1.0645e+04 *DWOC 1.0645e+04 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 6 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 1.0645e+04 *DWOC 1.0645e+04 *TRANZONE 139 *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 7 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 1.0645e+04 *DWOC 1.0645e+04 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INITREGION 8 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 *DGOC 1.1228e+04 *DWOC 1.1228e+04 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL INITREGION 9 *REFDEPTH 1.2880e+04 *REFPRES 5.5772e+03 140 *DGOC 1.1228e+04 *DWOC 1.1228e+04 *TRANZONE *WOC_SW 1.0 *GASZONE *NOOIL *INTYPE *IJK 1:61 1:120 1:129 1 62:92 1:120 1:129 2 93:108 1:120 1:129 3 109:139 1:120 1:129 4 140:155 1:120 1:129 5 156:171 1:120 1:129 6 172:187 1:120 1:129 7 188:218 1:120 1:129 8 219:248 1:120 1:129 9 *SW *IJK *INCLUDE WaterSaturationProfile.txt 141 *TEMP *CON 100 *mfrac_gas ’CH4’*con 1 ** ================ NUMERICAL CONTROL =================== *NUMERICAL *ITERMAX 100 *NORTH 100 *DTMAX 30 *UNRELAX -1 *RANGECHECK *OFF *NORM PRESS 600 *CONVERGE *MAXRES 1e-05 *RANGECHECK *ON tform zt *isothermal *minpres 0.0015 142 ** ============== GEOMECHANIC MODEL ====================== *GEOMECH ** geomechanics main keyword *GEOM3D *SOLVERG *AIMSOL *GCOUPLING 2 ** Coupling option # 2 *GPTOLMUL 0.0015 ** Tolerance for pressure in coupling option 2 *ITERMAXG 100 *GEOGRID *GCORNER 88 40 43 GCORNERS *INCLUDE GeoMesh.txt GEOROCK 1 **Overburrden *ELASTMOD 2088543.42730397 *FRICANGLE 30 *POISSRATIO 0.2 *COHESION 2900.7547601444 *BIOTSCOEF 0.7 143 GEOROCK 2 **Marble Falls *ELASTMOD 8702264.2804332 *FRICANGLE 30 *POISSRATIO 0.2 *COHESION 2900.7547601444 *BIOTSCOEF 0.7 GEOROCK 3 **Barnett *ELASTMOD 5801509.5202888 *FRICANGLE 30 *POISSRATIO 0.23 *COHESION 2900.7547601444 *BIOTSCOEF 0.7 GEOROCK 4 **Ellenburger *ELASTMOD 8702264.2804332 *FRICANGLE 30 *POISSRATIO 0.2 *COHESION 2900.7547601444 *BIOTSCOEF 0.7 GEOROCK 5 **Basement *ELASTMOD 6236622.73431046 *FRICANGLE 30 144 *POISSRATIO 0.27 *COHESION 2900.7547601444 *BIOTSCOEF 0.7 GEOROCK 6 **Fault Upper *ELASTMOD 5801509.52 *FRICANGLE 31 *POISSRATIO 0.25 *COHESION 1000 *BIOTSCOEF 0.7 GEOROCK 7 **Fault Lower *ELASTMOD 5801509.52 *FRICANGLE 31 *POISSRATIO 0.25 *COHESION 0 *BIOTSCOEF 0.7 NITERGEO 50 *GEOTYPE *IJK *INCLUDE GeomechanicalProperties.txt *STRESS3D 293.7018 323.07198 959.4247 0.0 0.0 0.0 *STRESSGRAD3D -0.198934 -0.2188274 -0.6498503 0.0 0.0 0.0 145 *STRESS_RETURN *SUBSTEP *RUN ** ===================RECURRENT DATA ================ *DATE 2005 01 01 *DTWELL 0.05 *INCLUDE IrvingWells.txt *INCLUDE IrvingHWKQ.txt *TIME 5844 146

Abstract (if available)

Abstract Propagation of production- and injection-induced stresses from reservoir layers to basement faults is not well understood, especially in reservoirs with hydraulically segregated production and injection units situated across low-permeability structural or stratigraphic boundaries. Having a comprehensive understanding of how flow-induced stresses propagate and induce fault reactivation has a significant impact in optimizing oil and gas extraction, wastewater injection, and CO2 sequestration strategies, which can promote environmentally safe subsurface operations. Numerical modeling and simulation is the only tool available to forecast stress evolution at field scales and quantify the risk of fault reactivation under different geologic and fluid flow scenarios. Numerical modeling requires the solution of the coupled equations of fluid flow, rock deformation, and fault failure over a computational mesh discretizing the geologic domain. This is challenging in fields with different injection and production intervals, hundreds of hydraulically fractured wells with time-varying flow schedules, and multiple curved and connected faults. We present a framework to model spatiotemporal evolution of production- and injection-induced stress based on a coupled formulation of multiphase flow, poroelasticity, and fault failure over a dual mesh where the flow mesh discretizes the flow region at higher resolution than the resolution of the geomechanics mesh discretizing both flow and no-flow (overburden and basement) regions. We develop a corner point geometry meshing algorithm that can model reservoirs with multiple curved and connected faults. This expands the applicability of many industry-standard simulation software limited to structured meshes and reduces the need for building unstructured meshes at high computational costs. We compare two methods to model seismic events from fault failure events in our simulations, namely, moment tensor and potency-magnitude scaling methods and show that the potency-magnitude scaling method provides a promising alternative to the moment tensor method for event magnitude calculation from quasi-static simulations. We use the proposed framework to study 2008-2015 seismicity around Azle, Dallas, Fort Worth and Irving cities in Texas, USA by modeling the effects of gas and water production from Barnett Shale and wastewater injection into Ellenburger on the stress states of the shale reservoir, the aquifer, and the granitic basement. We analyze stress propagation onto faults to discover new mechanisms of faults reactivation in such complex fields.
We also show that flow-induced anisotropy and heterogeneity in the stress state can be analyzed in ultra-low permeability reservoirs undergoing hydraulic fracturing to optimize the trajectories of new and infill wells and their hydraulic fracturing jobs.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

$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

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

PDF

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

PDF

Dynamics of CO₂ leakage through faults

$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

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

Modeling and simulation of complex recovery processes

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

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

PDF

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

PDF

Detailed properties of seismic waveforms and earthquake sources in Southern California

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

Asset Metadata

Creator Bubshait, Abdulrahman Abdulaziz (author)

Core Title Reactivation of multiple faults in oilfields with injection and production

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Degree Conferral Date 2022-05

Publication Date 04/18/2022

Defense Date 06/15/2021

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

Tag Azle,Barnett,CO₂,computational,Earthquake,Ellenburger,failure,fault,gas,induced,injection,Irving,mesh,Modeling,multiphase,numerical,OAI-PMH Harvest,oil,permeability,poroelasticity,production,reactivation,reservoir,seismic spatiotemporal,simulation,stresses

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jha, Birendra (committee chair), Ben-Zion, Yehuda (committee member), Ershaghi, Iraj (committee member)

Creator Email abubshai@usc.edu,dhaimb@gmail.com

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

Unique identifier UC111006497

Document Type Dissertation

Format application/pdf (imt)

Rights Bubshait, Abdulrahman Abdulaziz

Type texts

Source 20220418-usctheses-batch-928 (batch), 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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

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

Repository Email cisadmin@lib.usc.edu