University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content CONTINUUM MODELING OF RESERVOIR PERMEABILITY ENHANCEMENT AND ROCK DEGRADATION DURING PRESSURIZED INJECTION by Mahshad Samnejad A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PETROLEUM ENGINEERING) May 2020 Copyright 2020 Mahshad Samnejad Dedication To my precious parents, Mitra and Mohammad, and my beloved brother, Arshan ii Abstract In this dissertation, we focus on hydromechanical evolution of changes in subsurface reservoir systems during pressurized injection. Our main contribution is to extend the applicability of the available Effective Continuum Model (ECM) methods, previously used for modeling enhanced geothermal systems, to the domains of hydraulic fracturing and waste water injection by addressing the major deficiencies of the current models including: • Multi-mechanismFractures: TheprevailingECMapproachforhydraulicfrac- turing simulation is based on Linear Elastic Fracture Mechanics (LEFM). Gen- erally, LEFM gives reasonable hydraulic fracturing results predictions for hard brittle rocks, but fails to predict the fracture propagation path and pressure pat- terns in formations which can undergo ductile failure, such as poorly consolidated and unconsolidated sands as well as ductile shales. There are very few previous attempts to model inelastic mechanical behavior during pressurized injection, all lacking a corresponding fluid flow response to plastic deformation. Lack of a proper ductile model leads to overestimation of the success of a hydraulic frac- turing job. We address this shortcoming by adopting an elastoplastic material model for the porous medium, which addresses the multi-mechanism fracture propagation in iii the rock along both the linear elastic and nonlinear inelastic deformation paths. Crack growth in the medium is then translated into the development of contin- uum damage, which is related to the loss of solid surface inside each volumetric element creating a better conduction capability for fluid flow, i.e. permeability enhancement. • Anisotropy Consideration: Previous ECM approaches only consider isotropic conditions, which is a reasonable assumption to make for geothermal systems made of the structurally-isotropic granite. Fractures formed in such systems often follow a sharp path away from the wellbore in the direction perpendicular to the minimum principal stress. It has been observed, however, that sedimentary rocks can exhibit more complex forms of fracture propagation patterns due to the pre- existing planes of weakness in their structure. Lack of accurate predictions on fracturepropagationpatternsresultsinsub-optimalwellplacementandtrajectory design. We model the tendency of fractures to grow in preferred directions by allocating a set of pre-existing planar flaws placed in different orientations to represent a form of secondary anisotropy inside an isotropic medium. Our model incorporates the effects of both stress asymmetry and structural anisotropy, demonstrating a better match with field microseimic measurements and lab observations. In summary, we exhibit the superiority of our hydraulic fracturing model by fitting it into the available coupled fluid flow and geomechanics platforms which allows to capture the phenomena occurring at a microscale in a macroscopic framework. Unlike other common approaches, our method does not require prior specification of the crack path, has no mesh dependency, and is capable of incorporating heterogeneity while iv maintaining a reasonable computational time. We also identify areas of opportunity for future research and propose a scope for extending beyond the current work. v Acknowledgments First and foremost, I thank God, the creator and sustainer of the Universe, who has helped me throughout my personal life and academic journey. I would like to express my deep appreciation to my Ph.D. advisor, Dr. Fred Amin- zadeh, who has motivated and supported me throughout my graduate studies. His broad insight in the oil and gas domain leads to constantly seeking novel solutions to the challenging problems of the petroleum industry in innovative practice-oriented ways. I thank him for his belief in my work, benevolent pieces of advice, and thoughtful encouragements at difficult times during my Ph.D. enrollment. I am intensely grateful to Professor Iraj Ershaghi for his continuous support over the course of my graduate studies at the University of Southern California (USC). I have benefited from his extensive domain knowledge, profound understanding, and invaluable guidelines at every situation in the past years. My gratitude to him goes beyond words. My sincere thanks also goes to Dr. Birendra Jha for his helpful advice as my co- advisor and helpful suggestions, which improved the quality of my research. I also had the privilege to have Dr. Charles Sammis and Dr. Aiichiro Nakano as my committee members, which enabled me to benefit from their constructive comments. I owe gratitude to the Petroleum Engineering program at the University of Southern California, which provided me with various means to nurture my career in petroleum vi engineering and data science. I also thank Baker Hughes for giving me an internship opportunity which enabled me to put into practice what I had learned at school in the two domains. I am thankful to all my peers and former teammates at USC, who contributed to my research with their help and collaboration. My special thanks go to Dr. Atefeh Jahandideh, Dr. Arman Khodabakhshnejad, Dr. Shahram Farhadi, Magdalene Ante, Robert L. Walker, Cyrus Ashayeri, and Qianru Qi. I am appreciative of the presence of my admirable friends, with whom I share years of unforgettable memories of ups and downs. My special thanks go to Mehrdad Gharib Shirangi, Marjan Sherafati, Roya Ermagan, Shahrzad Hadian, Mohammad Saeed Abr- ishami, Mohammad Reza Chitgarha, Navid Azizan Ruhi, Sina Hasheminassab, Kiana Roshan Zamir, Samira Aghayee, Shahab Helmi, Sona Sehat, Zahra Abrishami, and many others not mentioned here. I would also like to sincerely thank my aunt and uncle, Mandana and Hamid Habibelahian, for their incredible support, kindness, and hospitality which always made me feel like home at their residence. My deepest gratefulness goes to my family–my precious parents, Mitra and Moham- mad, and my beloved brother, Arshan–who are my role models of humaneness and strength. If it was not for their unconditional love, I would have not been able to come this far. I appreciate each and every bit of inspiration, encouragement, motivation, guidance, and assistance from my brother, without which I would have never been able to endure at hard times. I am deeply indebted to my family for life. Lastly, I feel obliged to thank Iran, my splendid homeland, Iran, for being the unwavering foundation upon which I began all this journey in life. vii Publications • Samnejad, M, Aminzadeh, F, and Jha, B, 2018, A Novel Approach for Studying Hydraulic Fracturing Success Factors beyond Brittleness Indices, 52nd US Rock Mechanics/Geomechanics Symposium, Seattle, Washington, D.C. (Samnejad et al., 2018) • Samnejad, M, Walker, R, and Aminzadeh, F, 2018, Understanding Injection- induced Seismicity Effects of Fault Damage Zones: Beyond Poroelastic Models, AGU Fall Meeting, Washington, D.C. (Samnejad et al., 2018) • Walker,R,Samnejad,M,andAminzadeh,F,2018,ConsiderationofStressInduced Permeability Changes Within a Poroelastic Model Framework for Induced and Triggered Seismicity, AGU Fall Meeting, Washington, D.C. (Walker et al., 2018) • Samnejad,M,Aminzadeh,F,andJha,B,2017,SimulationofHydraulicFracturing- Induced Permeability Stimulation Using Coupled Flow and Continuum Damage Mechanics,SPEAnnualTechnicalConferenceandExhibition,SanAntonio,Texas (Samnejad et al., 2017b) • Samnejad, M, Aminzadeh, F, and Bubshait, A, 2017, Acoustoelastic Estimation of the In Situ Stresses from Sonic Logs for a Carbonate Reservoir, SPE Western Regional Meeting, Bakersfield, California (Samnejad et al., 2017a) viii • Samnejad, M, Aminzadeh, F, and An, Y, 2014, The Investigation Of An Inte- grative Approach For Fracture Density Inversion Using AzAVO Analysis, Pacific Section AAPG, SPE and SEPM Joint Technical Conference, Bakersfield, Califor- nia (Samnejad et al., 2014) ix Contents Dedication ii Abstract iii Acknowledgments vi Publications viii List of Tables xiii List of Figures xiv 1 Introduction 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 National Economy . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.3 Public Welfare . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.2 Problem Statement & Challenges . . . . . . . . . . . . . . . . . 10 1.2.3 Our Approach & Advantages . . . . . . . . . . . . . . . . . . . 11 1.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 Manuscript Organization . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2 Subsurface Multiphysics 18 2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.1 Mechanical Equilibrium . . . . . . . . . . . . . . . . . . . . . . 19 2.1.2 Fluid Flow through Porous Media . . . . . . . . . . . . . . . . . 19 2.1.3 Fracturing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Material Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 x 2.2.1 Mechanical Deformation . . . . . . . . . . . . . . . . . . . . . . 35 2.2.2 Continuum Damage Evolution . . . . . . . . . . . . . . . . . . . 48 2.2.3 Mechanical Degradation . . . . . . . . . . . . . . . . . . . . . . 50 2.2.4 Flow Enhancement . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3 ModelingWorkflowforPermeabilityEnhancementandMaterialDegra- dation during Pressurized Injection 59 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2 Implementation of an Integrative Workflow . . . . . . . . . . . . . . . . 60 3.3 Damage Evolution Model Validation . . . . . . . . . . . . . . . . . . . 65 3.4 Mechanical Degradation . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.4.1 Horizontal Cracks . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4.2 Vertical Cracks . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.5 Permeability Enhancement . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4 ReservoirSimulationofInjection-InducedPermeabilityEnhancement 79 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2 Mathematical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2.1 Damage Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2.2 Mechanical Model of the Damaged Medium . . . . . . . . . . . 82 4.2.3 Fluid flow Model inside Damaged Medium . . . . . . . . . . . . 86 4.3 Numerical Experiment at the RVE Level . . . . . . . . . . . . . . . . . 87 4.3.1 Example 1: Uniaxial Tensile Test . . . . . . . . . . . . . . . . . 88 4.3.2 Example 2: 2D Injection-Induced Stimulation . . . . . . . . . . 90 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5 Identifying Dominant Success Factors in Hydraulic Fracturing 99 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Conventional Brittleness Indices . . . . . . . . . . . . . . . . . . . . . . 100 5.3 Mathematical Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.3.1 Governing Equations of Reservoir Multi-physics Processes . . . 101 5.3.2 Fluid Conservation . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.3.3 Mechanical Equilibrium . . . . . . . . . . . . . . . . . . . . . . 102 5.3.4 Continuum Damage Mechanics Model for Element Behavior . . 102 5.4 Multi-Mechanism Crack Growth Modeling . . . . . . . . . . . . . . . . 105 5.4.1 Pre-existing Cracks Orientation . . . . . . . . . . . . . . . . . . 107 5.4.2 In-situ Stress Magnitude . . . . . . . . . . . . . . . . . . . . . . 108 5.4.3 In-situ Stress Anisotropy . . . . . . . . . . . . . . . . . . . . . . 109 5.4.4 Degree of Ductility . . . . . . . . . . . . . . . . . . . . . . . . . 110 xi 5.4.5 Yield Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.5 Operational Implications of Impact Factors . . . . . . . . . . . . . . . . 111 5.5.1 Pre-existing Cracks Orientation . . . . . . . . . . . . . . . . . . 112 5.5.2 In-situ stress magnitude . . . . . . . . . . . . . . . . . . . . . . 112 5.5.3 In-situ Stress Anisotropy . . . . . . . . . . . . . . . . . . . . . . 113 5.5.4 Degree of Ductility . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.5.5 Yield Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6 Fault Damage Zones Effects on In-situ Geomechanical Stresses 118 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.2 Overview of Seismicity . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.2.1 Mechanics of Fault Rupture . . . . . . . . . . . . . . . . . . . . 120 6.2.2 Mechanisms of Induced Seismicity . . . . . . . . . . . . . . . . . 122 6.2.3 The Role of Fault Damage Zones . . . . . . . . . . . . . . . . . 124 6.3 Governing Physics of Poroelascitiy . . . . . . . . . . . . . . . . . . . . 125 6.3.1 Geomechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6.3.2 Fluid Flow in Porous Media . . . . . . . . . . . . . . . . . . . . 126 6.4 Coupled Poroelastic Model of Fluid Injection . . . . . . . . . . . . . . . 127 6.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.5.1 The effect of hydromechanical irregularity in fault damage zones 132 6.5.2 The effect of the geometry of fault damage zones . . . . . . . . 134 6.5.3 The effect of fault transmissibility . . . . . . . . . . . . . . . . . 135 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7 Summary, Conclusions, and Future Work 140 Reference List 144 xii List of Tables 3.1 Measured Parameters of Uniaxial Tensile Test. . . . . . . . . . . . . . . 66 3.2 Calibration Parameters for Uniaxial Tensile Test. . . . . . . . . . . . . 67 3.3 Parameters of Triaxial Compression Test. . . . . . . . . . . . . . . . . . 73 4.1 Uniaxial Tensile Test Parameters. . . . . . . . . . . . . . . . . . . . . . 88 5.1 Brittle crack growth parameters. . . . . . . . . . . . . . . . . . . . . . . 106 5.2 Ductile crack growth parameters. . . . . . . . . . . . . . . . . . . . . . 106 6.1 Simulation Domain Parameters. . . . . . . . . . . . . . . . . . . . . . . 128 6.2 Fluid Properties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 xiii List of Figures 1.1 U.S. primary energy sources – data from the Monthly Energy Review of U. S. Energy Information Administration (2019). . . . . . . . . . . . . 1 1.2 Schematicexamplesofsubsurfaceactivitiespertainingtofluidflowthrough fractured rock is of importance: a) hydraulic fracturing b) waste disposal sites c) enhanced geothermal systems. . . . . . . . . . . . . . . . . . . . 2 1.3 2014 estimates of the economic impact of U.S. unconventionals including direct extraction, transport, and refining of unconventional oil and gas, as well as activities that support this production, such as oil field services and local services, data after Porter et al. (2016). . . . . . . . . . . . . 3 1.4 Hydraulic fracturing market is estimated to hit $66 billion by 2020, data after Grand View Research, Inc. (2016). . . . . . . . . . . . . . . . . . 4 1.5 Increasing seismicity rate from 2009, data after U.S. Geological Survey (USGS) (Rubinstein, 2015). . . . . . . . . . . . . . . . . . . . . . . . . 6 1.6 The coupled multiphysics involved in hydraulic fracturing include fluid flow through media, geomechanical deformation, and rock fracturing. . 8 1.7 Challenges associated with modeling the coupled process of fluid flow, geomechanical deformation, and rock fracturing. . . . . . . . . . . . . . 10 xiv 1.8 Advantages of Effective Continuum Models in representing the coupled process of fluid flow, geomechanical deformation, and rock fracturing dynamics of fluid flow and rock deformation in presence of fractures is central to designing and predicting hydraulic stimulation in these sub- surface applications (Shirangi et al., 2019; Jahandideh & Jafarpour, 2016). 12 1.9 Planar fractures are oversimplified and mesh dependent. a) from Was- antha & Konietzky (2017) b) from matthewpais.com . . . . . . . . . . 13 1.10 DFN models are also mesh dependent and require prior specification of the failure path. a) from Maxwell et al. (2016) b) from mve.com . . . . 14 1.11 ECM models have an advantageous ease of implementation and can be used to incorporate flow-related phenomena as well. a) from tu- freiberg.de b) from northwestern.edu . . . . . . . . . . . . . . . . . . . 15 2.1 The coupled multiphysics involved in hydraulic fracturing include fluid flow through media, geomechanical deformation, and rock fracturing. . 18 2.2 DamageinducedbyaxialloadF onacylindricalvolumeelementinterms of the fractional area δS D /δS occupied by micro-mechanical defects δS denoted by the normal vector n after Lemaitre (1996). . . . . . . . . . 30 2.3 Principal elements of the damage tensor in 3D (a) damaged configura- tion (b) strain-equivalent configuration denoted by tilde after Lemaitre (1996). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4 Microdefects in a 2D REV idealized as microcracks of a lower dimension (1D) (a) an individual crack characterized by its half-length,a, and angle with horizon, ψ (b) multiple crack families with various lengths and orientations distinguished by color. . . . . . . . . . . . . . . . . . . . . 34 2.5 The stress-strain curve for a linear elastic material under uniaxial loading. 36 xv 2.6 Experimental stress-strain curves obtained for various rock types demon- strating inelastic behavior: a) Jurassic Shale after Papanastasiou (1997) b)BereaSandstoneafterRutter&Glover(2012)c)HollingtonSandstone after Rutter & Glover (2012). . . . . . . . . . . . . . . . . . . . . . . . 39 2.7 Stress-strain curves plotted for two different types of material nonlinear- ity: a) nonlinear elasticity b) and elasto-plasticity. . . . . . . . . . . . . 39 2.8 Yield surface for a 2D stress space. . . . . . . . . . . . . . . . . . . . . 40 2.9 Stress-strain diagrams under: a) Perfect plasticity b) Hardening plasticity. 41 2.10 a) Isotropic hardening b) Kinematic hardening. . . . . . . . . . . . . . 42 2.11 Mohr-Coulomb criterion in principal stress space. . . . . . . . . . . . . 43 2.12 von Mises criterion in principal stress space. . . . . . . . . . . . . . . . 44 2.13 Drucker-Prager criterion in principal stress space. . . . . . . . . . . . . 45 2.14 Return mapping of the state of stress using an elastic predictor-plastic corrector algorithm in 2D. . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.15 Permeability enhancement attenuation in the inelastic region, plotted after Peach & Spiers (1996). . . . . . . . . . . . . . . . . . . . . . . . . 53 2.16 Permeability and porosity grow in different ways under strain. Porosity increases almost linearly with strain. Permeability enhancement, how- ever, attenuates on the onset of material plasticity. Porosity-based and empirical permeability models cannot model this transition of material flow behavior. Figure is plotted after Popp et al. (2001). . . . . . . . . 53 2.17 Directionality of permeability and damage: (a) Fracture and matrix are in series and the matrix permeability dominates flow (min permeability). b) Fracture and matrix are in parallel and the fracture permeability dominates flow (max permeability). . . . . . . . . . . . . . . . . . . . . 54 2.18 Schematics of the local stresses and applied stresses after Costin (1983). 57 xvi 3.1 An integrative workflow presented for the simulation of multi-mechanism hydraulic fracturing processes. . . . . . . . . . . . . . . . . . . . . . . . 61 3.2 Schematics of the pre-damaged samples with embedded fissures used in uniaxial tension tests. In the work by Jalbout (2006), ψ is set at 90 o for vertical cracks and 0 o for horizontal cracks. . . . . . . . . . . . . . . . . 67 3.3 Conceptual schematics of crack growth inside an REV during a uniax- ial tensile experiment. The incremental tensile load is applied in the z direction. The induced damage evolves depending upon the crack plane direction: a) Horizontal cracks induce damage in the z direction. b) Vertical cracks aligned perpendicular to they axis induce damage in the same direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.4 Stress-strain diagram for a damaging REV with horizontal fissures in response to increasing load compared to the linear behavior of intact medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.5 Horizontal crack growth inside REV in response to increasing load. Crit- ical crack length causes material failure prior to reaching intact rock strength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.6 Damage accumulation inside the sample reduces the elements of the stiff- ness tensor, changing the material’s mechanical behavior and causing deviation from linear elasticity. . . . . . . . . . . . . . . . . . . . . . . 71 3.7 Young’s Modulus drops due to the loss of effective solid surface inside each element. See Equation 2.61. Data after Jalbout (2006). . . . . . . 72 3.8 Damage accumulation caused by the growth of horizontal crack inside the REV starting from a negligible initial state to the critical value. . . 73 3.9 Stress-strain relationship corresponding to the growth of vertical cracks in an REV causing deviation from the intact linear response. . . . . . . 74 xvii 3.10 The growth of vertical cracks has a slow pace compared to that of hor- izontal fissures under incrementally increasing load, except for the late regions close to material failure. This is due to the fact that deviatoric stress components are the only driving mechanism for such deficiencies. See Equation 2.93. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.11 A deteriorating stiffness tensor component creating a non-linear response to axial load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.12 The loss of effective solid surface inside each element is translated into a decrease in Young’s modulus as in Equation 2.61. Data after Jalbout (2006). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.13 the damage component in the x direction evolves from a small value in the beginning to an ultimate critical value. . . . . . . . . . . . . . . . . 76 3.14 Porosity growth plotted against the axial strain ratio. Porosity grows linearly throughout the process, as it is only a function of the volumetric strain. See Equation 2.100. Data after Popp et al. (2001). . . . . . . . 77 3.15 Multi-mechanism crack growth causes an increase in permeability. This enhancement significantly slows down when the material yields at the end of brittle crack growth region. Plastic deformation only allows small permeability growth. Data after Popp et al. (2001). . . . . . . . . . . . 78 4.1 SchematicsofRVEwithmicrocracks(a)elementcontainingasinglecrack (b) element with two families of cracks. . . . . . . . . . . . . . . . . . . 81 4.2 Stress concentration around a hole in a 2D stressed plate. . . . . . . . . 82 xviii 4.3 The process of damage growth inside an RVE due to increasing load from left to right:(a) Concentrated load is lower than yield stress; there- fore, no crack growth happens. (b) Local plastic deformation leads to crack growth. (c) Concentrated stress exceeds material ultimate strength resulting in failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.4 Element constitutive law capturing the three stages of mechanical dam- age and failure: Elastic region maintaining the original damage, plastic hardening allowing microcrack extension, and formation of microfracture with post-failure material softening. In this figure, σ u is the ultimate strength of the material, σ y is the material yield stress, and σ r denotes the residual strength of the rock after failure. . . . . . . . . . . . . . . . 85 4.5 Permeability evolution with damage, whereD 0 is the onset of permeabil- ity enhancement, and Dc is when a macro-fracture forms. . . . . . . . . 87 4.6 Stress-strain response. . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.7 Crack growth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.8 Damage evolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.9 Permeability enhancement. . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.10 Model mesh geometry and boundary conditions. . . . . . . . . . . . . . 93 4.11 Reservoir state after an hour of injection. a) Pressure b) von Mises stress. 94 4.12 Effective plastic strain after an hour of injection around the wellbore. . 95 4.13 Microcrackgrowthafteranhourofinjectionaroundthewellbore. Cracks are preferentially oriented along y axis because of the applied compres- sion on the top boundary. . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.14 Damage evolution after an hour of injection around the wellbore. . . . . 96 4.15 Permeability enhancement after an hour of injection around the wellbore. 96 xix 4.16 Anisotropy in permeability enhancement. The two curves correspond to two elements in the mesh located at the co-ordinates shown in the legend. 97 4.17 Permeabilityenhancementinthexdirection,(a)cutline(b)permeability change with distance away from the wellbore. . . . . . . . . . . . . . . 97 4.18 Permeabilityenhancementintheydirection,(a)cutline(b)permeability change with distance away from the wellbore. . . . . . . . . . . . . . . 98 4.19 Comparisonofa)CrackLength, b)Permeability, andc)Damageresulted during pressurization of ductile (top row) and brittle (bottom row) mate- rial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.1 Schematic view of a volumetric element. . . . . . . . . . . . . . . . . . 103 5.2 Variation of scalar function f(r) with relative crack length r. . . . . . . 104 5.3 Pre-existing crack orientation effect on brittle crack growth. . . . . . . 107 5.4 Maximum principal stress effect on brittle crack growth. . . . . . . . . 109 5.5 In-situ stress anisotropy effect on brittle crack growth. . . . . . . . . . 110 5.6 Ductility effect on plastic crack growth. . . . . . . . . . . . . . . . . . . 111 5.7 Effect of the threshold pressure on plastic crack growth. . . . . . . . . . 112 5.8 Effect of the direction of pre-existing inclusions on the required injection pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.9 Effect of the magnitude of in-situ stress on the required injection pressure.114 5.10 Effect of stress anisotropy on the required injection pressure. Lower values on the x-axis indicates higher anisotropy. . . . . . . . . . . . . . 115 5.11 Effect of the rock ductility parameter on the required injection pressure. 116 5.12 Effect of the threshold pressure on the required injection pressure. . . . 117 xx 6.1 Schematic of 2-D fault planes (hatched surfaced) inside 3-D host rock blocks. It is often convenient to break the traction vector, t, into its normal,σ n , and shear,τ, components for a fault with a normal direction indicated by n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 6.2 Schematics of induced earthquakes mechanisms caused by fluid injection (a) Pore pressure elevation in waste disposal sites (b) Rise in the normal stress caused by a near-by energy production site. . . . . . . . . . . . . 123 6.3 Mohr Circle diagram of induced earthquakes mechanisms caused by fluid injection shows how the stress state perturbations caused by fluid injec- tion can bring an initially-stable fault system (solid semicircle) to a state of failure (dashed semicircle). (a) Pore pressure elevation in waste dis- posal sites (b) Rise in the normal stress caused by a near-by energy production site. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.4 Schematics of fault planes surrounded by highly permeable damage zones.124 6.5 Model Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.6 The evolution of a) Pore Pressure, b) Effective Normal Stress on the Fault Plane and c) Coulomb Stress after 5 days of pressurized injection in proximity of a sealing fault system with damage zone. . . . . . . . . 129 6.7 Induced seismicity potential for sealing faults with and without damage zones: a) The effective normal stress plotted against depth indicates more severe fault destabilization in presence of damage zones. b) Fault stress path is accelerated towards the failure line when damage zones exist. c) Damage zones cause a quicker drop of the effective normal stress, introducing safety concerns. . . . . . . . . . . . . . . . . . . . . 130 xxi 6.8 Induced seismicity potential for nonsealing faults with and without dam- age zones: a) More pronounced fault destabilization with fault damage zones b) Acceleration of fault stress path towards the failure line when damage zones exist c) Unsafe drop of the effective normal stress in pres- ence of fault damage zones. . . . . . . . . . . . . . . . . . . . . . . . . 131 6.9 Induced seismicity potential for sealing faults in presence of low and highly damaged zones: a) The effective normal stress plotted against depth indicates more severe fault destabilization in presence of highly damaged zones. b) The shifting of the stress path towards the failure line accelerates with the existence of highly damaged zones. c) Rapid drop of the effective normal stress, introducing safety concerns around highly damaged fault zones. . . . . . . . . . . . . . . . . . . . . . . . . 133 6.10 Induced seismicity potential for nonsealing faults in presence of low and highly damaged zones: a) More pronounced fault destabilization with highly damaged zones even when the fault is nonsealing b) Acceleration of fault stress path towards the failure line when damages are high c) Quick and unsafe drop of the effective normal stress in presence of highly damaged zones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.11 Induced seismicity potential for sealing faults in presence of small and largedamagedzones: a)Theeffectivenormalstressplottedagainstdepth indicates more severe fault destabilization in presence of larger damaged zones. b) The shifting of the stress path towards the failure line accel- erates with the existence of larger damaged zones. c) Rapid drop of the effective normal stress, introducing safety concerns around large dam- aged fault zones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 xxii 6.12 Induced seismicity potential for nonsealing faults in presence of small and large damaged zones: a) More pronounced fault destabilization with larger damaged zones b) Acceleration of fault stress path towards the failure line when damages are large c) Quick and unsafe drop of the effective normal stress in presence of large damaged zones. . . . . . . . 136 6.13 Induced seismicity potential for sealing and nonsealing faults without damaged zones: a) The effective normal stress plotted against depth indicates very similar stress drops for both systems when damage zones are not present. b) The shifting of the stress path towards the failure line shows no observable difference between sealing and nonsealing faults. c) The effective normal stress drops in a somewhat similar fashion for both sealing and nonsealing. . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.14 Induced seismicity potential for sealing and nonsealing faults in presence of damage zones: a) More severe fault destabilization for sealing faults b) More acceleration of the shifting of the stress path towards the failure line for sealing faults c) Quick and unsafe drop of the effective normal stress when more pressurization created by sealing faults. . . . . . . . . 138 7.1 Summary of contributions and comparison with previous work. . . . . . 141 7.2 Benefits of AI-assisted physics modeling. . . . . . . . . . . . . . . . . . 143 xxiii Chapter 1 Introduction 1.1 Motivation 1.1.1 National Economy As a fuel to power the engine of economy, energy is of utmost importance to the national security of any given country (Nandakumar, 2007). Reducing a country’s dependence on the imported energy helps mitigate economic losses caused by foreign oil supply disruptions (Brown & Huntington, 2010). As shown in Figure 1.1, Fossil fuels still dominate the energy market by providing more than 80% of the nation’s total energy demand. Another roughly 10% is provided by nuclear resources, which- similar to hydrocarbon production and geothermal sites-rely on operations including fluid injection into subterranean geologic formations, see Figure 1.2. Figure 1.1: U.S. primary energy sources – data from the Monthly Energy Review of U. S. Energy Information Administration (2019). 1 Figure 1.2: Schematic examples of subsurface activities pertaining to fluid flow through fractured rock is of importance: a) hydraulic fracturing b) waste disposal sites c) enhanced geothermal systems. United States has benefited economically from the boom in horizontal drilling and hydraulic fracturing technology, which rendered the development of domestic unconven- tional energy resources possible. U.S. production of liquid fuels surpassed the Middle East in 2013 (Yo & Neff, 14 July 2014), adding 169,000 jobs between 2010 and 2012 (Brown & Yücel, 2013). According to a report published by the United states Depart- ment of Labor (Pickenpaugh & Adder, 2018), Pennsylvania led a shale-related employ- ment growth of 121 percent from 2007 to 2016, while Ohio had a rise of 62 percent in the same years, compared with West Virginia’s 41-percent growth. Figure 1.3 shows the results of a report published by Harvard Business School highlighting the estimated economic benefits obtained from unconventional resources in 2014 (Porter et al., 2016). 2 2014 Estimate of the Economic Impacts of Unconventionals Figure 1.3: 2014 estimates of the economic impact of U.S. unconventionals including direct extraction, transport, and refining of unconventional oil and gas, as well as activities that support this production, such as oil field services and local services, data after Porter et al. (2016). It is expected for natural gas to witness a considerable growth in the next decade owing to wide applications across industries. It has an environment-friendly nature, making it preferable over other fuels. According to a report by Grand View Research, Inc. (2016), tight gas accounted for over 15% of the global hydraulic fracturing market share in 2015. Horizontal well fracturing is the main technology allowing access to these resources. Figure 1.4 shows the increasing global market revenue obtained by the old and new techniques of hydraulic fracturing well stimulation, namely, plug & perf and sliding sleeve respectively. Plug and perforation hydraulic fracturing market size witnessed at a rising value of over $15 billion in 2015. 3 Global Hydraulic Fracturing Market Revenue Figure 1.4: Hydraulic fracturing market is estimated to hit $66 billion by 2020, data after Grand View Research, Inc. (2016). The increase in demand for oil and natural gas drives the need for fracturing activ- ities and investments. Acumen Research and Consulting (2019) has reported that the hydraulic fracturing market size is estimated to hit $83 billion by 2026. According to Grand View Research, Inc. (2016) and given the abundant shale reserves exploration in North America, this region currently has the largest hydraulic fracturing market value. China is anticipated to have an exponential growth in hydraulic fracturing activities, leading to more demand of hydraulic fracturing analytical technology, equipment, and service owed in majority by American companies. 4 1.1.2 Summary Hydraulic fracturing has been growing in the past decade then as a pivotal tech- nology in various underground fluid extraction activities, such as unconventional oil and gas reservoirs and enhanced geothermal systems. In these operations, formation permeability is enhanced when fractures are created and/or stimulated by injecting a highly-pressurized fluid, which is critical to maintaining commercial fluid production. In all these activities, the fluid transport phenomenon both affects and is affected by rock deformation and fracturing processes. The success of investment decisions pertaining to the exploitation of unconventional resources depends strongly on the reliability of the models making predictions of post- stimulationperformance. However, duetoalackofmodelsbasedonaccurateknowledge of the reservoir and rigorous understanding of the governing physics, there is a tech- nology gap between the current models of stimulation and the field observations in the Exploration and Production (E & P) industry, causing the bypassing of significant oil reserves. 1.1.3 Public Welfare In certain other subterranean applications, such as underground waste water dis- posal, CO 2 sequestration, and waterflooding, injection-induced fractures are inhibited to avoid potential leakage, groundwater contamination, fault pressurization leading to induced seismicity. Oil and gas extraction, hydraulic fracturing of geothermal settings, nuclear waste disposal, and gas storage are among the activities which, on the one hand, interlink with the nation’s energy security, but are postulated to potentially impose seismicity risk on the other hand (Foulger et al., 2018) and possibly cause 5 other environmental impacts, such as groundwater contamination (Jabbari et al., 2017; Aminzadeh, 2018, 2019,?). Although it has been found that only a small portion of energy development activi- ties in the United States has contributed to the triggering of induced seismicity at levels noticeable to human beings (National Research Council, 2013), certain significant seis- mic events at a number of locations, such as the 1952 magnitude 5.7 earthquake in Oklahoma, may have been induced by waste water injection (Hough & Page, 2015; Langenbruch & Zoback, 2016; Hosseini et al., 2018). Some of the induced seismicity created in California oil fields are highlighted in Tiwari et al. (2014); Goebel et al. (2016); Hough et al. (2017). One may also find it crucial to look into the potential induced earthquake hazards associated with subsurface activities, which has been an area of public concern especially after observing the recent rise in seismic activity, see Figure 1.5. Figure 1.5: Increasing seismicity rate from 2009, data after U.S. Geological Survey (USGS) (Rubinstein, 2015). 6 1.1.4 Summary Considering the dramatic casualties caused by a large scale earthquake disaster as well as the damage to buildings and infrastructure-which can also result in severe environmentalconsequencesincludingair, water, andsoilpollution(Nolaetal.,2013)-it is necessary to understand the physical mechanisms involved in seismic events in order to take preventative measures. Despite this urge, scientists have failed to accurately predict and warn the public and the related agencies of the occurrence of induced earthquakes due to uncertainty about the subsurface rock systems and lack of predictive models developed for prescriptive purposes (National Research Council, 2013; Shirangi et al., 2019). Please see the University of Southern California’s Induced Seismicity Consortium (ISC.usc.edu) for work done on related topics since 2013. 1.1.5 Conclusion With emerging the computational capabilities, understanding and modeling the dynamics of fluid flow and geomechanics in presence of fracture initiation and propaga- tion have become possible. These phenomena, schematically shown in Figure 1.6, are central to the profit prediction and safe design of a variety of subsurface exploitation applications and have raised vast research interest in the development of physically realistic, mathematically rigorous, and numerically efficient models of fluid flow and geomechanical deformation in fractured porous media. In this work, aim to address this complex problem by introducing our representative model. 7 Figure 1.6: The coupled multiphysics involved in hydraulic fracturing include fluid flow through media, geomechanical deformation, and rock fracturing. 1.2 Research Objectives 1.2.1 Overview The energy industry is increasingly demanding appropriate computational models of rock fracturing processes for proper design, implementation, and monitoring of var- ious subsurface exploitation activities such as hydraulic fracturing, well stimulation, enhanced oil recovery, enhanced geothermal energy extraction, CO 2 sequestration, and nuclear waste disposal. The success of these operations highly relies on the accuracy of the models representing the multiphysics nature of such processes in geologic forma- tions. In this dissertation, we review the common methods that are used for describing reservoir fracturing processes, highlight their limitations, and present an integrative 8 simulation workflow to provide with a more representative model of such complex mul- tiphysical phenomena. The outcome of this work is essentially used for commercial fluid recovery predictions as well as the safe design of optimal development plans. Among the variety of approaches used in the literature for addressing the prob- lem of modeling fracturing processes in porous media, the use of Effective Continuum Models (ECM) has a number of advantages including ease of implementation with no prior specification of the crack path and the ability to include the variations of sub- surface rock properties from both mechanics and fluid flow points of view. ECMs are based upon Continuum Damage Mechanics principles (CDM) and have been adopted in subterranean applications with the initial purpose of representing geothermal set- tings. Such workflows are later integrated into computational frameworks bringing in the ability to perform dynamic studies of the fracture growth and the resulting flow enhancement. Despite the vast applicability of ECMs to geothermal systems, little to no promising result has been obtained in the past for the case of petroleum reservoirs. This is mainly due to the fact that hydrocarbon-bearing formations are inherently more complex than granite geothermal rocks, making their physics models more complicated to build. In this work, we introduce a rigorous effective continuum model of subsurface fracturing phenomena designed to address the specific complexities associated with petroleum reservoirs. Our approach, originally developed for hydraulic fracturing sim- ulation purpose, also has extensive potential application to the case of waste injection scenarios serving a rather different aim of fracturing prevention and control. 9 1.2.2 Problem Statement & Challenges Addressing certain practical challenges associated with hydraulic fracturing opera- tions, including country-specific restrictions related to potential environmental effects including noise and visuals impact, induced seismic events, land surface disturbance, methaneemissionandwatercontamination, amongothers. andhighoperationalexpen- ditures involved in the development of unconventional resources are beyond the scope of this work. Besides the above, modeling multiple interacting physical phenomena–fluid flow, rock deformation, and fracturing–calls for proper integration of different computational simulation platforms making the task of hydraulic fracturing simulation very challeng- ing. These phenomena occur at different scales inside a geometrically complex envi- ronment (Shirangi, 2016) where only sparse and mostly indirect measurements can be made for validation purposes, see Figure 1.7. Figure 1.7: Challenges associated with modeling the coupled process of fluid flow, geomechanical deformation, and rock fracturing. 10 This dissertation will address these technology challenges by developing a represen- tative modeling and simulation framework for capturing the multiphysical dynamics of the reservoir behavior in response to hydraulic fracturing. 1.2.3 Our Approach & Advantages Using a sequentially coupled numerical scheme of fluid flow, rock deformation and poromechanical damage we have developed a reservoir engineering simulation workflow. Ithelpsevaluatetheevolutionoffracturedrockpropertiesandstatevariablesasaresult of hydraulic stimulation. We regards the stimulated rock as a mechanically damaged continuum and model the evolution of rock properties by means of a synthetically modeled equivalent configuration. Our approach is essentially in contrast to other hydraulic fracturing simulation methodologies, in which the rock and fractures are considered two separate entities, and the focus is more attention is paid towards the characteristics of pre-defined frac- tures than the overall properties of the target rock. By contrast, we take a reservoir engineering standpoint and focus on the changes of flow and mechanical properties of the rock, namely, permeability and stiffness, as a result of hydraulic fracturing. We incorporate the principles of damage mechanics into the multiple governing physics of macroscale fracture formation processes. This allows to investigate the dynamics of pressure distribution and stress state during hydraulic fracturing activities. We revisit material constitutive relations with microcrack growth, and use experimental observa- tions to provide insight on flow behavior in presence of damage. We then implement our model inside a coupled poromechanical scheme for computational simulation of a subsurface fluid injection scenario so as to demonstrate the capability of our approach for hydraulic fracturing progressive results prediction. Figure 1.8 highlights the benefits 11 gained by using an Effective Continuum Modeling approach for modeling such complex processes. Figure 1.8: Advantages of Effective Continuum Models in representing the coupled process of fluid flow, geomechanical deformation, and rock fracturing dynamics of fluid flow and rock deformation in presence of fractures is central to designing and predicting hydraulicstimulationinthesesubsurfaceapplications(Shirangietal.,2019;Jahandideh & Jafarpour, 2016). 1.3 Literature Review Hydraulic stimulation has been growing as an important technology in the subsur- face energy applications such as shale gas production (Curtis, 2002) and the devel- opment of enhanced geothermal systems (Zimmermann & Reinicke, 2010). In cer- tain other applications, such as underground waste water disposal and waterflooding, 12 injection-induced fractures are inhibited to avoid leakage (Pruess, 2005), groundwater contamination, fault pressurization and seismicity (Bao & Eaton, 2016), and bypass- ing of oil reserves (Gadde et al., 2001). In all these activities, fluid flow processes both affect and are affected by rock deformation and fracture mechanics. Understanding and modeling the dynamics of fluid flow and rock deformation in presence of fractures is cen- tral to designing and predicting hydraulic stimulation in these subsurface applications (Jahandideh & Jafarpour, 2016; Shirangi et al., 2019). A variety of methods have addressed the problem of predicting fluid flow and geome- chanical response of hydraulically fractured reservoirs. For instance, classical models of planar fractures growing away from the wellbore have been widely used in the industry to obtain fracture dimensions under idealized conditions (Geertsma et al., 1969; Nord- gren & others, 1972; Rice & Cleary, 1976), as shown in Figure 1.9. Such traditional methods of bi-wing hydraulic fracture growth are applicable only under oversimplified conditions of homogeneous reservoir properties and planar fracture geometry (Nordgren & others, 1972). Figure 1.9: Planar fractures are oversimplified and mesh dependent. a) from Wasantha & Konietzky (2017) b) from matthewpais.com 13 More complex fracture geometries can be handled in Discrete Fracture Network (DFN) models (McClure et al., 2016), Distinct Element Models (DEM) (Cundall & Hart, 1992; Zangeneh et al., 2015), and Interface Element Models (IEM) with fractures modeled as zero-thickness interface elements using penalty parameters (Ferronato et al., 2008) or Lagrange multiplier (Jha & Juanes, 2014; Jha, 2016). Finite-Discrete Ele- ment Method (FDEM) is another method recently proposed to model both flow and microseismic response of hydraulic fractures under non-linear elastic fracture mechan- ics (Zhao et al., 2014). The models, however, require specification of the crack path in advance both mathematically and numerically. Some of the Cohesive Zone Element methods (Gens et al., 1989), also used to model hydromechanical fracturing, may intro- duce mesh dependency problems (de Borst et al., 2004) because they treat fractures as planes of discontinuity that can propagate only along the boundaries of the mesh elements, see Figure 1.10. Figure 1.10: DFN models are also mesh dependent and require prior specification of the failure path. a) from Maxwell et al. (2016) b) from mve.com 14 Efforts have been made to overcome the problem of mesh dependency and model large-scale discontinuities (Khodabakhshnejad, 2016; Zhou & Molinari, 2004). For example, Extended Finite Element Methods (XFEM), which do not require the mesh to conform to fracture paths, have been used to model flow and deformation in porous mediawithasinglefracture(Mohammadnejad&Khoei,2013;Khoeietal.,2016;Fuma- galli & Scotti, 2013). Such methods are useful in modeling and predicting results of fracturing experiments in synthetic materials. However, their applicability to large- scale petroleum reservoir simulations (Shirangi, 2014) with strong heterogeneity over multiple length scales is yet to be proven. The Equivalent Continuum Model (ECM) approach (Duarte Azevedo et al., 1998) treats the fractured rock as an overall configuration with a constitutive law that models the effect of fracture sliding or opening on flow and deformation field variables, see Figure 1.11. Figure 1.11: ECM models have an advantageous ease of implementation and can be used to incorporate flow-related phenomena as well. a) from tu-freiberg.de b) from northwestern.edu 15 An earlier example of such constitutive behavior is linear elastic intact rock with a viscoplastic model for the shear stress-shear displacement-dilatancy relation and a viscoelastic model for the normal stress-normal displacement relation (Barton et al., 1985). Continuum Damage Mechanics (CDM) theory, which uses a similar approach, is also used to model fracturing-induced damage. Recently, CDM has been applied successfully to model hydraulic stimulation in enhanced geothermal systems (Pogacnik et al., 2016) and to design fracturing in unconventional oil and gas reservoirs (Shen et al., 2016). In this dissertation, we use an ECM approach to represent the evolution in the poromechanicalpropertiesofafracturedrock, asitallowsustointegratethemulti-scale physical processes of flow, deformation and failure at the computational scale within the framework ofexisting coupled geomechanical simulators used in the petroleum industry. 1.4 Manuscript Organization This dissertation begins with an overview of the importance and applications per- taining to the modeling of fluid flow in fracturing porous medium, challenges, and a summary of available approaches. We also highlight the need for having a robust com- putational platform for simulating this phenomena, which is made possible through use of effective continuum methods. Next, we focus on formulating the mathematical representation of our multiphysics hydraulic fracturing model. We start with building a coupled geomechanics and fluid flow model using a sequential scheme, which is employed as a platform for our contin- uum model. We then combine the principles of brittle fractures with those of ductile fractures to build a realistic elastoplastic material model. Combined mechanisms of flow enhancement in the elastic and inelastic regions are also taken into account using 16 a permeability model based on continuum damage mechanics. In order to demonstrate thephysicscomplianceofourapproach, wevalidateourmulti-mechanismmodelagainst available experimental data capturing the behavior of the porous medium both in terms of mechanical deformation & failure and fluid flow enhancement. Once our model is calibrated and prepared, we challenge the available formation fracture potential estimations which are solely based on brittleness indices. We illus- trate the underlying physics of fracturing processes, which is affected by numerous factors beyond rock brittleness, including but limited to: pre-existing crack orienta- tion, formation depth and stress magnitudes, as well as the anisotropy of the in-situ stress field. In addition to these elements, the yield stress and the degree of plasticity of rocks dominate the fracturing process after loss of material linearity in relatively ductile settings. We then highlight the practical implications of our comprehensive approach by calculating the breakdown pressure required at the hydraulic fracturing operational sites under a variety of conditions. We also exhibit how we can deploy our volumetric element model inside a field-scale computational simulation framework leading to better hydraulic fracturing results pre- diction and optimization. This introduces another important potential of our approach for field-scale implementation, where factors such as reservoir heterogeneity and well placement & interference can be highly influential at a variety of scales. Another overlooked yet crucial application of our ECM is related to the issue of induced seismicity. In Chapter 6, we address this topic by starting with an overview of fault rupture mechanics and the different mechanisms of induced seismicity. We then computationally implement our model in adjacency to a fault system which is located in proximity of an injection well to demonstrate how the coupling between geologic continuum damage and subsurface pressurization could alter the stress components acting on fault planes imposing a higher risk of earthquakes. 17 Chapter 2 Subsurface Multiphysics As described in Chapter 1, downhole pressurized injection processes often involve simultaneous and interacting processes of geomechanics, fluid flow through porous media, and fracturing as in Figure 2.1. Figure 2.1: The coupled multiphysics involved in hydraulic fracturing include fluid flow through media, geomechanical deformation, and rock fracturing. In this chapter, we elaborate on each of the individual physical phenomena involved in such complex processes and the coupling between them. 18 2.1 Governing Equations 2.1.1 Mechanical Equilibrium In this work, we follow a continuum representation of the subsurface systems (Mase, 1970), in which the solid matrix, the embedded fractures, and the fluid are modeled as an overall configuration rather than discrete particles. Under the assumption of quasi- static motion with small deformations, the linear momentum balance of the solid/fluid system can be written as: O.σ +ρg = 0, (2.1) whereσ denotes the Cauchy total stress tensor andg is the gravitational acceleration vector. ρ is the bulk density obtained from ρ =φρ f + (1−φ)ρ s , (2.2) in whichφ is rock porosity (the ratio of the void space volume to the reference bulk volume), andρ f andρ s are fluid and solid densities respectively. These parameters are all referred in the deformed (current) configuration. 2.1.2 Fluid Flow through Porous Media Poroelasticity The coupling between geomechanics and fluid flow through porous media states that the rises and falls in the pore fluid pressure induce changes in the effective stress and leads to deformation of the porous material. The deformation of the porous medium also affects the fluid content and pressure. Such coupled processes are exhibited in 19 phenomena such as ground subsidence due to groundwater withdrawal. We reduce the Biot’s theory (Biot, 1941) to the case of single phase flow of fluids under isothermal con- ditions in a poroelastic medium to obtain Equation 2.3. For more details on multiphase fluid transport in porous media, see Parker (1989); Sherafati & Jessen (2018). dm dt +O.w =ρ f f, (2.3) where m is the fluid mass content per unit reference bulk volume of the porous medium, t denotes time, w is the fluid mass flux defined as the flow rate of fluid in terms of mass per unit area of the porous medium in unit time: w =ρ f v, (2.4) givenv is the seepage velocity of the fluid with respect to the solid skeleton under deformation. Darcy’s law (Bear, 1988) describesv as the following: v =− k μ (Op−ρ f g), (2.5) wherek is the rock permeability tensor, μ is the fluid dynamic viscosity, and p is the pressure of the pore fluid. Following Jha & Juanes (2014), we present ζ, the increment in fluid mass content as: ζ := ∂m ρ f0 , (2.6) where ∂m is the variation of fluid mass content compared to a reference state, i.e. m−m 0 , and ρ f0 is the reference density of the fluid. Coussy (1995) uses the theory of 20 poroelasticity introduced by Biot (1941) to obtain the incremental form of the coupled fluid flow and geomechanics equations as: δσ =C :ε−bδp1 (2.7) ζ =bε v + 1 M δp, (2.8) where C is the 4 th order drained elasticity tensor of the porous medium, 1 is the 2 nd rank identity tensor given as: 1 = [1, 1, 1, 0, 0, 0] T (2.9) and ε is the linearized strain tensor defined as the symmetric gradient of the dis- placement vector,u: ε = 1 2 (Ou +O T u). (2.10) ε v is the volumetric strain obtained from the trace of the strain tensor: ε v =tr(ε) (2.11) It should be noted that in this chapter we follow Jha & Juanes (2014)’s conven- tion of positive tensile stress, allowing us to segregate the volumetric and deviatoric components of the strain tensor as the following: ε = 1 3 ε v 1 +e (2.12) 21 Our sign convention may change for the purpose of convenience in the next chapters. Equation 2.12 leads to the expression of the increment in volumetric stress, σ v , as: δσ v =K dr ε v −bδp (2.13) It is useful to define the effective stress tensor,σ 0 , as the responsible factor for the deformation of the skeleton of the porous medium: σ 0 =C :ε (2.14) Following Equation 2.7, the effective stress during single-phase fluid poroelastic deformation can be defined as: σ 0 =σ +bδp1 (2.15) The coefficients used by Biot to explain the coupling between fluid flow and mechan- ics, the Biot modulus, M, and the Biot coefficient, b, are expressed in terms of rock and fluid properties as: b = 1− K dr K s , (2.16) 1 M =φ 0 c f + b−φ 0 K s , (2.17) where K d r is drained bulk modulus of the porous medium, K s is the bulk modulus of solid grains, andc f = 1 K f is the compressibility of the fluid phase. As in Equation 2.7, we observe that the fluid mass content variations can be attributed to two factors: the pore space volumetric deformation (expansion or contraction) and the changes in fluid pore pressure. Under the small strain assumption (Coussy, 1995), and by applying 22 linearization from a reference state to the current state, we can rewrite Equation 2.7 as: σ−σ 0 =C :ε−b(p−p 0 )1 (2.18) and get: 1 ρ f (m−m 0 ) =bε v + 1 M (p−p 0 ) (2.19) We then obtain the mass balance equation for the fluid in fixed-strain numerical schemese (Kim et al., 2009) by substituting Equation 2.19 into Equation 2.3: 1 M ∂p ∂t +b ∂ε v ∂t +O.v =f. (2.20) Similar to Equation 2.18, we linearize the changes of volumetric total stress com- pared to the reference state as the following: σ v −σ v0 =K dr ε v −b(p−p0), (2.21) which leads to track porosity changes due to the variations in the volumetric stress component and pore pressure. Using Equation 2.19 and keeping in mind thatm =ρ f φ, we get: ρ f ρ f0 φ−φ 0 = b K dr (σ v −σ v0 ) + ( b 2 K dr + 1 M )(p−p 0 ) (2.22) Equation 2.21 allows us to rewrite Equation 2.20 this time in terms of the total volumetric stress and pressure for fixed-stress numerical solutions as: (Kim et al., 2009): ( b 2 K dr + 1 M ) ∂p ∂t + b K dr ∂σ v ∂t +O.v =f. (2.23) 23 It is observed that in Equations 2.23 and 2.20, the individual rock and fluid com- pressibilities appear in form of combined properties indicating the poromechanical cou- pling between geomechanics and fluid flow. Poroplasticity While the same rule of linear momentum under conditions of quasi-static conditions and infinitesimal deformation applies to the poroplastic region as in Section 2.1.2, the fluid flow equations pointed out in that section are valid for situations where the porous medium experiences elastic strain only. In certain other circumstances, on the contrary, rock may undergo inelastic deformation, which requires additional formulations. The elastoplasticity formulation, as given in Coussy (1995) and Armero (1999), is based on an additive decomposition of the total strain and the variation of fluid mass content into elastic (reversible) and plastic (irreversible) parts (Simo, 1992): ε =ε e +ε p , (2.24) and ζ =ζ e +ζ p , (2.25) where the superscripts e and p denote elastic and plastic respectively. From Equa- tion 2.15, we know that: σ =σ 0 −bp1 (2.26) Equation 2.14 indicates that the effective stress is related to the deformation field through the elastic strain tensor, which gives the following when combined with Equa- tion 2.24: 24 σ 0 =C : (ε−ε p ) (2.27) The plastic effective stress tensor,σ 0 p , is defined as: σ 0 p =σ +βp1, (2.28) whereβ is the coefficient which accounts for permanent changes in the porous space. Similar to Equation 2.12, we also split the stress tensor into a mean part,σ v =tr(σ)/3, and a deviatoric part,s as in the following: σ =σ v 1 +s. (2.29) The first and second invariants of the stress tensor are defined as: I 1 = 3σ v (2.30) and J 2 = 1 2 s :s. (2.31) Then we express the volumetric stress in terms of the volumetric strain as: σ v =K dr (ε v −ε p v )−bp. (2.32) From the discussions in the previous session, we relate the pore pressure the fluid mass content changes as: p =M(ζ e −bε e v ) (2.33) 25 In this dissertation, the matrix of the porous medium is assumed to be isotropic with embedded microcracks. This allows us to split the elasticity tensor,C as follows: C = (K dr − 2G 3 )1⊗ 1 + 2GI, (2.34) whereI is the 4 th rank symmetric identity tensor and G is the shear modulus. C can be further broken into: C =K dr 1⊗ 1 + 2GI 0 , (2.35) whereI 0 is the 4 th order deviatoric operator tensor given by: I 0 =I− 1 3 1⊗ 1 (2.36) Note that the deviatoric components of stress and strain tensors can be related to each other using the equation above as: s =I 0 :σ (2.37) and e =I 0 :ε. (2.38) In an elastoplastic formulation, three components are needed to completely satisfy the problem: a yield condition, a flow rule based on a plastic potential, and a hardening or softening hypothesis. The yield condition specifies the state of stress at which the plastic flow initiates and is generally given in the following form: 26 f(σ,p,h) = 0 (2.39) where h is an internal state parameter. The plastic potential function is given in the form: g(σ,p,h) = 0 (2.40) The evolution of the internal variables is given by the flow rule expressed in terms of the plasticn: ˙ ε p = ˙ λ ∂g ∂σ , (2.41) where ˙ λ is the plastic multiplier and ∂g ∂σ is the gradient tensor of the plastic potential function with respect to the total stress. Similarly, the flow rule for plastic porosity or plastic fluid mass content variation is ˙ ζ p = ˙ λ ∂g ∂p . (2.42) When the plastic potential function is identical to the yield criterion (f ≡ g), the plasticity is described as associated. Otherwise, it is non-associated. The final component needed is a hardening hypothesis for evolution of the internal state variable h, which describes how the yield condition and flow rule are modified during plastic flow. When the yield condition and flow rule remain constant during plastic flow (e.g., no hardening), the material is referred to as perfectly plastic. Plastic consistency conditions and the complementary Kuhn-Tucker load- ing/unloading conditions are required to be imposed on the yield condition and plastic multiplier as the following: 27 ˙ λ≥ 0,f≤ 0 (2.43) and ˙ λf = 0. (2.44) As far as it relates to a coupled flow of fluid (Armero, 1999), the conservation of mass of an incompressible fluid flowing through a porous medium in absence of any source or sink reads as ˙ ζ− k μ O 2 p = 0, (2.45) where k is the absolute permeability of the porous medium, a 1D form of the per- meability tensor given in Equation 2.5, assumed to be constant. To calculate the accumulation term, we the increment in fluid mass is written in terms of elastic and plastic volumetric strains starting with: dm =d(ρ f φJ), (2.46) where J is the determinant of the deformation gradient, obtained by: J = 1 +ε v (2.47) under infinitesimal deformation assumption. Assuming 1 +ε v ≈ 1, the formulations above lead to: dm =ρ f (φdε v +dφ) +φdρ f , (2.48) 28 The relationship between the volumetric plastic strain and the plastic porosity can be established by partitioning the strain as: (1−φ)dε v = (1−φ)dε sv +dφ, (2.49) where dε sv = dε e sv + dε p sv . Similarly, the volumetric Cauchy total stress is also partitioned into the volumetric matrix stress and the fluid pressure to give: dσ v = (1−φ)dσ sv −φdp (2.50) We know that dε e sv =dσ sv /K s , and we combined Equations 2.49 and 2.50 to get: (1−φ)dε v = 1 K s (dσ v +φdp) + (1−φ)dε p sv +dφ (2.51) (1−φ)( dσ 0 v K s +dε p v ) = 1 K s (dσ 0 v −bdp +φdp) + (1−φ)dε p sv +dφ (2.52) ( 1−φ K dr − 1 K s )dσ 0 v + b−φ K s dp + (1−φ)(dε p v −dε p sv ) =dφ (2.53) From b = 1− K dr Ks , we get 1 Ks = 1−b K dr leading to: dφ = b−φ K dr (dσ 0 v +(1−b)dp)+(1−φ)(dε p v −dε p sv ) = b−φ K dr (dσ v +dp)+(1−φ)(dε p v −dε p sv ) (2.54) Thus, Equation 2.48 can be written as: dm =ρ f (φ( dσ v +bdp K d r +dε p v ) +dφ) +φdρ f (2.55) The stress-strain relations for both elastic and plastic materials are elaborated fur- ther in Section 2.2. 29 2.1.3 Fracturing We follow the Continuum Damage Mechanics (CDM) principles elaborate in Kachanov (1986) and Lemaitre (1996), which allow for the evolution of rock properties to be represented through use of an Effective Continuum Model (ECM). We consider that the mechanical and flow processes associated with hydraulic fracturing are occur- ring within the elements of the domain. Therefore, we incorporate the micromechanics of fracture initiation and growth using a damage mechanics approach which allows us to account for changes in the mechanical behavior as well as permeability evolution inside each Representative Volume Element (RVE). As described by (Lemaitre, 1996), a mechanically damaged continuum can be rep- resented by conceptualizing a body containing discontinuities at the microscale. A homogeneous damage variable, D, is defined simply as the effective surface density of microdefects (Figure 2.2). Figure 2.2: Damage induced by axial load F on a cylindrical volume element in terms of the fractional areaδS D /δS occupied by micro-mechanical defectsδS denoted by the normal vector n after Lemaitre (1996). 30 At any pointM and direction plane identified by normal vectorn, ifδS is the area of the intersection of the direction plane of interest with RVE, and δS D is the effective area of the intersections of all microcavities lying on δS, then the value of the damage variable at this point in this direction is the local surface density of the intersections of these cracks and cavities as: d(M,n) = δS D δS (2.56) It follows from the above definition that 0<d< 1. (2.57) This means that the damage variable can change from d = 0 for an intact material to d = 1 for a failed material. However, it should be noted that in most practical applications, rocks intrinsically have an initial set of cracks and inclusions in them. Also, failure occurs prior to reaching the maximum value of d = 1 due to unstable crack growth. Therefore, it is useful to characterize the initial damage state of rock, d 0 , and define a critical damage variable, d c , above which rupture occurs. d 0 <d<d c (2.58) Infact, thelossintheeffectivesurfacethatresiststheloadisconceptuallyequivalent to a situation where the material is subject to an excess load (Rabotnov, 1969). In 1D, this excess load can be obtained by: ˆ σ = σ 1−d , (2.59) 31 where, σ is the nominal stress (the Cauchy stress tensor in 3D), and ˆ σ is the actual stress resisted by a surface inside the damaged material. Hooke’s law can relate the actual load to the observed deformation denoted by strain, ε, as follows: ˆ σ =Eε, (2.60) whereE is known as Young’s modulus, a material property. Equation 2.60 is equiva- lentto applyingthe nominalload,σ, toa damagedmaterial witha deterioratedYoung’s modulus, ˆ E, as: σ = ˆ Eε, (2.61) where ˆ E is given by: ˆ E =E(1−d). (2.62) Analogous to the stress tensor definition, the 3D equivalent of this phenomenon can be represented using a symmetric 2 nd order tensor defined in Equation 2.63, corre- sponding to the reduction of the surface area of the tetrahedron drawn in Figure 2.3. A damaged configuration can be whereby substituted by a strain-equivalent configuration. Note that tilde denotes the equivalent normal vector and surface area of the equivalent configuration. For use in non-proportional loading problems, the one-dimensional definition of the damage variable can be extended to multiple dimensions and anisotropic cases by defining the damage tensor: 32 Figure 2.3: Principal elements of the damage tensor in 3D (a) damaged configuration (b) strain-equivalent configuration denoted by tilde after Lemaitre (1996). d =        d 1 0 0 0 d 2 0 0 0 d 3        (2.63) Inviewoftheexperimentalevidence, fracturesinsideanoriginally-isotropicmaterial areprecededbythedevelopmentofanisotropicdamage, introducedinEquation2.63. A generalizedanisotropicmodelofcontinuumdamagemechanicshassincebeendeveloped by Lemaitre (1996) and applied successfully to characterize crack initiation in various materials. In this work, we follow a similar approach by assuming the host material to be isotropic and idealizing the shape of all microdefects inside each REV as lower dimensional entities. Each individual crack is characterized by its half-length, a, and a normal vector,n, which determines the orientation of the crack with respect to horizon, ψ as in Figure 2.4. In general, multiple crack families can reside inside each volumetric element. In the following sections, it is shown that it is useful to capture the additive damage effects caused by all crack families using Equation 2.64, where m denotes the total number 33 Figure 2.4: Microdefects in a 2D REV idealized as microcracks of a lower dimension (1D) (a) an individual crack characterized by its half-length,a, and angle with horizon, ψ (b) multiple crack families with various lengths and orientations distinguished by color. of crack families and m i is the number of family members for a family i. The total damage tensor is then calculated as: d = m X i=1 m i d i , (2.64) whered i is the damage imposed by single crack from family i obtained from: d i =d i (n i ⊗n i ), (2.65) in whichn i denotes the normal vector associated with this family of cracks, and d i is the scalar damage corresponding to rock surface deterioration created as a result of 34 crack growth. In 3D, the cracks are 2D planes, and such scalar loss of surface can be accounted for as: d i = ( a i a e ) 2 , (2.66) where a e is the half-length of the volumetric element used for calculations. Under the assumption of no crack re-orientation in the porous media, the only crack characteristic changing as a result of pore pressure increase, i.e. tensile loading, would be the length, 2a, whichalsoimposesitsresultinggrowthofdamageasinEquations2.64through2.66. In Section 2.2, we elaborate on how such damage growth leads to changes in rock’s mechanical and flow characteristics. 2.2 Material Models 2.2.1 Mechanical Deformation Elasticity Foralinearelasticmaterialsubjectedtouniaxialstressloadingconditions, thestress strain curveσ =f(ε) is linear, according to Hook’s law, as shown in Figure 2.5 and the 1D equation below it. σ =f(ε) =Eε (2.67) In the general form of this equation, for a 3D object, the stress and strain are six-dimensional vectors which are related via in Voigt notation as: σ =C :ε (2.68) 35 Figure 2.5: The stress-strain curve for a linear elastic material under uniaxial loading. The number of material constants in the stiffness tensor, C, is 21, corresponding to the most general anisotropic linear elastic material model. The expanded form of this equation is                    σ 1 σ 2 σ 3 σ 4 σ 5 σ 6                    =                    C 11 C 12 C 13 C 14 C 15 C 16 C 22 C 23 C 24 C 25 C 26 C 33 C 34 C 35 C 36 C 44 C 45 C 46 symm C 55 C 56 C 66                                       ε 1 ε 2 ε 3 ε 4 ε 5 ε 6                    , (2.69) where σ 1 , σ 2 and σ 3 represent the axial stresses, where as σ 4 ,σ 5 and σ 6 are the shear stresses. The strain variables, ε 1 ,...,ε 6 , are defined similarly. The shear strains here represent the engineering shear strains which are the sum of symmetric components, i.e., ε 4 = 2ε 23 =ε 23 +ε 32 , etc. ItisusefultointroducetheLameparameters,λandμ, intermsofYoung’sModulus, E, Poisson’s ratio as: 36 λ = Eν (1 +ν)(1− 2ν) (2.70) and μ = E 2(1 +ν) (2.71) The stiffness tensor for an isotropic host rock can be written now in terms of the Lame parameters as the following: C =                    λ + 2μ λ λ 0 0 0 λ + 2μ λ 0 0 0 λ + 2μ 0 0 0 μ 0 0 symm μ 0 μ                    (2.72) Using the definitions above, Equation 2.69 can be written in inverted form as:                    ε 1 ε 2 ε 3 ε 4 ε 5 ε 6                    = 1 E                    1 −ν −ν 0 0 0 1 −ν 0 0 0 1 0 0 0 2(1 +ν) 0 0 symm 2(1 +ν) 0 2(1 +ν)                                       σ 1 σ 2 σ 3 σ 4 σ 5 σ 6                    , (2.73) 37 which implies the following: C −1 = 1 E                    1 −ν −ν 0 0 0 1 −ν 0 0 0 1 0 0 0 2(1 +ν) 0 0 symm 2(1 +ν) 0 2(1 +ν)                    (2.74) In future sections, we will explain how the damage process affects intact material properties described above. Plasticity Themajorityofcoupledfluidflowandgeomechanicssimulationsperformedintheoil and gas problems only consider elastic material behavior. While such linea models are simple to implement, they have limited practical scope due to the vast observations of material yield in formations including sandstones, shales, and carbonates (Papanasta- siou, 1997; Rutter & Glover, 2012). Figure 2.6 shows some examples for these nonlinear stress-strain curves. In our work, we have found it imperative to take into account such important observations of material nonlinearity. Plasticity refers to the occurrence of irreversible strain in the material. It is impor- tant not to associate all type of material nonlinearity to plastic behavior, as in the case of nonlinear elasticity the stress-strain relation becomes nonlinear but otherwise the behaviour follows that of linear elasticity, meaning that loading and unloading are not distinct except for the sign. By contrast, in elasto-plastic materials, an n irreversible 38 Figure 2.6: Experimental stress-strain curves obtained for various rock types demon- strating inelastic behavior: a) Jurassic Shale after Papanastasiou (1997) b) Berea Sand- stone after Rutter & Glover (2012) c) Hollington Sandstone after Rutter & Glover (2012). plastic strain of ε p remains even after complete unloading due to having been sub- jected to forces of a significant magnitude. The difference between two such materials is illustrated in Figure 2.7. Figure 2.7: Stress-strain curves plotted for two different types of material nonlinearity: a) nonlinear elasticity b) and elasto-plasticity. 39 Plasticitytheorymakesuseofthreefundamentalconcepts: theyieldcriterion, which defines the limit at which the material becomes plastic, and the flow rule describing the relationship between the stresses and strains after the material has become plastic, and the consistency condition which prevents stresses from exceeding the yield limit. Yield conditions define not only the set of permissible stresses, but also the condi- tions for which plastic deformations can continue to occur. In general, this function has the following scalar form: f(σ,α) = 0 (2.75) whereσ denotesthesixindependentstresscomponentsandαdefinesnspecificmaterial parameters. TheyieldfunctionpresentedinEquation2.75definesasurfaceinthestress space, which is illustrated in Figure 2.8 for the two-dimensional case. Conventionally, f < 0 is defined as the inside of the yield surface and f > 0 denotes outside. The boundary, f = 0, defines the elastic limit, i.e. the set of maximal permissible stresses. Figure 2.8: Yield surface for a 2D stress space. 40 Any state of stress corresponds to a point either inside or on the yield surface. No point outside the yield surface is admissible, and is returned back to the yield surface using return mapping algorithms or other Newton-Raphson nonlinear methods (Chen & Han, 2007). In the case of loading perfectly plastic material, only a redistribution between the different stress components take place such that the stress point can be imagined as sliding along the yield surface. On the contrary, if hardening is present, the stress point will still remain on the boundary defined by f = 0, but the boundary itself is allowed to shift or change size and shape as the loading progresses in accordance with the relevant hardening rules. When unloading, the stress point moves from the boundary of the yield surface to the inside and recovers its elastic properties. See the stress-strain relationship of perfectly plastic and hardening materials during loading and unloading in Figure 2.9. Figure 2.9: Stress-strain diagrams under: a) Perfect plasticity b) Hardening plasticity. Hardening is often described by a combination of two specific types of hardening, namely kinematic hardening and isotropic hardening corresponding with shifting and growing/shrinking of the yield surface respectively, see Figure 2.10. 41 Figure 2.10: a) Isotropic hardening b) Kinematic hardening. In the following sections, we outline the most commonly used yield criteria defined for specific materials and applications. It is worth mentioning that one needs to ensure that a rotation of the coordinate system does not influence the conditions at which the material yields. Therefore, yield criteria are usually defined in terms of the stress invariants rather than the stress components themselves. See Equations 2.30 and 2.31 for their corresponding definitions. Mohr-Coulomb: The Mohr-Coulomb criterion is widely used in geotechnical engi- neering to define the shear strength of soils and rocks at different effective normal stresses. This criterion represents a linear envelope that is obtained from a plot of the shear strength of a material versus the applied effective normal stress: τ +μ f σ 0 n −τ c = 0, (2.76) where τ is the shear stress, σ 0 n is the effective normal stress, and μ f and τ c can be considered as a material property. This criterion is shown in the principal stress space 42 in Figure 2.11. For more details on the practical aspects of the Mohr-Coulomb criterion in earthquake analysis, see Chapter 6. Figure 2.11: Mohr-Coulomb criterion in principal stress space. von Mises: The von Mises criterion is commonly used to model the behavior of met- als, which behave independent of the hydrostatic pressure applied, see Equation 2.77. When plotted in a three-dimensional principal stress space, the yield surface associated with this criterion depicts a cylinder parallel to the hydrostatic axis σ 1 = σ 2 = σ 3 as shown in Figure 2.12. q J 2 −κ = 0, (2.77) where J 2 is the second stress invariant defined in Equation 2.31 and κ is a material property. 43 Figure 2.12: von Mises criterion in principal stress space. Drucker-Prager: In this work, we focus on the Drucker-Prager elastoplastic model (Armero, 1999). This criterion was originally devised to model plastic deformation of soils, where the hydrostatic pressure dependence cannot be ignored. Drucker-Prager yield criterion is a suitable choice for modeling granular materials under load, e.g. rock systems undergoing severe deformation. It provides with a smooth yield approximation of the Mohr-Coulomb criterion given in Equation 2.76. We elaborate on this model and apply it to the case of porous media having pore pressure as the key component of the added load. Figure 2.13 shows this criterion plotted in the principal stress space. The Drucker-Prager yield criterion for perfect plasticity is described by the following: f(σ,p) =α f (σ v +βp) + q J 2 −γ, (2.78) and the plastic potential is given as: g(σ,p) =α g (σ v +βp) + q J 2 , (2.79) 44 where α f , α g , β, and γ are model constants. Figure 2.13: Drucker-Prager criterion in principal stress space. Setting the effective stressI 0 1 = 3(σ v +βp) with the effective normal tractionσ 0 n and √ J 2 with the shear traction, τ, the Drucker-Prager model would become analogous to the Mohr-Coulomb criterion stated above by exploring the evolution of the material point in I 0 1 − √ J 2 and σ 0 n −τ spaces respectively. We calculate the increments in plastic strain and plastic porosity using the gradients of the f and g functions: df = ∂f ∂σ :dσ + ∂f ∂p :dp, (2.80) in which ∂f ∂σ = 1 2 √ J 2 s + α f 3 1 (2.81) and ∂f ∂p =α f β (2.82) 45 The gradients of the plastic potential, g, can also be calculated similarly. We now substitute Equation 2.79 in Equation 2.41 to get: ˙ ε p = ˙ λ( α g 3 1 + 1 2 √ J 2 s) (2.83) We can also obtain the rate of increment in the irreversible variation of fluid content as ˙ ζ p = ˙ λα g β. (2.84) The stress-strain relation in the plastic region, as opposed to that of the elastic one, takes a nonlinear form urging us to follow an iterative procedure at each time step. If dt =t n+1 −t n denotes the time step taken when solving the system of equations at the time t n+1 , then the increment for any given variable x is by definition given as: dx =x n+1 −x n = ˙ x n+1 dt, (2.85) with the derivative of time, dt calculated at the newest time step, t n+1 . Similarly, we can calculate the plastic strain increment as the following: dε p =dλ ∂g ∂σ | n+1 (2.86) From Equation 2.83, we know that for the volumetric part, dε p v =dλα g (2.87) and for the deviatoric part, de p =dλ 1 2 q J n+1 2 s n+1 , (2.88) 46 and from Equation 2.84, we get: dζ p =dλ ∂g ∂p | n+1 =dλα g β. (2.89) The Kuhn-Tucker conditions, when applied to our solution strategy, can be written as: dλ≥ 0,f n+1 ≤ 0,dλf n+1 = 0. (2.90) To obtain results at any time step t n+1 , we implement an elastic predictor-plastic corrector algorithm, which is a type of return-mapping algorithm. For conciseness purposes, we skip the numerical details of the Newton-Raphson scheme used in our workflow. Instead, we refer interested readers to Armero (1999). We start with a trial set of solutions (dσ trial ,dp trial ) using the current estimate of the total strain, ε n+1 and keeping the internal variables, i.e. the plastic strain and the irreversible fluid content, constant such that dε p = 0 and dζ p = 0. The yield condition, given in Equation 2.78 is then evaluated at this linearized trial solution, which is obtained under the assumption of poroelastic behavior. Then, we check to see if f trial,n+1 (σ n +dσ trial ,p n +dp trial ) < 0. If that is the case, then the material behavior during this time step has been elastic and the trial solution is accepted as the solution for thet n+1 step. Otherwise, the behavior has been plastic and a plastic strain increment has to be computed in order to return the stress state to the yield envelope. See Figure 2.14 for schematics of this iterative scheme in two dimensions. 47 Figure 2.14: Return mapping of the state of stress using an elastic predictor-plastic corrector algorithm in 2D. 2.2.2 Continuum Damage Evolution In this section, we review the principle of continuum damage upon which we build our crack growth model. The obtained results are used to account for the changes in material behavior in terms of flow and mechanics as discussed in the previous sections. Brittle crack growth Classicalfracture mechanics principlesassumethat thebulk behavior ofthe medium can be described with linear elasticity. Griffith (1921) developed the concept of the stress intensity factor, K I , to predict the stress, σ, near the tip of a crack with half- length a. K I =σ √ πa (2.91) A pre-existing crack grows until K I becomes equal to a critical value K IC , fracture toughness, which is assumed to be a material property. 48 K I =K IC (2.92) For brittle rocks containing cracks under increasing local tensile loads, Shao et al. (2004) proposed a model for crack propagation, in which increasing the fluid pressure acts as an opening force on the crack plane extending its length by opposing normal compressive forces. K I = (σ n +p +f(r)σ t ) √ πa (2.93) The amount of macroscopic deviatoric stress is translated to the local tensile stress using a scalar factor f(r) , where r is the relative crack length with respect to the critical crack length, 2a c . For practical purposes, this function can be approximated through calibration with experimental rock failure data or numerical micromechanical analyses. This scalar function, in fact, accounts for the proportionality factor between the applied field stress the local stress concentration around microcracks (Lu et al., 2013). An example of this function used in a simulation workflow will be shown later in Figure 5.2. Ductile crack growth Asclassicalfracturemechanicsprinciplesfailtogiveaproperviewoffracturegrowth inabroadrangeofrocktypes, mostcommonhydraulicfracturingsimulationapproaches onlyconsiderbrittlecrackgrowth. Thisoftenleadstooverestimationsoffracturelength and stimulated reservoir volume if the rock exhibits ductility. Crack growth on the plastic yield surface, can be quantified by use of a constitutive relation between the plastic strain and the crack length (Shojaei et al., 2014). We use a simplified form of the proposed equation as follows: 49 a i =a i0 +χ a i a c i (a ic −a i0 ) N X k=0 |Δε p | k ε f , (2.94) where N is the total number of loading steps, ε f is the ultimate fracture strain, χ is a material parameter, a i denotes the crack half-length in ith direction, and|Δε p | k is the equivalent plastic strain increment at loading step k defined as: |Δε p | k = s 2 3 Δε p ij ε p ij (2.95) 2.2.3 Mechanical Degradation In order to model the influence of fracturing on material properties, we pursue the approach presented by Kachanov (1958), which ignores any mode shape of microscopic damage growth and all micro-damage geometries, postulating damage as a macroscopic state variable. The damage variable affects the macro-mechanical properties of the rock in the form of the exertion of an elevation of the nominal stress, which provides a measure of the actual state of stress acting within each material element. Using the tensorial representation of the damage variable introduced in Equa- tion2.63, wehaveestablishedtherelationshipbetweentheactualstressandthenominal stress in Equation 2.59, which can be rewritten in tensorial form as: ˆ σ =M(d) :σ (2.96) wherethesymbol(:)indicatesthetensorproductcontractedontwoindices. TheM(d), the damage effect tensor, is expressed in the principal coordinate system of damage as: 50 M(d) =                     1 1−d 1 0 0 0 0 0 0 1 1−d 2 0 0 0 0 0 0 1 1−d 3 0 0 0 0 0 0 1 √ (1−d 2 )(1−d 3 ) 0 0 0 0 0 0 1 √ (1−d 3 )(1−d 1 ) 0 0 0 0 0 0 1 √ (1−d 1 )(1−d 2 )                     (2.97) where d 1 , d 2 , and d 3 are the principal damage components as defined in Equa- tion 2.63. Equation 2.97, when combined with the hypothesis of elastic energy (Chow & Wang, 1989), implies that the constitutive damage equation of elasticity coupled with damage can be derived as ε = ˆ C −1 :σ (2.98) where ˆ C −1 is denoted as the following in the principal coordinate system of damage: ˆ C −1 = 1 E                     1 (1−d 1 ) 2 −ν (1−d 1 )(1−d 2 ) −ν (1−d 1 )(1−d 3 ) 0 0 0 1 (1−d 2 ) 2 −ν (1−d 2 )(1−d 3 ) 0 0 0 1 (1−d 3 ) 2 0 0 0 2(1+ν) √ (1−d 2 )(1−d 3 ) 0 0 symm 2(1+ν) √ (1−d 3 )(1−d 1 ) 0 2(1+ν) √ (1−d 1 )(1−d 2 )                     (2.99) Compare Equation 2.99 with Equation 2.74 written for intact rock. 51 2.2.4 Flow Enhancement The growth of damage inside the porous medium not only creates new void space (i.e. porosity) by nucleating vugs and microcracks, but it also provides with a more convenient path for the flow of fluid, which is perceived as an enhancement in rock permeability. The relationship between the changes in porosity during deformation can be obtained by the following equation after Tate (2005): dφ =−dε v (1 +φ 0 ), (2.100) whereε v is the volumetric strain defined in Equation 2.11, andφ 0 is the initial porosity. Note that the negative sign in Equation 2.100 is in compliance with the sign convention treating compressive strains as positive. It has been observed experimentally, that the intensity and trend of fluid flow enhancement depend on the damage and deformation evolution process. In fact, per- meability enhancement slows down significantly after the onset of material plasticity (Peach & Spiers, 1996). Figure 2.15 shows the different paths of porosity and permeability evolution under deformation. The vast majority of available methods,however, relate permeability to porosity and do not account for deviations in permeability enhancement after yield as shown in Figure 2.16. 52 Figure 2.15: Permeability enhancement attenuation in the inelastic region, plotted after Peach & Spiers (1996). Figure 2.16: Permeability and porosity grow in different ways under strain. Porosity increases almost linearly with strain. Permeability enhancement, however, attenuates on the onset of material plasticity. Porosity-based and empirical permeability models cannot model this transition of material flow behavior. Figure is plotted after Popp et al. (2001). 53 Permeability anisotropy, i.e. directionality, is also neglected in all porosity-based and empirical permeability models. The directionality of permeability and damage in most extreme cases of 2D flow in the directions perpendicular and parallel to the fracture are shown in Figure 2.17. Figure 2.17: Directionality of permeability and damage: (a) Fracture and matrix are in series and the matrix permeability dominates flow (min permeability). b) Fracture and matrix are in parallel and the fracture permeability dominates flow (max permeability). In the following section, we introduce a continuum damage model for studying the enhancement of rock permeability with the growth of damage taking into account both the neglected anisotropic evolution of permeability as well as the dual patterns observed in the elastic and plastic regions. Porosity-based Models The famous Kozeny-Carman model (Kozeny, 1927; Carman, 1937) is among the most widely used methods to link permeability and porosity in a general loading space (Walsh & Brace, 1984). This model is expressed as: 54 K = φ 3 Bτ 2 S 2 , (2.101) whereK denotes permeability,φ is porosity,τ is tortuosity, i.e. the ratio of real flow path to the straight path connecting flow-in to flow-out, S is the surface area per unit volumeofporousmedia, andB isaporeshapecoefficient. Undercertaincircumstances, this model provides reasonable simulation results; however, it it includes a number of hypothetical parameters which cannot be easily measured or calibrated. Empirical Models Semi-empirical equations have also been proposed to improve the estimation of rock permeability subjected to various loading conditions. lnK =a K +b K φ +c K φ 0.5 , (2.102) where a K , b K , and c K are clay content-related coefficients, and φ denotes porosity. Empirical methods are mostly lithology-dependent and only work for limited types of rock. For a review of permeability enhancement models during rock deformation, read Ma (2015). Damage-based Models As the crack grows inside each element, a better path is provided to the fluid to flow, creating an enhancement of permeability. After Shao et al. (2005), the relation between the overall permeability tensor and the crack half-length can be rearranged to obtain the changes of the permeability tensor as in (Lu et al., 2013): Δk = 1 6 aR(a)e(σ 0 t ,n) 3 (δ−n⊗n) (2.103) 55 where R(a) is the connectivity coefficient, given by: R(a) =t 1 ( a−a 0 a c −a 0 ) t 2 (2.104) In Eq. 2.104, a 0 is the initial crack half-length, t 1 and t 2 are constants representing the connectivity characteristics, and e(σ 0 t ,a) is the average aperture of microcracks denoted by: e(σ 0 t ,n) =πaσ 0 t /E 0 (2.105) whereσ 0 t is the local tensile stress acting on microcracks, andE 0 is the initial elastic modulus of the rock. The local tensile stress imposed on the crack surface, is formulated by Costin (1983) as: σ 0 t = 2 π [σ kk + (S i − S 0 r 0 a)F (γ)], (2.106) whereS i is theith component of the deviatoric stress tensor defined in Equation 2.109. If more than one deviatoric stress component exerts an influence on a given crack, the resulting influence is adjusted throughF (γ), which represents orientation factors. This is equivalent to assuming that the crack is oriented with its normal along one of the three components of the principal stresses, giving a nominal value to F (γ) accounting for the statistical average orientation and tensile field strength over a finite region. S 0 is the nominal threshold stress and r 0 is the minimal tensile stress influence distance. S 0 andr 0 may both be considered material properties. Figure 2.18 shows the schematic diagram of the parameters in the model presented by Costin (1983). 56 Figure 2.18: Schematics of the local stresses and applied stresses after Costin (1983). Shao et al. (2004) reformulated the work by Costin (1983) into the following simpler form: σ 0 t =σ n +p +f(a)n.S.n, (2.107) where σ n =n.σ.n (2.108) is the total normal stress applied to the crack surface, p is the pressure exerted by the fluid, andS denotes the tensor of deviatoric stress, defined as: 57 S =σ−tr(σ)/3I, (2.109) It should be noted that the damage parameter may have different definitions intro- duced by other authors. It could also be observed that after a number of algebraic rearrangements of Equation 2.103 and 2.65, the permeability enhancement equation can be rewritten in the following form: Δk =g(σ,a)a 2 d(a), (2.110) where a is the crack half-length, d(a) denotes the length-dependent damage parame- ter, and g(σ,a) represents all other factors accounting for material properties, cracks connectivity, and stress concentration around the crack tip. 2.3 Summary In this chapter, we went through the governing multiphysics phenomena of the subsurface systems during the complex coupled processes of fluid flow, geomechani- cal deformation, and fracturing. We also presented the mathematical expressions of the evolution of rock properties when fractured under geomechanical loading and pore pressure changes using a damage mechanics road map. In the following chapter, we demonstrate our implementation of this effective continuum approach inside a compu- tational simulation workflow for modeling hydraulic fracturing processes. 58 Chapter 3 Modeling Workflow for Permeability Enhancement and Material Degradation during Pressurized Injection 3.1 Introduction The current energy independence enjoyed by the United States, which has been ben- eficial both economically and geopolitically, can be directly attributed to the immense developments in unconventional oil and gas resources. This in turn has been made possible by hydraulic fracturing and horizontal drilling technology. While it is indis- putable that great strides have been made, significant uncertainties and safety concerns remain in the hydraulic fracturing design process in the area of optimized stimulation performance and, therefore, investment success (Aminzadeh, 2019,?). Stimulation success is largely governed by the union of two important factors: in- depth knowledge about the reservoir and the accuracy of the models representing the physics which control fracture initiation and growth leading to improved production through increases in the "Stimulated Reservoir Volume" or SRV. Several multi-physics 59 models have been developed to predict reservoir response to high-pressure fluid injec- tion, in-situ stress interaction with changes in pore pressure and the performance of induced fractures during fluid production. An outstanding challenge, however, is yet to be addressed to fill the discrepancy between simulation model predictions and field observations. From a multiphysics point of view, hydraulic fracturing is a very complex process involving multiple phenomena at different scales. The effects of fluid injection on frac- ture creation is governed by the interaction of various parameters of the porous medium such as the initial state of stress, rock elastic moduli, rock compressibility, Poisson ratio, etc. (Jha & Juanes, 2014). Intuitively, a comprehensive hydraulic fracturing model should accurately output the fracture path and dimensions after the stimulation job. Instead, most available approaches, such as Discrete Fracture Network (DFN) models (Zangeneh et al., 2015), require prior knowledge of the fracture path as input to the model, thereby enforcing a set of pre-determined fracture growth characteristics into the model. In Chapter 2, we elaborated on the details about our Effective Continuum Model (ECM). In this chapter, we deploy the model we have developed inside a computational simulation workflow. We implement a subsurface pressurized injection scenario and predict the results of hydraulic fracturing. We also introduce methods on how this work can be used for better and safer design of such jobs. 3.2 Implementation of an Integrative Workflow Figure 3.1 demonstrates a concise form of the integrative continuum workflow we have developed to numerically simulate multi-mechanism fracturing processes during hydraulic fracturing operations. 60 Figure 3.1: An integrative workflow presented for the simulation of multi-mechanism hydraulic fracturing processes. 61 Our workflow begins with a problem definition step, in which the model domain geometry is delineated. In accordance with the computational simulation platform which accommodates a coupled finite-volume fluid flow and finite-element geomechan- ics simulation, the model is discretized spatially into a set of Representative Element Volumes (REV)s, and the mesh geometry is specified. We also have to initialize the model with REV properties in terms of the material model, permeability, and dam- age state, as well as poroelastic and poroplastic coefficients discussed in Sections 2.1.2 and 2.1.2. It is also important to account for prescribed loads applied in form of initial and boundary conditions. We summarize our proposed procedure in form of the fellow chart depicted in Figure 3.1 and as follows: 1. Downhole pressurized injection leads to a change in the effective stress, which needs to be represented as a series of discrete load increments. The gradual pre- scription of stresses or boundary displacement values ensures maintaining quasi- static conditions. 2. Depending upon the simulation platform used such as CMG or COMSOL Multi- physics, afullycoupledorasequentialapproachissoughttosolvethemultiphysics problem of fluid flow and geomechanics at this given time step. This coupled sys- tem consists of the mechanical equilibrium equation (as in Equation 2.1): O.σ +ρg = 0, (3.1) and the fluid flow equation, described for poroelastic conditions in Equation 2.23 as: ( b 2 K dr + 1 M ) ∂p ∂t + b K dr ∂σ v ∂t +O.v =f. (3.2) 62 To account for plasticity, Equation 2.55 may be used to describe the flow of fluid as following: dm =ρ f (φ( dσ v +bdp K d r +dε p v ) +dφ) +φdρ f (3.3) 3. For each REV inside the mesh, the state parameters, i.e. pore pressure, p, and the stress tensor components,σ, are calculated. (See Sections 2.1.2 and 2.1.2.) 4. As we are considering both the elastic and inelastic modes of deformation, it is important to check for the yield criterion to guarantee proper calculations of strain and damage. For most subsurface applications, it is often convenient to use the Drucker-Prager criterion described for perfect plasticity by Equation 2.78 as the following: f(σ,p) =α f (σ v +βp) + q J 2 −γ (3.4) (See Section 2.2.1.) 5. If the material at a given REV is determined to yield, then the return algorithm is used. 6. We then adjust for the proper stress and strain values using the material model and the corresponding return method. (See Section 2.2.1.) 7. If the state of the material is determined to not meeting the yield criterion, it is still important to check for the value of plastic strain,ε p . Not meeting the yield criterion can be due to lack of excessive loads beyond the yield stress from the beginning, giving the values of plastic strains equal to zero. It can also occur after multiple iterations of the return leading to the adjustment of stress, in which case the plastic strains would get nonzero values. 63 8. If there is no plastic deformation, then the brittle model, which holds the under- lying assumption of elastic host behavior, would be sufficient to model the growth of crack inside each REV. Brittle crack growth is described as in Equation 2.93: K I = (σ n +p +f(r)σ t ) √ πa (3.5) (See Section 2.2.2.) 9. If there is nonzero plastic strain, then the principles of linear elastic fracture mechanics cannot be applied, and the ductile crack growth model needs to be used as expressed in Equation 2.95: |Δε p | k = s 2 3 Δε p ij ε p ij (3.6) (See Section 2.2.2.) 10. The brittle and ductile constitutive models enable us to accurately update the lengths, a i , for each family of cracks inside the REV. 11. The damage tensor is then calculated to capture the overall effect of the existence and growth of all microcracks in each REV using Equation 2.66: d i = ( a i a e ) 2 , (3.7) (See Section 2.1.3.) 12. The damage approach us then used to update the hydromechanical properties of the REV at this time step. The elastic stiffness tensor and rock permeability are 64 adjusted after examining the degradation caused by damage components as in Equation2.99: ˆ C −1 = 1 E                     1 (1−d 1 ) 2 −ν (1−d 1 )(1−d 2 ) −ν (1−d 1 )(1−d 3 ) 0 0 0 1 (1−d 2 ) 2 −ν (1−d 2 )(1−d 3 ) 0 0 0 1 (1−d 3 ) 2 0 0 0 2(1+ν) √ (1−d 2 )(1−d 3 ) 0 0 symm 2(1+ν) √ (1−d 3 )(1−d 1 ) 0 2(1+ν) √ (1−d 1 )(1−d 2 )                     (3.8) and Equation 2.103: Δk = 1 6 aR(a)e(σ 0 t ,n) 3 (δ−n⊗n) (3.9) 13. The updated continuum at this time step is used as new material parameters for the next time step. The same steps as the ones mentioned in this list are taken to reach a new equilibrium for any future step. If there are no more load increments, the simulation concludes here. In the next section, we show how our continuum model is validated via comparison with lab observations. 3.3 Damage Evolution Model Validation We verify the implementation of our coupled continuum model by exploring the response of a single REV to deformation under simple uniaxial loading conditions. We examine the evolution of damage inside the REV, which leads to changes in rock 65 Property Nomenclature Value Confining Pressure P c 0 MPa Fluid Saturation S f 0 Ultimate tensile strength σ ut 3.4 MPa Critical damage for ψ = 0 o d c 0.25 Critical damage for ψ = 90 o d c 0.33 Table 3.1: Measured Parameters of Uniaxial Tensile Test. properties. Thevariationsoccurringinthemechanicalandfluidflowbehaviorisstudied under two separate sections. The workflow adopted for both calibration sections is the same as depicted in Figure 3.1. 3.4 Mechanical Degradation Jalbout (2006) has carried out experimental analyses for studying the damage evolv- ing in underground rocks under tension and compression. We use the data obtained from a number of their uniaxial tension measurements for comparison and calibration with our model. The configurations under study have a pre-existing set of deficiencies oriented either in the horizontal or vertical direction. See Figure 3.2 for the schematics of these exper- iments. The initially-damaged samples are subject to incremental tensile loading, and their mechanical behavior is recorded until the point of failure. Table 3.1 shows a set of observed parameters in selected tests performed by Jalbout (2006). The parame- ters which we have used in our model to calibrate with their data are summarized in Table 3.2. We consider linear elastic fracture mechanics principles for modeling brittle crack growth inside the isotropic host rock. When the tensile load increment is applied, the pre-existing cracks grow inside the REV. Figure 2.2 and Equation 2.56 show how the 66 Figure 3.2: Schematics of the pre-damaged samples with embedded fissures used in uniaxial tension tests. In the work by Jalbout (2006), ψ is set at 90 o for vertical cracks and 0 o for horizontal cracks. Table 3.2: Calibration Parameters for Uniaxial Tensile Test. Property Nomenclature Value Fracture toughness K IC 1.5 MPa √ m Poisson’s ratio ν 0.2 Initial crack length 2a 0 1 mm Element length 2a e 10 cm induced damage is defined as the directional surface density of cracks based upon the normal vector associated with crack planes. Figure 3.3 indicates the conceptual model in a schematic form. 3.4.1 Horizontal Cracks Figures 3.4 through 3.8 show the results obtained for an REV containing horizontal cracks undergoing a uniaxial tensile test. Figure 3.4 shows the relationship between the axial tensile stress, σ z and the axial strain, z . The growth and coalescence of the embedded fissures causes a non-linear 67 Figure 3.3: Conceptual schematics of crack growth inside an REV during a uniaxial tensile experiment. The incremental tensile load is applied in the z direction. The induced damage evolves depending upon the crack plane direction: a) Horizontal cracks induce damage in thez direction. b) Vertical cracks aligned perpendicular to they axis induce damage in the same direction. elastic behavior which deviates from the linear line drawn for intact undamaged mate- rial. In Figure 3.5, we are plotting the growth of pre-existing cracks inside the medium with incremental tensile loading. It should be noted that the damaged material is expected to fail before reaching the ultimate tensile strength, σ ut , measured for intact rock. Thisisconfirmedboththroughmeasurementsandthroughmodelingbyobserving that the cracks reach their critical length when the tensile load reaches the portion (1−d c ) of the nominal ultimate strength (See Equation 2.59.). Figure 3.6 shows the impact of damage accumulation on material mechanical prop- erties represented by stiffness tensor, C. See Equation 2.69 for the extended form of the C tensor and Equations 2.97 and 2.99 to see how we can obtain deteriorated stiffness elements from the damage tensor. 68 Figure 3.4: Stress-strain diagram for a damaging REV with horizontal fissures in response to increasing load compared to the linear behavior of intact medium. The material degradation can also be expressed in terms of the decrease of the Young’s Modulus, E, as in Figure 3.7. In Equation 2.62, it is shown that the rise of damage in a medium deteriorates Young’s modulus causing material degradation. We have also plotted the accumulation of damage perpendicular to the z direction against the growth of horizontal cracks from the initial to the critical state in Figure 3.8. 3.4.2 Vertical Cracks We follow the same procedure for vertical cracks this time. The results obtained for an REV with vertical cracks during a uniaxial tensile test are shown in Figures 3.9 through 3.13. In Figure 3.9, the relationship between the axial tensile stress, σ z and the lateral strain, x is plotted. Due to propagation of the pre-existing deficiencies in the REV, deviation from the original linear behavior for intact rock occurs. 69 Figure 3.5: Horizontal crack growth inside REV in response to increasing load. Critical crack length causes material failure prior to reaching intact rock strength. Figure 3.10, shows the growth of vertical fissures in the REV with increasing tensile load. Failure occurs later than the previous horizontal crack case and earlier than the intact medium original strength, σ ut . Crack propagation in the case of vertical loading and vertical cracks occurs mainly due to the effect of deviatoric stresses as in Equation 2.93. The effect of damage accumulation on material degradation is illustrated in Fig- ure 3.11. In this case, the first element of the stiffness tensor, C 11 , is shown to be affected. Figure 3.12 demonstrates the material degradation due to the growth of damage in terms of the changes in Young’s Modulus, E. Refer to Equation 2.62 to see the relationship between the damage parameter and Young’s modulus. Vertical crack elongation causes damage accumulation perpendicular to the x or y direction depending upon the assumed orientation. We plot the evolution of the 70 Figure3.6: Damageaccumulationinsidethesamplereducestheelementsofthestiffness tensor, changing the material’s mechanical behavior and causing deviation from linear elasticity. damage parameter w.r.t the growth of crack length from the initial to the critical state in Figure 3.13. 3.5 Permeability Enhancement Popp et al. (2001) has set up an experimental procedure for examining the evo- lution of porosity and permeability during triaxial deformation using samples from the Gorleben salt dome. Their results suggest that significant dilatancy is induced in these rock samples during such compaction experiments due to the differential stresses involved. It is also illustrated that unlike porosity, permeability evolves in different stages corresponding to stress and strain: a drastic permeability enhancement of up to 5 orders of magnitude occurs after the dilation of the pore space by only a small amount (< 1%). This sharp permeability enhancement in the elastic region is then fol- lowed by an almost constant permeability during the inelastic strain hardening region. 71 Figure 3.7: Young’s Modulus drops due to the loss of effective solid surface inside each element. See Equation 2.61. Data after Jalbout (2006). Their work suggests that the enhancement of permeability is both a function of poros- ity, which linearly evolves with dilatancy, and crack linkage, which slows down under inelastic deformation (Peach & Spiers, 1996). See Figure 2.16. We use the data provided by Popp et al. (2001) to calibrate our permeability enhancement model. Figures 3.14 and 3.15 show the results. It is clearly observed that the porosity evolution has a somewhat linear trend with strain, which is confirmed by Equation 2.100 emphasizing that any increase in the volume (pore space expansion) is translated into a rise in porosity. By contrast, permeability enhancement almost ceases after the onset of plastic deformation. Table 3.3 summarizes the parameters used for this modeling. 72 Figure 3.8: Damage accumulation caused by the growth of horizontal crack inside the REV starting from a negligible initial state to the critical value. Table 3.3: Parameters of Triaxial Compression Test. Property Nomenclature Value Fracture toughness K IC 1 MPa √ m Young’s modulus E 35 MPa Poisson’s ratio ν 0.2 Ultimate fracture strain u 0.15 Ductile crack growth parameter χ 0.6 Power law strengthening coefficient K 0.56 Power law exponent n 0.6 Initial crack length 2a 0 1 mm Element length 2a e 10 cm 3.6 Summary We discussed the mathematics of our coupled continuum model in Chapter 2. In this chapter, we demonstrated the compliance of our implemented model at the REV level with experimental observations. We showed how the growth of damage inside an REV leads to the extension of microcracks causing material mechanical weakness and 73 Figure 3.9: Stress-strain relationship corresponding to the growth of vertical cracks in an REV causing deviation from the intact linear response. permeability increase. We now move forward in the next chapter towards deploying our workflow inside a computational simulation platform for modeling pressurized injection for productivity enhancement as well as the induced seismicity implications. 74 Figure3.10: Thegrowthofverticalcrackshasaslowpacecomparedtothatofhorizontal fissures under incrementally increasing load, except for the late regions close to material failure. This is due to the fact that deviatoric stress components are the only driving mechanism for such deficiencies. See Equation 2.93. Figure 3.11: A deteriorating stiffness tensor component creating a non-linear response to axial load. 75 Figure 3.12: The loss of effective solid surface inside each element is translated into a decrease in Young’s modulus as in Equation 2.61. Data after Jalbout (2006). Figure 3.13: the damage component in the x direction evolves from a small value in the beginning to an ultimate critical value. 76 Figure 3.14: Porosity growth plotted against the axial strain ratio. Porosity grows linearly throughout the process, as it is only a function of the volumetric strain. See Equation 2.100. Data after Popp et al. (2001). 77 Figure 3.15: Multi-mechanism crack growth causes an increase in permeability. This enhancement significantly slows down when the material yields at the end of brittle crack growth region. Plastic deformation only allows small permeability growth. Data after Popp et al. (2001). 78 Chapter 4 Reservoir Simulation of Injection-Induced Permeability Enhancement 4.1 Introduction Hydraulic fracturing techniques have emerged as a promising technology for the development of unconventional oil and gas resources, where access to invaluable hydro- carbon deposits has only been granted through the creation of flow-stimulating frac- tures in the formation. The success of hydraulic fracturing jobs depends on appropriate design parameters based on the rock properties, e.g. brittleness, as well as a number of other operational factors. We elaborate on these matters in Chapter 5. In this chapter, we introduce a simulation model of pressurized injection scenario with an emphasis on ductile fracturing processes to demonstrate the importance of using appropriate models in obtaining realistic and proper estimates of the resulting reservoir enhancement. The results are published in Samnejad et al. (2017b). We develop a coupled scheme for solution of the multiphysics problem of fluid flow, rock mechanical deformation, and the damage process. Our scheme allows us to model various phenomena occurring during hydraulic stimulation on a single computational platform. It also allows us to investigate the effects of natural fractures, initial rock 79 structural anisotropy, and anisotropic in-situ tectonic stresses on these phenomena. Our approach has several advantages over the existing methods. It fits easily into the currently available platforms of reservoir flow and geomechanical simulators, and overcomes some of the inherent limitations associated with previous methods such as the need for apriori knowledge of the fracture propagation path. Furthermore, it has the ability to incorporate multiple physical processes while maintaining the robustness in simulation of flow and mechanical response of hydraulically fractured rock with a relatively low computational cost. 4.2 Mathematical Model In this section, we begin by describing our physics model of the damaged medium and the different stages of damage evolution reflecting on the material constitutive law. We define the damage variable and relate its numerical value to the fracture parameters of medium. We then study the associated variation of rock permeability and elastic modulus. Wealsopresentthecoupledequationsusedforsimulatingthedamageprocess of hydraulic fracturing. 4.2.1 Damage Model To be able to model the initiation of fractures at the continuum scale, we first need to grow the pre-existing microcracks within each RVE. Each element, as shown in Figure 4.1, can have one or more set of identical fractures, which are characterized by their length, L, (in 2D) and normal direction, n, and are randomly distributed within the boundaries of the RVE. The elements of the 2D damage tensor are associated with the overall effect of microdefects in corresponding directions. Following (Lemaitre, 1996), in which damage 80 Figure 4.1: Schematics of RVE with microcracks (a) element containing a single crack (b) element with two families of cracks. is defined as the ratio of the lost surface to the total surface, one can calculate the components of the 2D damage tensor as follows: D i = X k m( L j −L 0 L c −L 0 ) 2 k (4.1) whereD i is the damage parameter ini direction,k is the number of microcrack families, m is the population of each family, L j is the crack length in the j direction for each family,L 0 is a defined standard length for microcracks, andL c is a critical length which defines accelerated coalescence of microcracks leading to element failure. Note that directionality can be easily incorporated into the model by changing the direction of the normal vector n and adding a directionality parameter. Supplementing the above definition of the damage tensor with an empirical damage- permeability relation allows us to monitor both the rock mechanical behavior and the flow behavior throughout the process of hydraulic fracturing. The following sections describe how we use the damage variable as the key to coupling mechanics and flow during and after hydraulic fracturing stimulation. 81 4.2.2 Mechanical Model of the Damaged Medium The existence of microdefects inside each element influences material constitutive behavior. Damage mechanics allows us to relate fracture growth occurring as a result of stress concentration near the crack tip, depicted in Figure 4.2, at the micro scale to macroscopic processes. Figure 4.2: Stress concentration around a hole in a 2D stressed plate. In our model, the mechanical failure of rock consists of three stages: linear elastic stage with initial damage inclusion, material hardening in the plastic stage where the damage grows, and fracture formation and post-failure softening. These stages are shown schematically in Figure 4.3. 82 Figure 4.3: The process of damage growth inside an RVE due to increasing load from left to right:(a) Concentrated load is lower than yield stress; therefore, no crack growth happens. (b) Local plastic deformation leads to crack growth. (c) Concentrated stress exceeds material ultimate strength resulting in failure. In an n-dimensional model space, fractures are regarded as entitites of a lower dimension n-1. For simplicity purposes, we describe our model below in 1D. However, a 2D model with linear microcracks is implemented and employed in the simulations in this chapter. In the next chapter, we present and use a 3D model with planar fractures. Elastic region: Under relatively small loads, i.e. within the elastic region, the damageparametercanmodelmaterialbehaviorbyconceptualizingalossintheeffective surface that resists the load, which is equivalent to a situation in which the material is subject to an excess load (Rabotnov, 1969): ˆ σ = σ 1−d , (4.2) 83 where, σ is the nominal stress (the Cauchy stress tensor in 3D), and ˆ σ is the effective stress resisted. By direct state coupling through the concept of effective stress, one obtains the elasticity modulus of the damaged material as: ˆ E =E(1−D), (4.3) where, E is the Young’s modulus of the undamaged material, and ˆ E is the damaged Young’s modulus. No crack extension is considered within the elastic region in our model. Initial linear elastic parameters of each element are maintained throughout the elastic deformation region. Plastic region: By contrast to the elastic region, in the plastic region, as the stress exceeds the yield strength of damaged material, concentrated stress in the vicinity of the crack tips leads to crack extension and damage growth (Inglis, 1913). In our model, microcrack growth is associated only with plastic deformation. Shojaei et al. (2014) provides a way forward for modeling plastic strain effect on crack extension. Figure 4.3 shows this effect of stress concentration around crack tips schematically. Similar to material elastic constants, other material parameters, such as the yield strength, need to be modified in order to account for the existence of microcracks inside the material. ˆ σ y = σ y K T , (4.4) where ˆ σ y is the actual damaged yield stress, σ y is the yield strength of intact rock, and K T is the stress concentration factor. Failure and post-failure: For a crack to extend in an unstable manner, the stress intensity factor, K I , needs to reach a critical value known as the fracture toughness, 84 which is a material property. Thus, if this critical value, K IC , is reached, then a macroscopic fracture forms. K IC =Y s π L c 2 ˆ σ u (4.5) where Y is a geometry factor, L c is the critical crack length, and ˆ σ u is the ultimate strength of cracked rock. Following (Zhu & Tang, 2004), we model material softening after failure by incorporating idealized material residual strength. Again, it should be noted that these values need to be modified for a material containing cracks. Thus, material parameters are allowed to change as a function of fracture intensity and other fracture properties. At the macroscopic level, all the phenomena described above are captured by the material constitutive law shown in Figure 4.4. Figure 4.4: Element constitutive law capturing the three stages of mechanical damage and failure: Elastic region maintaining the original damage, plastic hardening allowing microcrack extension, and formation of microfracture with post-failure material soften- ing. In this figure, σ u is the ultimate strength of the material, σ y is the material yield stress, and σ r denotes the residual strength of the rock after failure. The material first goes through an elastic deformation stage during which damage does not evolve. When the yield stress is exceeded, plastic deformation is observed at 85 the macroscopic level while crack growth is happening at the microscopic level. When stress exceeds the ultimate strength, the medium fails due to unstable crack growth and formation of a macro-fracture. After the fracturing point, a softening behavior is implemented. 4.2.3 Fluid flow Model inside Damaged Medium Pressure effect can be modeled as an extra tensile load on the crack planes caus- ing material weakness. This is incorporated using the concept of effective stress first introduced by Terzaghi et al. (1996). σ 0 =σ +αpI, (4.6) where σ 0 is the effective stress tensor, α is the Biot coefficient, and I is the iden- tity matrix. Pogacnik et al. (2016) utilized a qualitative approach for matching the experimental data Wang & Park (2002) of permeability enhancement with a calibrated function. In our work, we use a similar sigmoidal function, shown below, honoring the mechanical model of rock as well as experimental observations of permeability enhance- ment. k i =k mat + 1 (1 + exp (−a(D j −b)) )(k fr −k mat ) (4.7) where k i is permeability in a direction, k mat is the matrix permeability, k fr is the fracture permeability, D j is the damage parameter in the other direction, and a and b are the calibration constants. Figure 4.5 shows the evolution of permeability with the directional damage parameter schematically. Permeability is assumed to be constant in the elastic region; which is, however, allowed to vary spatially as a function of the initial fracture density in the rock and the fracture orientation. Permeability enhancement 86 Figure 4.5: Permeability evolution with damage, where D 0 is the onset of permeability enhancement, and Dc is when a macro-fracture forms. begins as the stress exceeds the material yield strength, and reaches a certain maximum at failure, which is maintained afterwards in the post-failure region. 4.3 Numerical Experiment at the RVE Level We validate our model by conducting two simulations of hydraulic stimulation in a homogeneous isotropic reservoir. We assume that the reservoir is saturated with water at hydrostatic pressure and deforms under applied load following linear elasticity. We assume an initial distribution of cracks in the reservoir that is used to initialize the damage variable, which evolves during hydraulic stimulation. The elastic moduli, yield strength and fracture toughness evolve with damage following the mechanical model of damaged medium presented above. The permeability evolves following the flow model of damaged medium presented above. 87 Table 4.1: Uniaxial Tensile Test Parameters. Parameter Nomenclature Value Unit Youngs modulus E 1.0E+10 Pa Poissons ratio ν 0.2 - Matrix permeability k mat 5 md Fracture permeability k fr 100 md Fracture orientation θ π/6 - Initial crack length a 0 1 cm Critical crack length a c 4.5 cm 4.3.1 Example 1: Uniaxial Tensile Test Here we simulate a simple case of uniaxial extension in which the rock fails under tensile mode similar to what happens during hydraulic fracturing. Table 4.1 shows the parameters that we used for simulating a 2D tensile extension before failure, in order to validate our damage evolution model. Figs. 4.6 through 4.9 show simulation results. We performed the simulation for three different pore pressure values of 0, 10 MPa, and 30 MPa. Figure 4.6 shows the stress-strain relationship. The material goes through a linear elastic deformation, which is followed by plastic strain hardening. As the pore pressure increases, smaller tensile loads are sufficient for the element to yield and, ultimately, to fail. Figure 4.7 shows how a single crack inside RVE extends when tensile loads increase. In the elastic region, no crack growth happens, while plasticity extends the microcracks up to the failure point. Same extension behavior is expected for other cracks with other initial length and orientation. Increasing the pore pressure causes the crack to reach unstable critical growth earlier, i.e. at smaller loads. Figure 4.8 demonstrates the evolution of damage as the cracks inside RVE grow. For the case of multiple crack sets, combined effect of all cracks needs to be considered. 88 Figure 4.6: Stress-strain response. Additional pore pressure accelerates the damaging process. It is also interesting to see the directionality that crack orientation imposes on the evolution of damage tensor. Figure4.9showstheresultofdamageevolutiononpermeabilityenhancement, which models hydraulic stimulation. No permeability enhancement is observed within the elastic region. However, strain hardening, which happens afterwards, helps enhance permeability by causing the microcracks to grow. Ultimately, when a microfracture forms, permeability of the element reaches a maximum value. Again, by increasing pore pressure, the element approaches the failure envelope at the in-situ stress. Note that these results are plotted for a tilted microcrack; therefore, crack extension favors permeability enhancement in one direction, more than the other, as plotted in Fig- ure 2.17. 89 Figure 4.7: Crack growth. 4.3.2 Example 2: 2D Injection-Induced Stimulation Here, wesimulateinjection-inducedstimulationaroundawellborein2Dtoquantify the dynamics of permeability enhancement during and after stimulation. We consider a square domain of length 20 m with a circular well bore of radius 0.2 m. The initial permeability is 10 md, and microcracks are distributed inside the domain with an initial length of 0.01 m. Water is injected in the wellbore for 1 hour to increase pressure and stress in the surrounding rock and induce fracture growth and permeability enhance- ment. The domain is discretized with triangular elements as shown in Figure 4.10. The left and right boundaries of the domain are fixed in the normal direction, the bottom boundary is fixed in both directions, and a compression of 20 MPa is applied on the top boundary to create a compressive initial stress state. The maximum principal stress is in the y-direction. 90 Figure 4.8: Damage evolution. The simulation results are shown in Figure 4.11. Reservoir pressure and Von Mises stress increase as a result of injection-induced overpressure and shear stress around the wellbore. As the material yield strength is exceeded, plasticity occurs in proximity of the wellbore, whereas more remote locations of the reservoir continue to deform elastically. The plastic strain, plotted in Figure 4.12, is larger in the y-direction compared to x-direction because of the maximum principal stress direction. This results in extension of microcracks, shown in Figure 4.13, from an initial length of 0.01 m for the entire model to a maximum length of approximately 0.12 m. The damage parameter, also, evolves with microcrack growth as shown in Fig- ure 4.14. The permeability, depicted in Figure 4.15, increases as a result of damage evolution within each mesh element. The area of permeability enhancement is expected to grow in time as pressurization is continued. 91 Figure 4.9: Permeability enhancement. The spatial pattern of hydraulic stimulation around the well becomes anisotropic due to the stress anisotropy imposed on the medium through the boundary conditions. Figure 4.16 shows where permeability evolution with pressure at two mesh elements adjacent to the wellbore are plotted. Permeability is enhanced significantly in the areas exposed to higher Von Mises stress along the y-axis compared to the areas along the x-axis direction. In order to better understand the decline in permeability enhancement away from the wellbore and its sensitivity to in-situ stress anisotropy, we analyze changes of per- meability with distance from the wellbore along x and y directions in Figures 4.17 and 4.18 respectively. In both cases, permeability starts from a maximum at the wellbore, and decreases with distance to a minimum given by the initial permeability, 10 md. The enhancement is significantly more in the y direction due to the stress anisotropy. Finally, it is interesting to also compare the results obtained using the ductile model with those of the brittle assumption. As is seen in Figure 4.19, the models assuming brittlecrackgrowthcanoverestimatethepermeabilityenhancementresultsofhydraulic 92 Figure 4.10: Model mesh geometry and boundary conditions. fracturing jobs by orders of magnitude, resulting in tremendous wastes of investment capital. 4.4 Conclusions We presented a coupled simulation framework to model the processes of fluid flow, rock deformation, rock failure relevant for injection-induced hydraulic stimulation of petroleum reservoirs. Our framework uses continuum damage mechanics to capture the evolution of poromechanical properties such as permeability and stiffness and state variables such as stress and strain as the induced fractures grow with time. The frame- work does not require apriori determination of the fracture paths in the mesh because we capture the effect of microscopic crack growth on flow and deformation processes in terms of a macroscopic constitutive behavior that allows for elastic deformation, plastic 93 Figure 4.11: Reservoir state after an hour of injection. a) Pressure b) von Mises stress. deformation with hardening, tensile failure, and post-failure softening. We demonstrate the usefulness of our framework by simulating injection-induced stimulation around a wellbore in a 2D reservoir. 94 Figure 4.12: Effective plastic strain after an hour of injection around the wellbore. Figure 4.13: Microcrack growth after an hour of injection around the wellbore. Cracks are preferentially oriented along y axis because of the applied compression on the top boundary. 95 Figure 4.14: Damage evolution after an hour of injection around the wellbore. Figure 4.15: Permeability enhancement after an hour of injection around the wellbore. 96 Figure 4.16: Anisotropy in permeability enhancement. The two curves correspond to two elements in the mesh located at the co-ordinates shown in the legend. Figure 4.17: Permeability enhancement in the x direction, (a) cut line (b) permeability change with distance away from the wellbore. 97 Figure 4.18: Permeability enhancement in the y direction, (a) cut line (b) permeability change with distance away from the wellbore. Figure 4.19: Comparison of a) Crack Length, b) Permeability, and c) Damage resulted during pressurization of ductile (top row) and brittle (bottom row) material. 98 Chapter 5 Identifying Dominant Success Factors in Hydraulic Fracturing 5.1 Introduction The success of hydraulic fracturing jobs is often related to rock brittleness indices, which are taken as the sole impact factor determining fracturing results. Indeed, hydraulic fractures play a principal role in producing from low-permeability reservoirs; however, brittleness is not the only parameter contributing to productivity of unconven- tional resources. Under a variety of circumstances, brittleness indices are insufficient to explainrockfracture potentialandpermeabilityenhancementduringhydraulicstimula- tion. For better prediction and design, it is imperative to identify and understand other factors affecting fracture creation and propagation, and to build models that include the effect of these factors on flow enhancement. To numerically model permeability enhancement after injection, we can regard fractured rock as a damaged continuum, which allows simulation of the deformation and fracturing response of the reservoir using material constitutive laws for brittle and ductile regions. We outline a coupled flow-geomechanical simulation framework that fits into avail- able reservoir simulation platforms and does not require pre-specified fracture paths. We develop the fracture growth mechanisms for the coupled simulation framework by 99 analyzing the effect of rock properties and in-situ stresses on the fracture length at dif- ferent injection pressures. Based on these mechanisms, we propose factors that quantify the success of hydraulic fracturing jobs beyond the simplified rock brittleness indices. In this chapter, we examine reservoir post-stimulation flow behavior by numerically simulating rock permeability enhancement from hydraulic fracturing while considering variable rock mechanical properties and reservoir in-situ conditions. We model the rock and evolution in its poromechanical properties using the principles of Continuum Damage Mechanics (CDM), which has proven promising for geothermal applications (Pogacnik et al., 2016). The extension of our methodology (Samnejad et al., 2017b) exhibits the benefits of applying CDM approaches to model hydraulic stimulation of unconventional reservoirs and quantify various factors affecting the success of hydraulic fracturing jobs. Our technique has the advantage of fitting into commercial coupled simulators of flow and geomechanics, overcome the mesh dependency problem existing in previous methods, and honor principles of fracture mechanics for brittle rocks and plastic damage growth for rock ductility. The results are published in Samnejad et al. (2018). 5.2 Conventional Brittleness Indices 5.3 Mathematical Formulation Our goal here is to demonstrate how we can build a comprehensive hydraulic frac- turing model that can later be fitted into coupled flow and geomechanical platforms. Therefore, it is useful to discuss the numerical simulation applications of our model. If we regard this problem from a reservoir engineering standpoint, the domain is divided into volumetric elements and each element is undergoing changes in the stress state due 100 to pressurized fluid flow across the element. Within each element, continuum damage mechanics principles can be used to describe the variations of reservoir parameters. 5.3.1 Governing Equations of Reservoir Multi-physics Pro- cesses In order to simulate hydraulic fracturing processes, multiple physical phenomena need to be modeled. Well stimulation jobs include fracture initiation and propagation in a subsurface environment under tectonic stresses, leading to enhancement of rock permeability. Therefore, we present this complex problem as a coupled flow and geome- chanics simulation framework, incorporating fracturing processes within each volume element of the grid. 5.3.2 Fluid Conservation In presence of geomechanical deformation, Darcy flow of a single-phase compressible fluid in the reservoir can be described as: ( α−φ K s + φ K f ) ∂p ∂t +α ∂ v ∂t +O.[ k μ (−Op +ρ f g)]−Q p = 0 (5.1) where α is the Biot coefficient, φ is porosity, K s and K f are bulk moduli for the solid and fluid phases respectively, p is pressure,t is time, v is the volumetric strain,k is the rock permeability tensor, μ is the dynamic fluid viscosity, ρ f is the fluid density, g is the gravitational acceleration vector, and Q p is the source term due to injection in the volume element. It should be noted that when solving this equation to model hydraulic fracturing processes, the rock properties cannot be assumed constants and should be updated 101 as functions of stress at each time step. This is in contrast to traditional reservoir simulations with constant permeability, porosity, and Youngs modulus. 5.3.3 Mechanical Equilibrium For quasi-static deformation, the mechanical equilibrium equations can be written as the following: O.σ 0 +ρg = 0 (5.2) in which σ is the total stress tensor, and ρ is the bulk density. The total stress tensor can be related to the effective stress tensor, σ 0 , with the following equation: σ 0 =σ +αpI (5.3) whereI is the identity matrix. The effective stress tensor is related to the strain tensor using the material constitutive laws. We assume that normal stresses are positive under tension. 5.3.4 Continuum Damage Mechanics Model for Element Behavior We use a continuum damage mechanics workflow for modeling the change of hydromechanical properties of the rock, namely, permeability and stiffness, which influ- ence the coupled behavior of flow and geomechanics. In our workflow, a pre-existing set of mico-inclusions within each 2D representative volume element can extend under bi- axial loads caused by fluid injection. Linear extensions of such micro-inclusions lead to 102 increase in flow permeability and change in material stiffness. Figure 5.1 schematically shows a single crack inside a 2D element. Figure 5.1: Schematic view of a volumetric element. To capture the overall effect of existence of N different families of cracks in the element, withm k members in thek t h family, we define a damage parameter that sums up the relative lost surfaces associated with individual cracks of each family, D = N X k=1 m k a a c 2 k (n⊗n) (5.4) where n denotes the normal vector to the crack plane, and a c is the critical crack half-length above which unstable crack growth occurs. It is important to note that damage parameter evolution, which depends on crack extension, controls permeability enhancement. Therefore, if we capture crack extension properly, we would be able to forecast reservoir post-stimulation performance. We use principles of brittle fracture mechanics and ductile damage growth to trace the evolution of crack lengths. For 103 simplicity purposes, we assume one crack family per element with one member in the family, i.e. N = 1, m k = 1. In contrast to common approaches, which overlook non- brittle inelastic crack growth, we include multiple mechanisms of brittle and ductile crack propagation. Refer to Section 2.2.2 for a detailed discussion about how we model the evolution of the crack length in accordance with rock deformation and failure. In our multi- mechanism damage growth model, cracks grow in a brittle manner when material undergoes elastic deformation as in Equation 2.93. Figure 5.2 shows the evolution of the f function we have chosen for this equation, aiming to attenuate sliding effects as the critical crack length is approached. The ductile fracture growth matching with material plasticity is obtained from Equation 2.94. Figure 5.2: Variation of scalar function f(r) with relative crack length r. We use a simplified power law constitutive model of plasticity to relate strain to stress, as shown in Equation 10 for uniaxial deformation. 104 σ−σ y0 =kσ y0 n p (5.5) where σ y0 is the yield stress, k is the power law multiplier, and n is the power law exponent. We found it more convenient in our applications to work with threshold pressure, P t , as the onset of plastic yielding. See Equation 5.3 to see how it relates to stress. We demonstrate the effect of certain material properties on the crack growth beyond the elastic limit. We neglect the changes of the final fracture strain, f , with in-situ stresses and rearrange Equation 2.94 to solve for a i explicitly: a i = a 0 1− χ a ic (a ic −a i0 ) P N k=0 |Δp| k f (5.6) One may note that the initial crack length after the onset of plastic yielding is, in fact, the maximum crack length achieved during linear elastic deformation. Incorporating rock ductile behavior introduces nonlinearity to the problem. It enables us to go beyond existing linear elastic fracture mechanics approaches and quan- tify fracturing processes in formations that exhibit brittle fracturing only partially. The next sections provides an insight on anumber of important impact factors that are often overlooked in common analyses of fracturing processes. 5.4 Multi-Mechanism Crack Growth Modeling We use the material presented in the discussion section above to illustrate the evo- lution of brittle crack lengths within each volume element with increasing borehole pressure to better understand the parameters that influence the creation of fractures in the reservoir, and ultimately, permeability enhancement as a result of fluid injection. 105 Table 5.1: Brittle crack growth parameters. Property Nomenclature Value Biot coefficient α 0.9 Crack orientation ψ 60 o Fracture toughness K IC 0.7MPa √ m Initial crack half-length a 0 1mm Critical crack half-length a c 50cm Young’s modulus E 100GPa Poisson’s ratio ν 0.2 Horizontal Stress σ 1 25MPa Vertical Stress σ 2 35MPa Table 5.2: Ductile crack growth parameters. Property Nomenclature Value Biot coefficient α 0.9 Crack orientation ψ 45 o Power law exponent n 0.7 Power law multiplier k 2 Proportionality factor χ 1 Fracture Strain f 0.05 Initial crack half-length a e max @ brittle Critical crack half-length a c 50cm Threshold pressure p t 30MPa Horizontal Stress σ 1 25MPa Vertical Stress σ 2 40MPa Table 5.1 summarizes the reference parameters that we have used for brittle fracture growth, based on data in part provided by Gowd & Rummel (1980) for sandstone. Table 5.2 is a similar list of parameters when the increasing pressure in the wellbore corresponds to a rising stress field above yield strength of the rock. We identified five important factors that can control the growth of cracks in the reservoir in addition to simplified brittleness indices: alignment of inclusions, confine- ment load, load biaxiality (triaxiality in 3D), degree of ductility, and the onset pressure of yield. Sections 3.1 through 3.5 elaborate on these by showing plots of the extension 106 of crack lengths, a, with respect to the injection pressure considering different impact factors. 5.4.1 Pre-existing Cracks Orientation Rocks undergo numerous geological processes, which lead to preferred orientation of theirinclusionsandflaws. Theseflawsbearportionsoftheremotestressfielddepending on their orientation as shown in equations 2.107 through 2.109. Therefore, understand- ing the initial crack orientation helps better quantify hydraulic fracturing job results. The effect of initial crack orientation is shown in Figure 5.3, where the propagated crack half-length is plotted against the injection pressure. Considering the maximum principal stress being in the vertical direction (see Table 5.1), the initial cracks are expected to prefer the vertical alignment to the horizontal alignment (ψ≥ 45 o ). Figure 5.3: Pre-existing crack orientation effect on brittle crack growth. As expected, the higher the number of cracks aligned normal to the minimum prin- cipal stress direction, the faster is crack propagation, i.e., the cracks reach a certain 107 value of half-length at lower injection pressures. In this work, we have not included the effect of fracture re-orientation with respect to the minimum principal stress, but it can be studied using models as presented in Ashby & Sammis (1990). 5.4.2 In-situ Stress Magnitude Besides transition to ductile behavior, in-situ stresses play an important role in the design and outcome of hydraulic fracturing (Samnejad et al., 2017b). Higher confining stresses impose larger compression on the crack faces against each other thereby pre- venting further growth. Therefore, hydraulic fracturing in deep formations can be very challenging due to very high downhole pressures required for growth. Figure 5.4 shows the impact of maximum principal stress on brittle crack growth. To avoid redundancy, another figure for minimum principal stress is not included since the impact is similar. In general, lower magnitudes of in-situ stresses allow easier fracturing in brittle rocks. As can be seen on the plot, 15 MPa of stress drop results in significant crack growth stimulation. This observation can have important implications when identifying suitable fracturing intervals. Rocks behaving the same in the lab may not exhibit similar fracturing behavior in the field, if placed at different depths, i.e. under different loads 108 Figure 5.4: Maximum principal stress effect on brittle crack growth. 5.4.3 In-situ Stress Anisotropy We analyze the stress anisotropy effect in Figure 5.5 by looking into the stress anisotropy ratio, λ, defined as the ratio of the two principal stress magnitude values: λ = σ min σ min (5.7) Stress biaxiality (increasingλ) has negative effects on crack growth, and presence of larger stress anisotropy could contribute to better fracture development. Note that if the minimum principal stress becomes comparable to the maximum principal stress, the compressional forces of crack closure become dominant, and higher pumping pressures will be required to initiate noticeable crack growth. 109 Figure 5.5: In-situ stress anisotropy effect on brittle crack growth. 5.4.4 Degree of Ductility It has been shown experimentally that permeability can improve up to an order of magnitude during ductile deformation. This motivates us to understand the crack growth mechanism associated with plasticity. Here, we demonstrate that although a crackpropagatesmoreslowlyontheyieldsurfacecomparedtothelinearelasticdomain, the process of irreversible fracturing continues. In Figure 5.6, the closer the n exponent is to unity, i.e. linearity, the higher is the crack growth inside the medium. Note that a bottomhole pressure of 30 MPa is assumed to be the onset of plastic yielding. At lower pressures, linear elastic material models are sufficient to explain deforma- tion and fracturing phenomena. After the onset of plastic yielding, depending upon the degree of rock ductility, fracture growth process slows down to a certain extent. 110 Figure 5.6: Ductility effect on plastic crack growth. 5.4.5 Yield Stress If the rock yields after large effective loads are imposed, then plastic strain only happens at very high pressures, and a substantial portion of the crack growth occurs in a brittle manner leading to higher ultimate crack lengths. This can be shown by plotting the effect of threshold pressure, pt, below which rock behaves elastically, on crack growth as shown in Figure 5.7. 5.5 Operational Implications of Impact Factors In this section, we demonstrate how the factors analyzed in the previous section can dictate operational conditions that are required to obtain desirable fracture growth. More specifically, we are analyzing the relative extension of created fractures with respect to their original length, a i , under a given set of reservoir and operational con- ditions. We study how the required breakdown pressure differs depending upon the aforementioned factors. 111 Figure 5.7: Effect of the threshold pressure on plastic crack growth. 5.5.1 Pre-existing Cracks Orientation Fractures oriented in favorable (optimum) directions, i.e. perpendicular to the min- imum principal stress, can grow more rapidly compared to other fractures. This can be seen in Figure 5.8, where even at relatively lower pressures, fractures of extended length are formed, primarily because of their preferred direction. Note that we assume an initial crack length of 1 mm and a critical length of 50 cm (Table 5.1). 5.5.2 In-situ stress magnitude Hydraulic fracturing operations in highly confined media can be very challenging because the required pumping pressure to create fractures with desired length increases with the in-situ stress magnitude. Figure 5.9 shows this phenomenon, where only increasing values of pressure can cancel out the effect of high overburden loads. A 112 Figure 5.8: Effect of the direction of pre-existing inclusions on the required injection pressure. similar trend is expected to be observed for the minimum principal stress, which is not included here to avoid redundancy. 5.5.3 In-situ Stress Anisotropy We emphasize on the observation made in section 3.3 by showing how increasing stress biaxiality imposes higher required injection pressures to achieve desired fracture length. Figure 5.10 shows the dependency. 5.5.4 Degree of Ductility To better demonstrate the effect of rock ductility on required injection pressures, we have chosen an early transition of the rock fracturing behavior from brittle to ductile, at an injection pressure of 30 MPa, when not much brittle crack extension has taken place. Figure 5.11 shows that for a predominantly ductile rock, even high pressurization 113 Figure 5.9: Effect of the magnitude of in-situ stress on the required injection pressure. cannot cause significant crack coalescence. This phenomenon is less pronounced when the rock exhibits less ductility (higher n value). 5.5.5 Yield Stress Earlytransitiontoductilebehaviorpreventsrapidcrackextension. Weshowthatfor small values of the threshold pressure, i.e. the pressure at which the rock yields, crack extension occurs very slowly due to ductile transition. On the contrary, if the major portion of crack growth happens in the brittle region, i.e. a higher threshold pressure, rock breaks at even low injection pressure. Figure 5.12 depicts this phenomenon. 5.6 Conclusions We identified and quantified the impact of a number of rock mechanical and in-situ parameters on fracture propagation during hydraulic fracturing. These effects have 114 Figure 5.10: Effect of stress anisotropy on the required injection pressure. Lower values on the x-axis indicates higher anisotropy. not been addressed in existing models of fracturing operation success prediction. We developed a multi-mechanism model of crack growth that incorporates the effects of these factors on the fracture length and the required injection pressure, both of which determine the extent of success in a given hydraulic fracturing operation. Our model is advantageous over existing DFN model-based approaches because it does not require prior knowledge of the location and path of created fractures since it predicts the growth of fractures. Our model is capable of incorporating multiple mechanisms of fracture propagation arising from multiple independent rock properties and in-situ stress conditions. Since we use a continuum damage mechanics approach, ourmodeldoesnotencountermeshdependencyproblemsduringnumericalsimulations. We conclude that the initial alignment of pre-existing inclusions, absolute and rel- ative magnitudes of in-situ stresses, degree of ductility, and plastic yield stress are 115 Figure 5.11: Effect of the rock ductility parameter on the required injection pressure. important factors that contribute to growth of fractures and, therefore, to success of hydraulic fracturing jobs. Future work is aimed at implementing our model into a numerical simulator to demonstrate the applicability of our method to reservoir-scale simulations and highlight the robustness of workflow. 116 Figure 5.12: Effect of the threshold pressure on the required injection pressure. 117 Chapter 6 Fault Damage Zones Effects on In-situ Geomechanical Stresses 6.1 Introduction Attributed to anthropogenic causes, the recent rise in earthquake activity has become an important topic highlighting the need for prescriptive frameworks of hazard mitigationwhicharebaseduponaccuratepredictivemodels. Fluidinjectionoperations, which are pivotal in a variety of subsurface exploitation applications including nuclear waste disposal, oil and gas production, and geothermal energy extraction, are regarded as activities which can potentially trigger seismicity on critically-stressed faults. It has been proven that the injection-induced pore pressure diffusion anomalies in fault damage zones can significantly alter the stress perturbation patterns acting on fault systems. Therefore, it is crucial to understand the coupled nature of the fluid flow and geomechanical processes both in the intact host rock and in the damaged zones located adjacent to fault planes. Despite the importance of the role of such geologic complexities in determining the fault in-situ stresses and the associated seismicity risk, previous modeling efforts addressing fault system multiphysics lack any representation of damage zones inside the subterranean environment. In this chapter, we compute the spatiotemporal evolution of the pore pressure and the stress state acting on fault systems using a poroelastic simulation model which 118 takes into account the anomaly created by fault damage zones. We illustrate the stress- dominating role of such zones by assigning a set of distinct hydromechanical properties to fault damaged zones and the intact formation. We deploy a continuum workflow in our computational scheme, allowing to integrate multi-scale physical processes of flow and deformation. Our results indicate that the presence of damaged zones in proximity to fault planes can perturb the effective stress state acting on the fault system significantly. Fault damage zones result in a more pronounced reduction of the effective normal stress, leading to an accelerated shift of the in-situ state of stress towards the rupture envelope. Theextentofsucheffectsdependsuponthegeometryoffaultdamagezones, theamount of hydromechanical irregularity in them, and the transmissibility of the fault itself. We illustrate that overlooking the existence of such complexities in the subsurface rock systemsinourmodelsresultsinunderestimatingthechanceoftheoccurrenceofinduced seismicity, which may lead to dramatic casualties suffered by the public. Finally, we expect to shed light on the key for safe fluid injection operations in terms of well placement, flow rate, and pressure under various in-situ conditions by demonstrating the influence of including fault damage zones in our predictive models. In this chapter, we aim to shed light on the physics of induced seismicity by linking geomechanicalmodelswiththoseoffluidflowthroughporousmediahonoringthemech- anisms that may cause stress state changes leading to potentially seizable earthquakes. The result of this work is anticipated to help mitigate seismicity risk while also taking advantage of the economic benefits of various subterranean exploitation activities. The results are presented in Samnejad et al. (2018). 119 6.2 Overview of Seismicity 6.2.1 Mechanics of Fault Rupture Toassesstheriskofinducedearthquakes,itisnecessarytounderstandthemechanics ofrupture, happeningonexistingfaultplanes. Earthquakesinvolvethereleaseofelastic strain energy, which is typically stored when the tectonic stress field imposes traction on fault surfaces, as depicted in Figure 6.1. Several methods have been proposed to estimate and characterize the in-situ tectonic stresses in terms of magnitude and direction (Samnejad et al., 2017a). Figure6.1: Schematicof2-Dfaultplanes(hatchedsurfaced)inside3-Dhostrockblocks. It is often convenient to break the traction vector, t, into its normal, σ n , and shear, τ, components for a fault with a normal direction indicated by n. The Mohr-Coulomb failure criterion (Coulomb, 1776) is among the most common methods for describing the failure condition to initiate slip. It states that a fault system remains in equilibrium as long as the applied shear stress, τ, is lower than a critical limit value, τ c : τ <τ c (6.1) 120 According to the Mohr-Coulomb failure criterion, the value of the critical shear stress, which represents the strength of the contact, is given by: τ c =τ 0 +μσ n (6.2) whereτ 0 is the cohesive strength of the sliding surfaces being negligible for faults under typicalcrustalconditions(Zoback,2007)andμisthecoefficientoffrictionlyingbetween 0.6 and 1.0 for almost all rock types. σ n denotes the effective normal stress applied on fault surfaces. Considering that in addition to waste disposal sites, a considerable number of energy development sites also perform operations which involve injecting a type of fluid into geologic formations and/or producing from underground reservoirs, pore pressure is an important element of any subsurface analysis pertaining to induced earthquakes. In presence of underground pore pressure, the effective normal stress component can be obtained using the following equation: σ n =S n −P (6.3) where S n is the total tectonic normal stress and P is the pore pressure. The friction coefficient,μ, which is the slope of the failure enveloped in Figure 6.3, can be expressed in terms of the angle of internal friction, φ: μ =tan(φ) (6.4) 121 6.2.2 Mechanisms of Induced Seismicity Generally, two main mechanisms can be postulated to induce seismicity as a result of fluid injection: 1. The rise in the pore pressure caused by fluid injection causes a reduction in the effective normal stress. The critical stress state created can potentially trigger movement along a near-by fault leading to induced earthquakes (J van der Elst et al., 2013). 2. Change of the gravitational loading conditions caused by fluid injection could also perturb the stress state in the surrounding rock and bring faults with certain orientations closer to the failure envelope even without an elevation in the pore pressure values (Yang et al., 2017). Figure 6.2 schematically depicts the triggering mechanisms of potential seismicity caused by fluid injection, and the corresponding Mohr Circle representations are plotted in Figure 6.3. 122 Figure 6.2: Schematics of induced earthquakes mechanisms caused by fluid injection (a) Pore pressure elevation in waste disposal sites (b) Rise in the normal stress caused by a near-by energy production site. Figure 6.3: Mohr Circle diagram of induced earthquakes mechanisms caused by fluid injection shows how the stress state perturbations caused by fluid injection can bring an initially-stable fault system (solid semicircle) to a state of failure (dashed semicircle). (a) Pore pressure elevation in waste disposal sites (b) Rise in the normal stress caused by a near-by energy production site. 123 6.2.3 The Role of Fault Damage Zones Field observations often illustrate that large-scale faults consist of fault core planes surrounded by damaged medium possessing a higher fracture intensity compared to the host rock. Fault damage zones form as a result of localized slip along the fault planes over geologic time and contain a network of subsidiary fractures (Caine et al., 1996) making them highly conductive for fluid flow. Figure 6.4 shows schematics of a geologic setting containing fault damage zones. Figure 6.4: Schematics of fault planes surrounded by highly permeable damage zones. As seen in Section 6.2.2, the seismicity risk associated with subsurface exploitation activities involving the withdrawal or injection of fluids is directly related to the in-situ stress state. Although major faults are normally avoided during site selection for such projects, the stress state acting on the fault planes can still be perturbed by remote operations due to pore pressure diffusion. The effect of pore pressure diffusion can be intensified in presence of fault damage zones, as the local fracture density is significantly higher in such areas than adjacent media; thus, the hydraulic communication between the injection well and the fault is facilitated by the highly conductive medium. 124 Fault damage zones are recognized as areas of high importance. A number of studies in the past have focused on characterizing fault damage zones by modeling (Steacy & Sammis, 1992), local measurements (Johri et al., 2014), or lab observation (Faulkner et al., 2011; Misra et al., 2015). Recently, a number of authors have investigated potential well production and injection-induced activation of a fault system through performing numerical analysis on the hydraulic communication between the fault planes and the wells (Chang & Segall, 2016; Asaithambi & Jha, 2017). However, our understanding of the contribution of fault damage zones in creating fault instability is limited due to the complexity of the multiphysics involved. In this work, we follow our previous continuum damage model approach (Samnejad et al., 2017b) to incorporate the effect of fault damage zones in poroelastic pressure diffusion models and illustrate the amplification of the in-situ stress state perturbations in presence of such damaged areas. 6.3 Governing Physics of Poroelascitiy The injection of fluid into the reservoir for waste disposal purposes is accompanied by both pore pressure and stress state changes inside the geologic settings. We employ a coupled computational model to describe the interactions between the fluid flow and geomechanical deformation which concurrently happen in the subsurface environment. The governing physics pertaining to these coupled processes is described by the linear momentum balance equation for geomechanics and Darcy’s law for fluid flow through porous medium as described below. 125 6.3.1 Geomechanics Under quasi-static conditions, the geomechanics equation can be written as the following equation: ∇.σ +ρg = 0 (6.5) in which σ is the total stress tensor, and ρ is the bulk density. The total stress tensor can be related to the effective stress tensor,σ 0 , using Equation 6.6. σ 0 =σ +αpI (6.6) whereI is the identity matrix. The effective stress tensor,σ 0 , is related to the strain tensor,, using the material constitutive laws. We assume that the normal stresses are positive under tension. 6.3.2 Fluid Flow in Porous Media In presence of geomechanical deformation, Darcy flow of a single-phase compressible fluid in the reservoir can be described as: ( α−φ K s + φ K f ) ∂p ∂t +α ∂ v ∂t +∇.[ k μ (−∇p +ρ f g)]−Q p = 0 (6.7) where α is the Biot coefficient, φ is porosity, K s and K f are the bulk moduli for the solid and fluid phases respectively,p is pressure,t is time, v is the volumetric strain,k is the rock permeability tensor, μ is the dynamic fluid viscosity, ρ f is the fluid density, g is the gravitational acceleration vector, and Q p is the source term representing fluid injection. 126 In the following sections, we conduct a number of case studies by building a coupled poroelastic model of the reservoir and fault systems to illustrate how damage zones can alter the pressurization patterns in the subsurface, which can potentially introduce additional risk of induced seismicity. 6.4 Coupled Poroelastic Model of Fluid Injection Our coupled fluid flow and geomechanical model has a 2D domain with dimensions of 2km by 1km in the horizontal (x) and vertical (y) directions respectively as shown in Figure 6.5. Figure 6.5: Model Setup. We include a faulted flat reservoir layer with 100m of thickness inside impermeable rock. The reservoir medium is assumed to be homogeneous and isotropic except for the damagezonewhichislocatedadjacenttothefaultplanes, whichhasdifferentproperties thanthehostrock. Weinitializethemodelundernormalfaultstressregime,andsubject it to an overburden stress of 18 MPa. The lateral boundary conditions of roller type allow for vertical displacement on the left and right boundaries. The bottom of the 127 domain is fixed with a zero displacement (Dirichlet) boundary condition. The damage zone is assumed to be extending over a span of 200m from each side in the base model. The injection well is placed at the far-left boundary of the reservoir. Initially, the model is assumed to be fully saturated and at a hydrostatic equilibrium pressure of 22 MPa. Tables 6.1 and 6.2 summarize the domain parameters and fluid properties used for our coupled model respectively. Table 6.1: Simulation Domain Parameters. Property Nomenclature Value Reservoir Rock Density ρ r 2600 kg/m 3 Reservoir Porosity φ 0.1 Reservoir Permeability κ 25 md Biot Coefficient α 0.9 Bulk Modulus K 1e10 N/m 2 Shear Modulus G 4e8 N/m 2 Fluid Density ρ f 1000 kg/m 3 Fluid Viscosity μ f 1 cp Damage Zone Density ρ r 2600 kg/m 3 Damage Zone Porosity φ 0.2 Damage Zone Permeability κ 50 md Table 6.2: Fluid Properties. Property Nomenclature Value Density ρ f 1000 kg/m 3 Viscosity μ f 1 mPa.s 128 6.5 Results and Discussion In this section, we use the model described in Sections 6.3 and 6.4 to perform a poroelastic numerical study examining how fluid injection into a faulted reservoir perturbs the pore pressure and stress components resulting in potential induction of earthquakes. The results are shown in Figure 6.6. Figure 6.6: The evolution of a) Pore Pressure, b) Effective Normal Stress on the Fault Plane and c) Coulomb Stress after 5 days of pressurized injection in proximity of a sealing fault system with damage zone. Tobetteranalyzeandillustratetheeffectoffaultdamagezones, weplottheeffective normal and shear stresses at the center of the fault plane with respect to depth and time for a period of 5 days of downhole pressurization. The effect of the presence of fault damage zones on the stress patterns is shown in Figures 6.7 and 6.8 for sealing and nonsealing faults respectively. According to Equation 6.2, the effective normal stress acts as a stabilizing force for the fault system. 129 Figure 6.7: Induced seismicity potential for sealing faults with and without damage zones: a) The effective normal stress plotted against depth indicates more severe fault destabilization in presence of damage zones. b) Fault stress path is accelerated towards the failure line when damage zones exist. c) Damage zones cause a quicker drop of the effective normal stress, introducing safety concerns. Therefore, a drop in this stress component leads to more chance of induced seismicity. Thiscanbeinducedthroughtheinjectionoffluidsintotheformation, whichdestabilizes the fault system by reducing the magnitude of the effective normal stress. It can be seen in Figures 6.7-a and 6.8-a that when the damage zone exists, the fault-destabilizing pressure effects are more pronounced compared to a situation where no damage zone is considered to exist around the fault planes. It can also be seen that nonsealing faults act in a safer manner by allowing the fluid to flow through them and preventing severe pressurization. Equation 6.2 also establishes the relationship between the effective shear and normal stress components both at equilibrium stable conditions and at failure. Figure 6.7-b compares the temporal stress paths acting on sealing fault planes with and without 130 Figure 6.8: Induced seismicity potential for nonsealing faults with and without damage zones: a) More pronounced fault destabilization with fault damage zones b) Acceler- ation of fault stress path towards the failure line when damage zones exist c) Unsafe drop of the effective normal stress in presence of fault damage zones. damage zones. One may observe that the initially-stable stress conditions starts to move towards the failure line once the fluid injection begins, and the process of fault destabilization is accelerated in presence of fault damage zones. Again, Figure 6.8-b indicates that this phenomena has a more attenuated form when the fault is nonsealing and allows a certain degree of pressure relief. Another interesting observation, which is essential to prescribing safe injection con- ditions, is the relationship between the changes of the effective normal stress with time as fluid is being injected in nearby wells. Figures 6.7-c and 6.8-c illustrate the process of the reduction of the fault stabilizing force, i.e. the effective normal stress, with time for situations where the area adjacent to the fault is damaged and where no damage zone exists. It should be noted that the presence of fault damage zones creates a more severe drop in the effective normal stress. This implies that injection wells need to shut 131 off sooner than planned for safety purposes when such zones are in place. Both sealing and nonsealing faults respond similarly to the presence of fault damage zones, with nonsealing faults exhibiting a less pronounced pressurization in general due to their higher transmissibility. We conduct three other case studies to investigate the effect of three important parameters on the potential for induced seismicity caused by fault damage zones, namely, the amount of irregularity, the extent of the damage zone, and fault sealing. 6.5.1 The effect of hydromechanical irregularity in fault dam- age zones Fault damage zones are created during geologic time as a result of persistent strain and rupture on and around the fault. The amount of irregularities created depends on the dynamics of the fault system and can vary from fault to fault. Figures 6.9 and 6.10 compare the potential for induced seismicity for sealing and nonsealing faults respectively with fault damage zones having low and high effects translatingtolow(30md)andhigh(100md)permeabilityinthenumericalsimulations. The degree of the anomaly observed in the effective normal stress variation with depth depends on the amount of irregularity which exists in the damage zone as com- pared to the host rock. When such irregularities are high, significant destabilization is expected whereas little to no effective normal stress reduction occurs when the anoma- lies are small, Figures 6.9-a and 6.10-a. The shifting of the stress path towards the failure line, which is depicted in Fig- ure 6.9-b and 6.10-b, exhibits a similar trend for both sealing and nonsealing faults, suggesting that the presence of fault damage zones can, in some cases, be more influ- ential than the transmissibility of the fault itself. Furthermore, for a low-anomaly 132 damage zone, such shifting does not accelerate after a period of injection. Conversely, when such anomalies are high, regardless of fault transmissibility, the stress continues to move towards the dangerous envelope of failure as the injection proceeds. From Figures 6.9-c and 6.10-c, one may infer that the presence of highly damaged zones around the fault systems leads to significantly faster reduction of the effective normal stress, as the fluid is provided with highly conductive paths towards the fault system resulting in rapid pressurization regardless of the fault being sealing or nonseal- ing. Figure 6.9: Induced seismicity potential for sealing faults in presence of low and highly damaged zones: a) The effective normal stress plotted against depth indicates more severe fault destabilization in presence of highly damaged zones. b) The shifting of the stress path towards the failure line accelerates with the existence of highly damaged zones. c) Rapid drop of the effective normal stress, introducing safety concerns around highly damaged fault zones. 133 Figure 6.10: Induced seismicity potential for nonsealing faults in presence of low and highly damaged zones: a) More pronounced fault destabilization with highly damaged zones even when the fault is nonsealing b) Acceleration of fault stress path towards the failure line when damages are high c) Quick and unsafe drop of the effective normal stress in presence of highly damaged zones. 6.5.2 The effect of the geometry of fault damage zones Fault damage zones can vary in size in accordance with the geologic setting they are created by. In our model, we change the extent of the fault damage zone from 200 m to 400 m to investigate how the extent of such zones affects the perturbations created in the in-situ stress state. Figures 6.11-a and 6.12-a imply that a larger fault damage zone results in a higher reduction in the effective normal stress. A large damage zone is seen to allow for rapid pressure diffusion from the injection well towards the fault shifting the stress path towards the failure line compared to the sealing faults. Similar to previous scenarios, higher fault transmissibility acts as a pressure-releasing factor, Figures 6.11-b and 6.12-b. This phenomenon is also confirmed in Figures 6.11-c and 6.12-c. 134 Figure 6.11: Induced seismicity potential for sealing faults in presence of small and large damaged zones: a) The effective normal stress plotted against depth indicates more severe fault destabilization in presence of larger damaged zones. b) The shifting of the stress path towards the failure line accelerates with the existence of larger damaged zones. c) Rapid drop of the effective normal stress, introducing safety concerns around large damaged fault zones. 6.5.3 The effect of fault transmissibility Fault transmissibility is among the factors which can affect the pressure diffusion patterns inside the formation. Sealing faults, which have zero to low transmissibility, causemorepressurizationbycreatingabarrierforthefluid. Onthecontrary, nonsealing faults allow for the fluid to flow through them, which releases a part of the pore pressure increase. It can be observed in Figures 6.13 through 6.14 that while such effects are not distinctively visible when no damage zone exists, the behavior of sealing and nonsealing faults becomes more distinguishable in presence of fault damage zones. Figure 6.13-a showsverysimilareffectivenormalstressbehaviorforbothsealingandnonsealingfaults when no damage zone exists. Plots in Figure 6.13-b also indicate very similar stress 135 Figure 6.12: Induced seismicity potential for nonsealing faults in presence of small and large damaged zones: a) More pronounced fault destabilization with larger damaged zonesb)Accelerationoffaultstresspathtowardsthefailurelinewhendamagesarelarge c) Quick and unsafe drop of the effective normal stress in presence of large damaged zones. patterns for sealing and nonsealing faults, which is also confirmed in Figure 6.13-c. This emphasizes the importance of taking fault damage zones into account in addition to fault transmissibility, as pore pressure perturbations, and consequently, in-situ stress patterns, can differ significantly for sealing and nonsealing fault systems when such zones are in place, Figure 6.14-a. The shift of the in-situ stress components towards the failure line plotted in Figure 6.14-b also indicates a different behavior of sealing and nonsealing faults, with sealing faults exhibiting more failure tendency due to more severe pressurization, which is also confirmed in Figure 6.14-c. 136 Figure 6.13: Induced seismicity potential for sealing and nonsealing faults without damaged zones: a) The effective normal stress plotted against depth indicates very similar stress drops for both systems when damage zones are not present. b) The shiftingofthestresspathtowardsthefailurelineshowsnoobservabledifferencebetween sealingandnonsealingfaults. c)Theeffectivenormalstressdropsinasomewhatsimilar fashion for both sealing and nonsealing. 6.6 Conclusion The injection of fluids into subsurface formations create pore pressure perturba- tions, which is transmitted to faults in forms of shifts in the in-situ stress components destabilizing the fault system inducing a higher potential for seismicity. Predicting such stress perturbations, which is crucial for ensuring safe injection pressure limits, can be a challenging task, partially, due to the geologic complexity of the subsurface environment. We developed a poroelastic model of fluid injection scenarios taking into account the effect of fault damage zones, which have proven to play an important role in controlling the pressure diffusion patterns inside target formations. We analyze the spa- tial and temporal changes of the important Coulomb stress components acting on the 137 Figure 6.14: Induced seismicity potential for sealing and nonsealing faults in presence of damage zones: a) More severe fault destabilization for sealing faults b) More accel- eration of the shifting of the stress path towards the failure line for sealing faults c) Quick and unsafe drop of the effective normal stress when more pressurization created by sealing faults. fault planes and demonstrate the importance of taking fault damage zones by studying the effect of three main factors affecting the perturbed stress patterns. The include the magnitude and extent of damage anomalies as well as fault transmissibility. The results of our analysis illustrate that it is important to characterize and model the degree of damage to the rock in proximity of the fault systems as it has significant impact on the shifting of the effective stress components. The extent of fault damage zones also influences the pressure diffusion patterns inside the reservoir, which can potentially lead to more severe shifting of the stress components towards the Coulomb failure line when such zones extend to a large area around the faults. Lastly, fault transmissibility can partially alleviate the pressurization by allowing the fluid to flow to another compart- ment of the reservoir. However, one may need to examine the geologic setting before 138 injecting into compound formations in order to avoid undesired flow of fluid into risky areas for the case of stacked reservoirs. 139 Chapter 7 Summary, Conclusions, and Future Work Knowledgeofthecomplexphysicsofporousrockfracturingandsimultaneousflowof fluid is critical to a variety of subsurface applications, including petroleum production, nuclear waste injection, and CO 2 sequestration. Despite their significance, the available approaches in this domain are oversimplified due to a variety of challenges associated with the modeling of such processes. The lack of proper understanding of the related physics and associated complexities often times led to unrealistic predictions. This, in turn,cangiverisetosub-optimalfielddevelopment,and/orpotentialinducedseismicity, or other environmental hazards. We present a robust computational platform for simulating these phenomena by pursuing an effective continuum method. Unlike current methods, our model has the capability to demonstrate the evolving propagation path of created fractures and their growth within reasonable computation time while also providing with an estimate of the resulting permeability enhancement during well stimulation operations. Unlike current models, we take into account the effect of various degrees of ruck ductility which may be observed when hydraulic fracturing jobs are in progress. We also demonstrate how the current methods fail to account for a number of important parameters when estimating the outcome of hydraulic fracturing jobs. We also illustrate how our comprehensive model can be deployed to quantify the effect of those overlooked impact factors. 140 We demonstrate the validity of our approach by calibrating our model with the availableexperimentaldatawhichcapturesthebehavioroftheporousmediaundergoing the stimulation processes both in terms of mechanical deformation & failure and fluid flow enhancement. The results obtained using our approach shows the superiority of our model in providing with more accurate predictions of the economic success of well stimulation investments as well as the required operational specifications which need to be under- taken in order to attain the desired results. When implemented at the field scale, we show that our model is not only appli- cable in modeling hydraulic fracturing processes, but it also can help avoid induced earthquakes by more accurately predicting the changes of stresses acting on active fault systems. This is another critical application of our methodology in ensuring safety in the design of waste disposal sites, which helps avoiding the risk of potential induced seismicity. It provides a broader insight about the stresses acting on fault systems. A summary of our contributions is shown in Figure 7.1. Figure 7.1: Summary of contributions and comparison with previous work. 141 We have identified a number of areas where there is opportunity for future advances beyond the current work. Some of such future research topics are summarized below. • This dissertation incorporates the complexities imposed by the in-situ stress anisotropy on the fracture propagation geometry. One may also consider the structural anisotropy created by sedimentary strata in geologic formations using a non-isotropic stiffness tensor, e.g. transverse isotropic material. • This work provides a good understanding of the propagation pattern of fractures, in terms of length and direction in response to downhole injection operations. The ultimate recovery outcome, however, also depends on the type of the proppants used (Aminzadeh, 2019) to keep fractures open as well as certain production settings. The investigation of the physics of post-stimulation fracture closure could be a valuable addition to the current body of knowledge. • Throughout this work, we have focused on rock fracturing processes coupled with geomechanical deformation and single-phase flow of fluid. Another layer of com- plexity during these phenomena can be added through the consideration of the presence of multi-phase fluid flow (Sherafati & Jessen, 2018). • We have calibrated our model quantitatively with available experimental data and qualitatively with observed field data. Further calibration at the field scale is possible through linking of mechanical damage processes with seismic velocities (Hamiel et al., 2009), which is outlined in the literature and can be exhibited in field microseismic measurements. • Last but not least, the digital transformation occurring in the oil and gas indus- try urges scholars to integrate data science and machine learning (ML) techniques with their legacy research methods (Shirangi et al., 2020; Popp & Shirangi, 2020). 142 Specifically, data-driven methods can be employed in conjunction with computa- tional simulations for smart sweet spot identification (Maity & Aminzadeh, 2015; Aminzadeh, 2019), development of proxy models (Alenezi & Mohaghegh, 2016), rig operations monitoring and control (Samnejad et al., 2020), field-scale develop- ment (Shirangi, 2019), and optimum design prescription (Shirangi et al., 2019). Such applications of AI-assisted physics modeling are summarized schematically in Figure 7.2 Figure 7.2: Benefits of AI-assisted physics modeling. 143 Reference List Acumen Research and Consulting (2019) Hydraulic fracturing market size hit $83 billion by 2026. Alenezi F, Mohaghegh S (2016) A data-driven smart proxy model for a comprehensive reservoir simulation In 2016 4th Saudi International Conference on Information Technology (Big Data Analysis) (KACSTIT), pp. 1–6. Aminzadeh F (2019) Hydraulic Fracturing and Well Stimulation Wiley. Aminzadeh F (2018) Hydraulic fracturing, an overview. Journal of Sustainable Energy Engineering 6:204–228. Aminzadeh F (2019) Hydraulic Fracturing, pp. 1–17 Springer New York, New York, NY. Armero F (1999) Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions. Computer Methods in Applied Mechanics and Engineering 171:205 – 241. Asaithambi R, Jha B (2017) Coupled Flow and Geomechanical Modeling of Reser- voir Seismicity: Effects of Hydraulic Communication and Well Type. SPE Western Regional Meeting, 23-27 April, Bakersfield, CA . Ashby M, Sammis C (1990) The damage mechanics of brittle solids in compression. Pure and Applied Geophysics 133:489–521. Bao X, Eaton DW (2016) Fault activation by hydraulic fracturing in western canada. Science 354:1406–1409. BartonN,BandisS,BakhtarK(1985) Strength, deformationandconductivitycoupling of rock joints In International journal of rock mechanics and mining sciences & geomechanics abstracts, Vol. 22, pp. 121–140. Elsevier. Bear J (1988) Dynamics of Fluids in Porous Media Dover Civil and Mechanical Engi- neering Series. Dover. 144 Biot MA (1941) General theory of three-dimensional consolidation. Journal of Applied Physics 12:155–164. Brown SP, Huntington HG (2010) Estimating us oil security premiums. Resources for the Future, Washington, DC . Brown SP, Yücel MK (2013) Shale gas and tight oil boom: Us states economic gains and vulnerabilities. Caine JS, Evans JP, Forster CB (1996) Fault zone architecture and permeability struc- ture. Geology 24:1025–1028. Carman PC (1937) Fluid flow through granular beds. Trans. Inst. Chem. Eng. 15:150–166. Chang KW, Segall P (2016) Seismicity on basement faults induced by simultaneous fluid injection–extraction. Pure and Applied Geophysics 173:2621–2636. Chen W, Han D (2007) Plasticity for Structural Engineers J. Ross Publishing Classics. J. Ross Pub. Chow CL, Wang J (1989) Crack propagation in mixed-mode ductile fracture with con- tinuum damage mechanics. Proceedings of the Institution of Mechanical Engineers, Part C: Mechanical Engineering Science 203:189–199. Costin LS (1983) A microcrack model for the deformation and failure of brittle rock. Journal of Geophysical Research: Solid Earth 88:9485–9492. Coulomb CA (1776) Essai Sur Une Application Des regles De Maximis et Minimis a Quelques problemes De Statique, Relatifs a L’architecture Paris : De l’Imprimerie Royale. Coussy O (1995) Mechanics of Porous Continua Wiley. Cundall PA, Hart RD (1992) Numerical modelling of discontinua. Engineering compu- tations 9:101–113. Curtis JB (2002) Fractured shale-gas systems. AAPG bulletin 86:1921–1938. de Borst R, Gutiérrez MA, Wells GN, Remmers JJ, Askes H (2004) Cohesive-zone mod- els, higher-order continuum theories and reliability methods for computational failure analysis. International Journal for Numerical Methods in Engineering 60:289–315. Duarte Azevedo I, Vaz L, Vargas E (1998) A numerical procedure for the analysis of the hydromechanical coupling in fractured rock masses. International journal for numerical and analytical methods in geomechanics 22:867–901. 145 Faulkner DR, Mitchell TM, Jensen E, Cembrano J (2011) Scaling of fault damage zones with displacement and the implications for fault growth processes. Journal of Geophysical Research: Solid Earth 116. Ferronato M, Gambolati G, Janna C, Teatini P (2008) Numerical modelling of regional faults in land subsidence prediction above gas/oil reservoirs. International journal for numerical and analytical methods in geomechanics 32:633–657. Foulger GR, Wilson MP, Gluyas JG, Julian BR, Davies RJ (2018) Global review of human-induced earthquakes. Earth-Science Reviews 178:438 – 514. Fumagalli A, Scotti A (2013) A numerical method for two-phase flow in fractured porous media with non-matching grids. Advances in Water Resources 62:454–464. Gadde PB, Sharma MM et al. (2001) Growing injection well fractures and their impact on waterflood performance In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Geertsma J, De Klerk F et al. (1969) A rapid method of predicting width and extent of hydraulically induced fractures. Journal of Petroleum Technology 21:1–571. Gens A, Carol I, Alonso E (1989) An interface element formulation for the analysis of soil-reinforcement interaction. Computers and Geotechnics 7:133–151. Goebel THW, Hosseini SM, Cappa F, Hauksson E, Ampuero JP, Aminzadeh F, Saleeby JB (2016) Wastewater disposal and earthquake swarm activity at the southern end of the central valley, california. Geophysical Research Letters 43:1092–1099. Gowd T, Rummel F (1980) Effect of confining pressure on the fracture behaviour of a porous rock In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, Vol. 17, pp. 225–229. Elsevier. Grand View Research, Inc. (2016) Hydraulic fracturing market analysis by technology (plug & perf, sliding sleeve), material, application (shale gas, tight gas, tight oil, coal bed methane (cbm)) and segment forecasts to 2024. Griffith AA (1921) Vi. the phenomena of rupture and flow in solids. Philosophical trans- actions of the royal society of london. Series A, containing papers of a mathematical or physical character 221:163–198. Hamiel Y, Lyakhovsky V, Stanchits S, Dresen G, Ben-Zion Y (2009) Brittle defor- mation and damage-induced seismic wave anisotropy in rocks. Geophysical Journal International 178:901–909. 146 Hosseini SM, Goebel THW, Jha B, Aminzadeh F (2018) A probabilistic approach to injection-induced seismicity assessment in the presence and absence of flow bound- aries. Geophysical Research Letters 45:8182–8189. Hough S, Tsai V, Walker R, Aminzadeh F (2017) Was the mw 7.5 1952 kern county, california, earthquake induced (or triggered)? Journal of Seismology . Hough SE, Page MT (2015) A century of induced earthquakes in Oklahoma? Bulletin of the Seismological Society of America 105. Inglis CE (1913) Stresses in a plate due to the presence of cracks and sharp corners. Trans Inst Naval Archit 55:219–241. J van der Elst N, M Savage H, Keranen K, A Abers G (2013) Enhanced remote earthquake triggering at fluid-injection sites in the midwestern united states. Sci- ence 341:164–167. JabbariN, AminzadehF, BarrosFPJ(2017) Hydraulicfracturing andtheenvironment: risk assessment for groundwater contamination from well casing failure. Stochastic Environmental Research and Risk Assessment 31:1527 – 1542. Jahandideh A, Jafarpour B (2016) Optimization of hydraulic fracturing design under spatially variable shale fracability. Journal of Petroleum Science and Engineer- ing 138:174–188. Jalbout A (2006) Etude experimentale de l’influence de l’endommagement induit sur le comportement mecanique d’une roche fragile dissertation, Ecole Polytechnique Universitaire de Lille, Lille, France. Jha B (2016) Evolution of stress transfer mechanisms during mechanical interaction between hydraulic fractures and natural fractures. Journal of Sustainable Energy Engineering 4:296–309. Jha B, Juanes R (2014) Coupled multiphase flow and poromechanics: A computa- tional model of pore pressure effects on fault slip and earthquake triggering. Water Resources Research 50:3776–3808. Johri M, Zoback MD, Hennings P (2014) A scaling law to characterize fault-damage zones at reservoir depths In AAPG Bulletin. American Association of Petroleum Geologists. Kachanov L (1958) Time of the Rupture Process under Creep Conditions. Izvestiia Akademii Nauk SSSR, Otdelenie Teckhnicheskikh Nauk 8:26 – 31. 147 Kachanov L (1986) Introduction to continuum damage mechanics J. Ross Publishing Classics. Springer Netherlands. Khodabakhshnejad A (2016) An extended finite element method based modeling of hydraulic fracturing Ph.D. diss., Viterbi School of Engineering, Los Angeles, CA. Khoei A, Hosseini N, Mohammadnejad T (2016) Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model. Advances in water resources 94:510–528. Kim J, Tchelepi HA, Juanes R et al. (2009) Stability, accuracy and efficiency of sequen- tial methods for coupled flow and geomechanics In SPE reservoir simulation sympo- sium. Society of Petroleum Engineers. Kozeny J (1927) Uber kapillare leitung des wassers im boden, vol. 136 (2a). Akademie der Wissesschaften, Wien p. 271. Langenbruch C, Zoback MD (2016) How will induced seismicity in oklahoma respond to decreased saltwater injection rates? Science Advances 2. Lemaitre J (1996) A course on damage mechanics Springer, Berlin, Heidelberg. Lu Y, Elsworth D, Wang L (2013) Microcrack-based coupled damage and flow mod- eling of fracturing evolution in permeable brittle rocks. Computers and Geotech- nics 49:226 – 244. Ma J (2015) Review of permeability evolution model for fractured porous media. Jour- nal of Rock Mechanics and Geotechnical Engineering 7:351 – 357. Maity D, Aminzadeh F (2015) Novel fracture zone identifier attribute using geophysical and well log data for unconventional reservoirs. Interpretation 3:T155–T167. Mase G (1970) Schaum’s Outline of Theory and Problems of Continuum Mechanics [including 360 Solved Problems] McGraw Hill, New York, NY. Maxwell SC, Chorney D, Goodfellow SD (2016) Microseismic geomechanics of hydraulic-fracture networks: Insights into mechanisms of microseismic sources. The Leading Edge . McClureMW,BabazadehM,ShiozawaS,HuangJetal.(2016) Fullycoupledhydrome- chanical simulation of hydraulic fracturing in 3d discrete-fracture networks. SPE Journal 21:1–302. Misra S, Ellis S, Mandal N (2015) Fault damage zones in mechanically layered rocks: The effects of planar anisotropy. Journal of Geophysical Research: Solid Earth 120:5432–5452. 148 Mohammadnejad T, Khoei A (2013) Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method. International Journal for Numerical and Analytical Methods in Geome- chanics 37:1247–1279. Nandakumar J (2007) Sino-Indian cooperation in the search for overseas petroleum resources: Prospects and implications for india. International Journal of Energy Sector Management 1:84–95. National Research Council (2013) Induced Seismicity Potential in Energy Technologies The National Academies Press, Washington, DC. Nola IA, Jelinić JD, Žuškin E, Kratohvil M (2013) Earthquakes – a historical review, environmental and health effects, and health care measures. Archives of Industrial Hygiene and Toxicology 64:327 – 337. Nordgren R et al. (1972) Propagation of a vertical hydraulic fracture. Society of Petroleum Engineers Journal 12:306–314. PapanastasiouP(1997) Theinfluenceofplasticityinhydraulicfracturing. International Journal of Fracture 84:61–79. Parker JC (1989) Multiphase flow and transport in porous media. Reviews of Geo- physics 27:311–328. Peach CJ, Spiers CJ (1996) Influence of crystal plastic deformation on dilatancy and permeability development in synthetic salt rock. Tectonophysics 256:101–128. Pickenpaugh GC, Adder JM (2018) Shale gas production and labor market trends in the u.s. marcellus-utica region over the last decade. Pogacnik J, Elsworth D, O’Sullivan M, O’Sullivan J (2016) A damage mechanics approach to the simulation of hydraulic fracturing/shearing around a geothermal injection well. Computers and Geotechnics 71:338 – 351. Popp T, Shirangi MG (2020) Prescriptive data analytics to optimize casing exits In IADC/SPE International Drilling Conference and Exhibition (SPE-199612-MS), Houston, Texas. IADC, SPE. Popp T, Kern H, Schulze O (2001) Evolution of dilatancy and permeability in rock salt during hydrostatic compaction and triaxial deformation. Journal of Geophysical Research: Solid Earth 106:4061–4078. Porter ME, Gee DS, Pope GJ (2016) America’s unconventional energy opportunity. 149 Pruess K (2005) Numerical studies of fluid leakage from a geologic disposal reservoir for co2 show self-limiting feedback between fluid flow and heat transfer. Geophysical research letters 32. Rabotnov YN (1969) Creep rupture In Applied mechanics, pp. 342–349. Springer. Rice JR, Cleary MP (1976) Some basic stress diffusion solutions for fluid-saturated elas- tic porous media with compressible constituents. Reviews of Geophysics 14:227–241. Rubinstein J (2015) Increasing rate of earthquakes beginning in 2009. Rutter E, Glover C (2012) The deformation of porous sandstones; are byerlee friction and the critical state line equivalent? Journal of Structural Geology 44:129 – 140. Samnejad M, Aminzadeh F, An Y (2014) The investigation of an integrative approach for fracture density inversion using azavo analysis In Pacific Section AAPG, SPE and SEPM Joint Technical Conference, Bakersfield, California. AAPG, SPE, and SEPM. Samnejad M, Aminzadeh F, Bubshait A (2017a) Acoustoelastic Estimation of the In Situ Stresses from Sonic Logs for a Carbonate Reservoir In SPE Western Regional Meeting (SPE-185644-MS), Bakersfield, California. Society of Petroleum Engineers. Samnejad M, Aminzadeh F, Jha B (2017b) Simulation of hydraulic fracturing-induced permeability stimulation using coupled flow and continuum damage mechanics In SPE Annual Technical Conference and Exhibition (SPE-187250-MS), San Antonio, Texas. Society of Petroleum Engineers. Samnejad M, Aminzadeh F, Jha B (2018) A novel approach for studying hydraulic fracturing success factors beyond brittleness indices In 52nd US Rock Mechan- ics/Geomechanics Symposium (ARMA 18-0187), Seattle, Washington. American Rock Mechanics Association. Samnejad M, Shirangi MG, Ettehadi Osgouei R (2020) A digital twin of drilling fluids rheologyforreal-timerigoperations InOffshore Technology Conference (OTC-30738- MS), Houston, Texas. SPE. Samnejad M, Walker R, Aminzadeh F (2018) Understanding injection-induced seismic- ity effects of fault damage zones: Beyond poroelastic models In AGU Fall Meeting (abstract S21A-07), Washington, D.C. American Geophysical Union. Shao JF, Zhou H, Chau KT (2005) Coupling between anisotropic damage and perme- ability variation in brittle rocks. International Journal for Numerical and Analytical Methods in Geomechanics 29:1231–1247. 150 Shao JF, Lu Y, Lydzba D (2004) Damage modeling of saturated rocks in drained and undrained conditions. Journal of engineering mechanics 130:733–740. Shen X, Standifird W, Evans S et al. (2016) Optimizing multistage hydraulic-fracturing designbasedon3dcontinuumdamagemechanics InSPE/IADC Middle East Drilling Technology Conference and Exhibition. Society of Petroleum Engineers. Sherafati M, Jessen K (2018) Coarse-Scale Modeling of Multicomponent Diffusive Mass Transfer in Dual-Porosity Models In SPE Annual Technical Conference and Exhibi- tion, Dallas, Texas. Society of Petroleum Engineers. Shirangi MG DLJ (2016) A general method to select representative models for decision making and optimization under uncertainty. Computers and Geo- sciences 96:109 – 123. Shirangi MG (2014) History matching production data and uncertainty assessment with an efficient tsvd parameterization algorithm. Journal of Petroleum Science and Engineering 113:54–71. Shirangi MG (2019) Closed-loop field development with multipoint geostatistics and statistical performance assessment. Journal of Computational Physics 390:249–264. Shirangi MG, Ettehadi Osgouei R, Aragall R, Furlong E, May R, Dahl T, Samnejad M, Thompson C (2020) Digital twins for drilling fluids: Advances and opportunities In IADC/SPE International Drilling Conference and Exhibition (SPE-199681-MS), Houston, Texas. IADC, SPE. Shirangi MG, Wilson T, Furlong E et al. (2019) Prescriptive analytics for completion optimization in unconventional resources In SPE Western Regional Meeting. SPE- 195311-MS. Shojaei A, Taleghani AD, Li G (2014) A continuum damage failure model for hydraulic fracturing of porous rocks. International Journal of Plasticity 59:199–212. Simo J (1992) Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Computer Methods in Applied Mechanics and Engineering 99:61 – 112. Steacy SJ, Sammis CG (1992) A damage mechanics model for fault zone friction. Journal of Geophysical Research: Solid Earth 97:587–594. Tate RLIII (2005) Encyclopedia of Soils in the Environment: Volume 1-4, Vol. 170 Elsevier. 151 Terzaghi K, Peck RB, Mesri G (1996) Soil mechanics in engineering practice John Wiley & Sons. Tiwari A, Walker R, Aminzadeh F (2014) Hydraulic fracturing and induced seismicity in california. Hydraulic Fracturing Journal . U. S. Energy Information Administration (2019) Monthly energy review. Walker R, Samnejad M, Aminzadeh F (2018) Consideration of stress induced per- meability changes within a poroelastic model framework for induced and triggered seismicity In AGU Fall Meeting (abstract S33C-0592), Washington, D.C. American Geophysical Union. Walsh JB, Brace W (1984) The effect of pressure on porosity and the transport prop- erties of rock. Journal of Geophysical Research: Solid Earth 89:9425–9431. Wang JA, Park H (2002) Fluid permeability of sedimentary rocks in a complete stress– strain process. Engineering geology 63:291–300. Wasantha P, Konietzky H (2017) Hydraulic fracture propagation under varying in-situ stress conditions of reservoirs. Procedia Engineering 191:410 – 418 ISRM European Rock Mechanics Symposium EUROCK 2017. Yang H, Liu Y, Wei M, Zhuang J, Zhou S (2017) Induced earthquakes in the develop- mentofunconventionalenergyresources. Science ChinaEarth Sciences60:1632–1644. Yo M, Neff S (14 July 2014) Liquid fuels and natural gas in the americas. Zangeneh N, Eberhardt E, Bustin R (2015) Investigation of the influence of natural fractures and in situ stress on hydraulic fracture propagation using a distinct-element approach. Canadian Geotechnical Journal 52:926–946. Zhao Q, Lisjak A, Mahabadi O, Liu Q, Grasselli G (2014) Numerical simulation of hydraulic fracturing and associated microseismicity using finite-discrete element method. Journal of Rock Mechanics and Geotechnical Engineering 6:574–581. Zhou F, Molinari JF (2004) Dynamic crack propagation with cohesive elements: a methodologytoaddressmeshdependency. International Journal for Numerical Meth- ods in Engineering 59:1–24. Zhu W, Tang C (2004) Micromechanical model for simulating the fracture process of rock. Rock Mechanics and Rock Engineering 37:25–56. Zimmermann G, Reinicke A (2010) Hydraulic stimulation of a deep sandstone reser- voir to develop an enhanced geothermal system: Laboratory and field experiments. Geothermics 39:70–77. 152 Zoback MD (2007) Reservoir Geomechanics Cambridge University Press. 153

Abstract (if available)

Abstract In this dissertation, we focus on hydromechanical evolution of changes in subsurface reservoir systems during pressurized injection. Our main contribution is to extend the applicability of the available Effective Continuum Model (ECM) methods, previously used for modeling enhanced geothermal systems, to the domains of hydraulic fracturing and waste water injection by addressing the major deficiencies of the current models including: ❧ Multi-mechanism Fractures: The prevailing ECM approach for hydraulic fracturing simulation is based on Linear Elastic Fracture Mechanics (LEFM). Generally, LEFM gives reasonable hydraulic fracturing results predictions for hard brittle rocks, but fails to predict the fracture propagation path and pressure patterns in formations which can undergo ductile failure, such as poorly consolidated and unconsolidated sands as well as ductile shales. There are very few previous attempts to model inelastic mechanical behavior during pressurized injection, all lacking a corresponding fluid flow response to plastic deformation. Lack of a proper ductile model leads to overestimation of the success of a hydraulic fracturing job. ❧ We address this shortcoming by adopting an elastoplastic material model for the porous medium, which addresses the multi-mechanism fracture propagation in the rock along both the linear elastic and nonlinear inelastic deformation paths. Crack growth in the medium is then translated into the development of continuum damage, which is related to the loss of solid surface inside each volumetric element creating a better conduction capability for fluid flow, i.e. permeability enhancement. ❧ Anisotropy Consideration: Previous ECM approaches only consider isotropic conditions, which is a reasonable assumption to make for geothermal systems made of the structurally-isotropic granite. Fractures formed in such systems often follow a sharp path away from the wellbore in the direction perpendicular to the minimum principal stress. It has been observed, however, that sedimentary rocks can exhibit more complex forms of fracture propagation patterns due to the pre-existing planes of weakness in their structure. Lack of accurate predictions on fracture propagation patterns results in sub-optimal well placement and trajectory design. ❧ We model the tendency of fractures to grow in preferred directions by allocating a set of pre-existing planar flaws placed in different orientations to represent a form of secondary anisotropy inside an isotropic medium. Our model incorporates the effects of both stress asymmetry and structural anisotropy, demonstrating a better match with field microseimic measurements and lab observations. ❧ In summary, we exhibit the superiority of our hydraulic fracturing model by fitting it into the available coupled fluid flow and geomechanics platforms which allows to capture the phenomena occurring at a microscale in a macroscopic framework. Unlike other common approaches, our method does not require prior specification of the crack path, has no mesh dependency, and is capable of incorporating heterogeneity while maintaining a reasonable computational time. We also identify areas of opportunity for future research and propose a scope for extending beyond the current work.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Reactivation of multiple faults in oilfields with injection and production

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

PDF

An extended finite element method based modeling of hydraulic fracturing

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

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

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

$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

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

PDF

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

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

Geothermal resource assessment and reservoir modeling with an application to Kavaklidere geothermal field, Turkey

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

PDF

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

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

PDF

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

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale

Asset Metadata

Creator Samnejad, Mahshad (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 02/04/2020

Defense Date 11/08/2019

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

Tag anisotropy,computational reservoir simulation,continuum damage mechanics,coupled geomechanics,energy sustainability,faults,fluid flow through porous media,hydraulic fracturing,induced seismicity,material degradation,national economy,OAI-PMH Harvest,permeability enhancement,predictive modeling,pressurized injection,Public welfare,rock brittleness,subsurface multiphysics,well stimulation

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Aminzadeh, Fred (committee chair), Ershaghi, Iraj (committee member), Nakano, Aiichiro (committee member)

Creator Email mahshadsamnejad@gmail.com,samnejad@usc.edu

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

Unique identifier UC11673846

Identifier etd-SamnejadMa-8139.pdf (filename),usctheses-c89-264266 (legacy record id)

Legacy Identifier etd-SamnejadMa-8139.pdf

Dmrecord 264266

Document Type Dissertation

Rights Samnejad, Mahshad

Type texts

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

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

Repository Name University of Southern California Digital Library

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