University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content A COUPLED POROMECHANICS-DAMAGE MECHANICS APPROACH TO MODEL FRACTURING IN MULTIPHASE FLOW by Ahmed Khaled Bubshait A Dissertation Presented to FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA in Partial Fulllment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY In Petroleum Engineering Ph.D. Committee Prof. B. Jha (Advisor), Prof. I. Ershaghi, and Prof. R. Ghanem December 2019 Copyright 2019 Ahmed Bubshait Dedication This dissertation is dedicated to my parents, Khaled and Huda, for their unlimited love, dedication and support in every day of my life, and to my lovely wife Faten for her love, compassion, encouragement and patience over the many years we have been together, to my daughters Aljoud, Deema and Nouf and son Khaled for bringing the greatest joy to my life, and to my parents-in-law, Khalid and Wafa, for their sincere support and encouragement. i Acknowledgements Primarily, I would like thank Allah for being able to complete this thesis with success. This thesis has been made possible by many people who have helped me with consider- able eort as follow: My deepest appreciation goes to my advisor, Prof. Birendra Jha who I have been working for him since 2016. For three years, I have overcome a lot of diculty in my research with which I would have been lost. He rst motivated me for coupled ow and geomechanics as the Ph.D. research topic. He also encouraged to build my own code where I successfully made the code. Whenever I faced a problem, he always encouraged me with warm advice. He guided me on the right track, showing me valuable and critical suggestions and providing me with the freedom for various ideas. This was a great opportunity and privilege for me to work for Prof. Birendra Jha. My thanks also go to Prof. Fred Aminzadeh , my co-adviser. I have worked for him for four years, when I took the rst step in my PhD program. I would like to express my heartfelt thanks to Prof. Iraj Ershaghi. His classes on advanced reservoir characterization and unconventional reservoir have become critical resources for this thesis. Also, he listened to my research challenges and gave me fruitful advice on naturally fractured formation at a turning point of my research. My sincere thanks go to another member of my Ph.D. committee, Prof. Roger Ghanem, who made valuable comments regarding numerical coupling between dierent processes. I gratefully acknowledge the funding source from Saudi Aramco, that has made my Ph.D. work possible. I also would like to express my deep appreciation and thank the President and CEO of Saudi Aramco, MR. Amin AL-Nasser, Upstream Senior VP , Dr Mohammed Al-Qahtani and MR. Nasir Al Naimi Vice President Petroleum Engineering Development for their support, advise and encouragement. I am very grateful to my manager, MR. Ali Al-Shahri, my superior, MR. Arafat Shurei, and my advisor, MR. Khaled Kilany, who have helped me for the initial study on coupled ow, geomechanics and fracture approach. Needless to say, I am thankful to all of my colleagues at USC and the sta in the Mork Family Department of Chemical Engineering and Materials Science. ii Abstract Injection of water or gas in a naturally fractured low permeability reservoir often leads to dynamic changes in permeability and elastic moduli of the reservoir due to simul- taneous propagation of multiple pre-existing fractures under injection-induced stress. Wellbore injectivity, reservoir productivity, and spatial distribution of uids within the reservoir depend on these properties. Numerical modeling of the evolution in the prop- erties require solution of the coupled physical processes of uid ow, rock deformation, and fracture propagation, which is challenging when multiple fractures are propagating simultaneously due to nonlinearity of the governing equations and presence of multi- ple length scales. We present a numerical framework to model the dynamic evolution in permeability and elastic moduli based on a coupled formulation of multiphase ow, poroelasticity, and fracture-induced damage processes. We solve the uid ow problem using the nite volume method and the quasi-static deformation problem using the - nite element method. We use the continuum damage mechanics method to propagate the fractures and evolve the properties dened at the representative elementary volume scale. We generalize the xed stress method of solving the coupled ow-deformation problem to model fracturing-induced anisotropy in permeability, elastic moduli and in- duced stresses. The proposed framework can model complex fracture networks with randomly distributed fractures and does not require an explicit representation of frac- tures in the computational mesh. This reduces the computational burden associated with modeling fracture paths and ow through discrete fractures. As a result, fracture- specic multiphase ow functions are not required. We use the proposed framework to study the eects of variable injection/production rate, multiple fracture networks, boundary stress anisotropy, and multiphase ow properties on evolution in the proper- ties. Our work shows that the eect of fracture propagation can be included in existing, sequentially coupled ow-geomechanics models of fractured reservoirs relatively easily by implementing a damage-based solver module within the sequential loop. iii Table of Contents Page Dedication ................................................................................. i Acknowledgements ........................................................................ ii Abstract.................................................................................... iii List of Figures ............................................................................. vii List of Tables .............................................................................. xix Nomenclature .............................................................................. xx CHAPTER 1 INTRODUCTION ........................................................................ 1 1.1 Naturally Fractured Reservoirs ................................................. 3 1.1.1 Denition of a Fracture ................................................. 3 1.2 Theory of Poromechanics ....................................................... 4 1.2.1 Stress-Strain Relation and Poroelastic Properties of a Rock......... 5 1.3 Fracture Mechanics ............................................................. 6 1.3.1 Grith Cracks ........................................................... 7 1.4 Fracture Modeling ............................................................... 7 1.4.1 Dual Porosity and Dual Permeability Models ......................... 8 1.4.2 Cohesive Zone Elements ................................................ 9 1.4.3 Extended Finite Element Method (XFEM) ........................... 9 1.4.4 Discrete Fracture Network .............................................. 10 1.4.5 Lorenz Coecient ....................................................... 11 CHAPTER 2 MATHEMATICAL MODELING AND GOVERNING EQUATIONS ............... 13 2.1 Poroelastic Problem ............................................................. 13 2.1.1 Single-phase Poroelasticity ............................................. 14 2.1.2 Multiphase Poroelasticity ............................................... 18 2.2 Fracture Dynamics and Propagation ........................................... 20 2.2.1 Stress Intensity Factor (SIF) ........................................... 20 2.2.2 A Model for Crack Length.............................................. 22 iv v 2.3 Continuum Damage Mechanics Fundamentals ............................... 23 2.3.1 Damage in Mechanical Properties...................................... 23 2.3.2 Damage in Flow Properties............................................. 25 2.3.3 Discretization of the Coupled Problem ................................ 27 2.3.3.1 Flow Problem ................................................. 27 2.3.3.2 Mechanics Problem ........................................... 28 CHAPTER 3 SEQUENTIALLY COUPLED SIMULATOR ........................................... 32 3.1 Coupling Methods ............................................................... 32 3.1.1 Fully Coupled or Simultaneous Solution Method ..................... 32 3.1.2 Iteratively Coupled or Sequential Solution Methods ................. 33 3.1.3 One-way Coupled Sequential Solution Methods ...................... 34 3.2 Proposed Methods ............................................................... 36 3.2.1 Extended Fixed Stress Method with Explicit Damage (CFM) ...... 36 3.2.2 Extended Fixed Strain Method ........................................ 37 3.2.3 Extended Fixed Stress Method with Implicit Damage (CFMD) .... 37 3.2.4 One-way Coupled Method for Multiphase Geomechanics ............ 37 3.3 Implementation .................................................................. 40 3.3.1 Fracture Distribution in Space ........................................ 41 3.3.2 Fracture Evolution in time ............................................. 41 CHAPTER 4 SIMULATOR VALIDATION ............................................................ 42 4.1 Uniaxial Consolidation .......................................................... 42 4.1.1 Eect of boundary condition ........................................... 42 4.2 Model Sensitivity ................................................................ 45 4.2.1 Time Step Renement Eect ........................................... 46 4.2.2 Eect of Mesh Renement.............................................. 46 CHAPTER 5 REPRESENTATIVE NUMERICAL SIMULATIONS.................................. 51 5.1 Two-way Coupled Single-phase Flow-Geomechanics Simulator .............. 51 5.1.1 Physical Model .......................................................... 52 5.1.2 Coupling and Sequential Iteration In uences ......................... 53 5.1.2.1 Coupling Methods ............................................ 53 5.1.2.2 Adaptive Coupling ............................................ 54 5.1.2.3 Damage in Flow and Mechanical Properties ............... 55 5.1.3 Eect of Biot Coecient................................................ 56 5.1.4 Eect of Young's Modulus ............................................. 64 5.1.5 Case I: Injection Well at the Center of the Domain (Base Case) .... 67 5.1.6 Case II: Eect of Variable Injection Rate ............................. 73 5.1.7 Case III: Injecting at Constant rate for Dierent Times ............. 76 5.1.8 Eect of Fracture Orientation .......................................... 79 5.1.9 Multiple Fracture Families.............................................. 82 5.1.10 Boundary Stress Anisotropy ............................................ 84 5.2 One-way Coupled Simulator for Multiphase Flow and Geomechanics....... 85 5.2.1 Comparison with Single-phase Simulator.............................. 85 vi 5.2.2 Physical Model of an Injector-Producer Pair.......................... 86 5.2.3 Eect of Damage in Flow properties .................................. 87 5.2.3.1 Relative Permeability Eect ................................. 93 5.2.3.2 Fluid Viscosity Eect......................................... 95 5.2.4 Eect of Mechanical Properties ........................................101 5.2.4.1 Young's Modulus Eect ......................................107 5.2.4.2 Biot Coecient Eect ........................................109 5.2.5 Eect of Damage in Flow and Mechanics Properties .................115 5.2.6 Eect of Well Placement................................................121 5.2.7 Eect of Multiple Fracture Families ...................................127 CHAPTER 6 CONCLUSION AND FUTURE WORK ...............................................144 6.1 Conclusion .......................................................................144 6.2 Recommendation ...............................................................146 References..................................................................................149 List of Figures 1.1 The three modes of crack propagation: (I) opening mode, (II) sliding mode, and (III) tearing mode [43].............................................. 7 1.2 Application of the cohesive zone element method in modeling fracture propagation. Cohesive elements are inserted along the boundaries of bulk elements lying on the two sides of the potential fracture [53]. ............... 9 1.3 A crack cutting the simulation mesh [12] ..................................... 10 1.4 Discrete fracture networks in two and three dimensions [3]. Fracture- fracture intersections are zero-dimensional nodes in 2D and one-dimensional lines in 3D. ....................................................................... 11 1.5 Lorenz Plot [68]. ................................................................. 12 2.1 Left: 3D domain of a naturally fractured reservoir containing microc- racks and larger fractures. Right: A two-dimensional conceptual model of a horizontal cross-section of the fractured reservoir showing multiple pre-existing microcracks. The model domain containing the randomly distributed, non-interacting microcracks is considered to be undamaged before the wells become active. ................................................. 14 2.2 Stress eld in a REV around a crack of half-length a k induced by injec- tion in a horizontal well of half-length L w in a domain under tectonic compressions h and H . is the fracture direction measured from the x-axis. t is the tension at the crack tip due to the normal projections of the eective stress and the deviatoric stress tensors. p is pore pres- sure induced by injection. The distance between the well and the crack is characterized with lengths r, r 1 , and r 2 , and angles , 1 , and 2 . The REV size is exaggerated with respect to the well size. ...................... 17 2.3 Bi-linear nite element shape functions for the four nodes of a quadrilat- eral element [49] ................................................................. 29 3.1 Iteratively coupled methods for coupled ow and mechanics. (A) Drained, (B) Undrained (C) Fixed-strain and (D) Fixed-stress methods. p = p j+1 p j , m = m j+1 m j , _ v = (@ v =@t) j+1 (@ v =@t) j , and _ v = (@ v =@t) j+1 (@ v =@t) j . ...................................................... 33 3.2 OC: The one-way coupled poromechanics-damage mechanics sequential solution algorithm. Here,C 1 = (A e =t) M 1 + T C dr 1 andC 2 = 0. See text for denitions of other quantities. Time level information (n or n + 1) is shown only for pressure, displacement, stress, and the fracture length. Damage Mechanics Method (DMM) is used to update the poromechanical properties. ................................................ 35 vii viii 3.3 CFM: The coupled ow and mechanics sequential solution algorithm. Here,C 1 = (A e =t) M 1 + T C dr 1 andC 2 = (A e =t)( T C dr 1 ). See text for denitions of other quantities. Time level information (n or n + 1) is shown only for pressure, displacement, stress, and the frac- ture length. Damage Mechanics Method (DMM) is used to update the poromechanical properties. .................................................... 38 3.4 CFMD: The coupled poromechanics-damage mechanics sequential solu- tion algorithm with implicit damage. ......................................... 39 4.1 Physical model of Terzaghi's uniaxial consolidation problem. ............... 43 4.2 Pressure, x-displacement, and y-displacement from our numerical simu- lation at the last time step. ..................................................... 43 4.3 Comparison of analytical and simulated pressure evolution at the top element. ......................................................................... 44 4.4 Comparison of analytical and simulated displacement prole across the domain at the last time step. ................................................... 44 4.5 Experimental Solution: 2D displacement xing all boundaries ............. 45 4.6 Experimental Solution: 2D displacement xed bottom boundary and rolling left and right boundaries. ............................................... 45 4.7 Time Step Renement Eect: Evolution of the normalized fracture length (a k ) d at a cell with dimensionless time (t d ) for dierent time step sizes. As the time step is decreased, the nonlinear behavior is captured better, and the fracture length solution slowly converges towards the true, unknown solution. .......................................................................... 47 4.8 Time Step Renement Eect: Evolution of the normalized permeability log(k x ) d at a cell with dimensionless time (t d ) for dierent time step sizes. 48 4.9 Time Step Renement Eect: Evolution of the normalized pressure (p d ) at a cell with dimensionless time (t d ) for dierent time step sizes. ......... 48 4.10 Mesh Renement Eect: Time evolution of the normalized fracture length (a k ) d for dierent mesh sizes. .................................................. 49 4.11 Mesh Renement Eect: The normalized permeability log(k x ) d evolution comparison for dierent mesh sizes. ........................................... 49 4.12 Mesh Renement Eect: The normalized pressure enhancement (p d ) for dierent mesh sizes. ............................................................ 50 5.1 Coupling Methods: Locations of the two fractures selected for detailed investigation. The two fractures are located in the left and right of the injector as shown above. ....................................................... 53 5.2 Coupling Methods: Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. All three coupling techniques are shown dierent behaviors .................................................. 54 5.3 One-way Coupling (OC) and two-way Coupled Flow Mechanics Damage (CFMD) methods. CFMD uses the Converged decision box to iterate over the ow and mechanics solutions within a sequential loop. OC does not use that. ..................................................................... 55 5.4 Adaptive Coupling: The four examined fractures locations for the single- phase simulator. ................................................................. 56 5.5 Adaptive Coupling: The normalized fracture propagation of the three methods for location(1) and location(2). ...................................... 57 ix 5.6 Adaptive Coupling: The normalized fracture propagation of the three methods for location(3) and location(4). ...................................... 57 5.7 Adaptive Coupling: (a) Normalized pressure p d , (b) Normalized frac- ture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . These plots show the comparison between the three methods at location (4). .................................................................. 58 5.8 Single-phase Flow. Eect of Damage at Location(1): Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the loca- tion. ............................................................................. 58 5.9 Single-phase Flow. Eect of Damage at Location(1): (a) Normalized sti- ness in x and y directions, C dr;ij =E i , (b) normalized Biot coecient in x and y directions, d , (c) normalized Biot modulus,M and (d) normalized permeability, log(kx;d). ......................................................... 59 5.10 Single-phase Flow. Eect of Damage. Snapshots of the permeability (k x ) dierence eld over 8 time steps, log(k x;d ) = log(k (x;d;F +M) ) log(k (x;d;F ) ). F +M indicates coupled ow and mechanics. Higher per- meability dierence (log(k x;d )) in north and south area of the domain. .................................................................................... 59 5.11 Single-phase Eect of Damage: Snapshots of the pressure dierence eld over 8 time steps, p d =p (d;F +M) p (d;F ) . Pressure is higher withF +M in north and south of the domain compared to the F case. ................. 60 5.12 Single-phase Eect of Damage: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N =S (N;F +M) S (N;F ) . More stress acts in north and south area. .................................................. 60 5.13 Single-phase Eect of Damage: Comparison of Lorenz Coecient (LC) between damage in ow and damage in ow and mechanics properties for permeability (k x ) in x-direction. .............................................. 61 5.14 Single-phase Eect of Damage: DamageD yy evolution iny-direction. The damage D xx equals to zero. .................................................... 61 5.15 Single-phase Eect of Damage: Snapshots of the stress ( xx ) for dam- age in ow case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. ........................ 62 5.16 Single-phase Eect of Damage: Snapshots of the stress ( yy ) for dam- age in ow case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. ........................ 62 5.17 Single-phase Eect of Damage: Snapshots of the stress ( xx ) for dam- age in ow and mechanics case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. ...... 63 5.18 Single-phase Eect of Damage: Snapshots of the stress ( xx ) for damage in ow and mechanics eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. ................... 63 5.19 Single-phase - Biot Coecient Eect: Pressure (p d ) and Fracture length (a k;d ) enhancement processes for dierent Biot coecients ( = 0:1 and = 0:8). ......................................................................... 64 5.20 Single-phase - Biot Coecient Eect: Snapshots of the permeability (k x ) dierence eld over 8 time steps, log(k x;d ) = log(k (x;d;0:8 )log(k (x;d;0:1) ).Higher permeability dierence (log(k x;d )) in north and south area of the domain for = 0:8 ....................................................................... 65 x 5.21 Single-phase - Biot Coecient Eect: p d =p (d;0:1) p (d;0:8) . Snapshots of the pressure dierence eld over 8 time steps. Pressure diuses more of damage ow and mechanics case in north and south of the domain. More pressure accumulates in that case. ............................................ 65 5.22 Single-phase - Biot Coecient Eect: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N = S (N;0:8) S (N;0:1) . More stress acts in north and south area. ........................................... 66 5.23 Single-phase - Biot Coecient Eect: Biot Coecient Eect: Comparison of Lorenz Coecient (LC) between damage in ow and damage in ow and mechanics properties for permeability (k x ) in x-direction. ............. 66 5.24 Single-phase - Young's Modulus Eect: Pressure (p d ) and Fracture length (a k;d ) enhancement processes for two Young's Modulus (E = 68 GPa and E = 136 GPa) . ................................................................. 67 5.25 Single-phase - Young's Modulus Eect: Snapshots of the dierence in k x permeability eld over 8 time steps. The permeability dierence is dened as log(k x;d ) = log(k (x;d;E=136) ) log(k (x;d;E=68) ). Higher permeability dierence in north and south area of the domain for E = 136 GPa. ....... 68 5.26 Single-phase - Young's Modulus Eect: Snapshots of the pressure dier- ence eld over 8 time steps, p d =p (d;E=68) p (d;E=136) . Pressure diuses more uniformly and higher for the case of E = 68 GPa. .................... 68 5.27 Single-phase - Young's Modulus Eect: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N = S (N;2E) S (N;E) . More stress acts in north and south area. ........................................... 69 5.28 Single-phase - Young's Modulus Eect: Comparison of Lorenz Coecient (LC) between Youngs Modulus (E = 68 GPa and E = 136 GPa) for permeability (k x ) in x-direction. .............................................. 69 5.29 Case I: Distribution of the initial fracture direction (t = 0) in the do- main. All the fractures are oriented along the general northwest-southeast direction. All the fractures belong to one single family i.e. N = 1, and each mesh element contains one fracture i.e. m k = 1=A e . ................... 70 5.30 Case I: Snapshots of the pressure eld over approximately three months with a time step of t = 11 day between the snapshots. Pressure diuses uniformly for the rst three time steps such that the pressure contours are elliptic in shape with the well lying at the center along the north-south long axis of the ellipse. Later the contours start to orient along the NW- SE direction following the permeability enhancement in that direction because of injection-induced fracture growth. ................................ 70 5.31 Case I: stresses (MPa) at two dierent time steps showing the creation of anisotropy and heterogeneity in the stress eld, which is homogeneous isotropic at t = 0. Anisotropy and shear, resulting from the eect of in Eq. (2.19), leads to an increase of the deviatoric stress S and its crack- normal projectionS N , which aects crack propagation through Eq. (2.34). The heterogeneity in the stress eld results from the heterogeneous distri- bution of , which makesS N and subsequentlyk andC dr heterogeneous. .................................................................................... 71 xi 5.32 Case-I: the change in the permeability eld k max at the same 8 time steps as in Figure 5. The change in permeability increases in the beginning (changing color from blue to red along the NW-SE corridor) and decreases later (red to yellow) as the induced stresses relax because of pressure diusion and fracture growth processes. ...................................... 71 5.33 Case-I: Locations of the two fractures selected for detailed investigation. The two fractures are located in NW and NE corners of the domain as shown by the circles. The fracture lengths are dierent in space because the plots are at the end of Case I simulation. ................................. 72 5.34 Case I: (a) The deviatoric stress S N is positive (tensile) at Loc (2) and negative (compressive) at Loc (1) because of the dierence in signs of , 1 , and 2 in Eq. (2.31). At both locations,S N evolves towards zero as the stresses relax with fracture propagation. The volumetric stress v at both locations are negative (compressive) because of near-wellbore expansion. This is also evident from the negative sign of v in the rst equation of Eq. (2.19) for both positive and negative signs of . (b) Normalized stiness in the x direction, C dr;1111 =E 0 , decreases with time because of fracture growth, (c) 11 increases with time suggesting reduction in the drained bulk modulus, (d) normalized permeability k max =k 0 vs. normal- ized stiness C dr;1111 =E 0 curves show the two implications of fracturing: increase in permeability and decrease in stiness. ............................ 73 5.35 Case I and II: cumulative injection prole in the constant and variable rate injection cases. The variable rate well starts o at an injection rate (slope of the dotted curve above) lower than that of the constant rate well. Later it gradually increases and surpasses the constant rate. The total injected volume is the same in both cases. .............................. 74 5.36 Case II: evolution of pressure over three months with a time step of t = 11 day during variable rate injection. Compare with Case I above. ......... 74 5.37 Case II: Changes in the permeability eld, k max , at dierent times with a constant time step of t = 11 day. Compare with Case I above. ........ 75 5.38 Case II: fracture length (a k ) and permeability (k max ) enhancement pro- cesses are dependent on the location and the well injection rate. CR and VR denote constant and variable rate Cases 1 and 2, respectively. Loc (1) and Loc (2) denote the NW and NE corner locations selected earlier... 75 5.39 Case III: the cumulative injected water volume (bbl) from the single hor- izontal well. ...................................................................... 77 5.40 Case III: The normalized fracture propagated length (a n;d ) for location(3) and location(4) with time. ...................................................... 77 5.41 Case III: the pressure behaviors at the injector and the observation loca- tion. ............................................................................. 78 5.42 Case III: (a) Normalized pressurep d , (b) Normalized fracture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . Two selected locations were investigated. Location (3) is located to the left of the horizontal injector and Location (4) in the right. ....................... 78 xii 5.43 Single-phase - Eect of Fracture Orientation: the eect of the fracture direction on evolution in the deviatoric tension S N , the local stress ratio xx = yy , and the maximum permeability k max at Location (1) at early times. S N is non-monotonic with because of the sin 2 term in Eq. (2.30). Growth in k max is proportional to the growth in S N mag- nitude. k max vs. S N shows that as the magnitude of deviatoric stress increases, the permeability increases for all . ............................... 79 5.44 Single-phase - Eect of Fracture Orientation: the eect of the fracture direction on evolution in the deviatoric tension S N , the local stress ratio xx = yy , and the maximum permeability k max at Location (2) at early times. k max is proportional to S N magnitude. ......................... 80 5.45 Single-phase - Eect of Fracture Orientation ( = 0 ): Evolution of pres- sure over three months with a time step of t = 11 day. Compare with Case I above . ................................................................... 80 5.46 Single-phase - Eect of Fracture Orientation ( = 90 ): Evolution of pressure over three months with a time step of t = 11 day. Compare with Case I above . ............................................................. 81 5.47 Single-phase - Eect of Fracture Orientation ( = 135 ): Evolution of pressure over three months with a time step of t = 11 day. Compare with Case I above . ............................................................. 81 5.48 Single-phase - Eect of Fracture Orientation - Location (A): Pressure (p) and Fracture length (a k ) enhancement processes are dependent on the location and the well injection rate. ........................................... 82 5.49 Distribution of the fracture direction for the two fracture families present in the domain. The fracture orientations are randomly distributed. ....... 83 5.50 Multiple Fracture Families for single-phase: Snapshots of the change in permeability eld for a fracture family with fracture directions distributed over the entire range of 0-180 . ................................................ 83 5.51 Boundary Stress Anisotropy for single-phase: Evolution of the normalized permeability k max =k 0 and the deviatoric tension S N with the stress ratio (SR) at the two locations. The plot shows that as SR increases i.e. the boundary stress anisotropy decreases the permeability increases at Location (1) and decreases at Location (2). S N decreases with SR at both locations. .................................................................. 84 5.52 Compare one-way coupled simulators for single-phase and multiphase (MRST) codes at location one. ................................................ 85 5.53 One-way coupling multiphase simulator : Pressure enhancement pro- cesses for 8 time steps. ......................................................... 86 5.54 Left: Two horizontal wells injector and producer placed in a naturally fractured reservoir containing microcracks and larger fractures. Right: A two-dimensional conceptual model of the two horizontal wells cross- section of the fractured reservoir showing the injector and producer inter- secting multiple pre-existing microcracks. The model domain containing the randomly distributed, non-interacting microcracks is considered to be undamaged before the well becomes active. ................................... 87 5.55 Multiphase - Damage in Flow properties Eect: (a) Pressure (p d ) and Fracture length (a k;d ) evolution are dependent on the location. ............ 88 xiii 5.56 Multiphase - Damage in Flow properties Eect: the upper two gures illustrate water saturation (%) and permeability change (logk x ) d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. 89 5.57 Multiphase - Damage in Flow properties Eect: The normalized pressure (p d ) enhancement at nine time steps. ......................................... 89 5.58 Multiphase - Damage in Flow properties Eect: The normalized perme- ability (logk x ) d ) growth in x-direction when fractures are perpendicular to horizontal wells . There is no fracture propagation in y-direction. ...... 90 5.59 Multiphase - Damage in Flow properties Eect: The normalized stress distribution ( xx ) d in x-direction at nine time steps. ........................ 90 5.60 Multiphase - Damage in Flow properties Eect: The normalized stress distribution ( yy ) d in y-direction at nine time steps. ........................ 91 5.61 Multiphase - Damage in Flow properties Eect: The normalized devia- toric stress (S N ) d at nine time steps. ......................................... 91 5.62 Multiphase - Damage in Flow properties Eect: the oil saturation (S o ) at nine time steps. .............................................................. 92 5.63 Multiphase - Damage in Flow properties Eect: the fractures propagation distribution for the last time step (t d = 100). ................................ 92 5.64 Multiphase - Damage in Flow properties Eect: Comparison of Lorenz Coecient (LC) between Model(1) and Model(2)for permeability (k x ) in x-direction. ...................................................................... 93 5.65 Multiphase - Relative Permeability Eect: Relative permeability curve for the second case. The residual oil saturation is 20% and the initial water saturation is 20%. ....................................................... 94 5.66 Multiphase - Relative Permeability Eect: watercut at the producer for both cases. ...................................................................... 94 5.67 Multiphase - Relative Permeability Eect for Location(1): The examined fractures location for multiphase ow-geomechanics simulator. .............. 95 5.68 Multiphase - Relative Permeability Eect: (a) Normalized pressurep d , (b) Normalized fracture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . The plots illustrate the comparison between the two cases for the middle fracture. This fracture is located between the injector and producer. ...................................................... 96 5.69 Multiphase - Relative Permeability Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ...... 96 5.70 Multiphase - Relative Permeability Eect: The dierence in normalized pressure at nine time steps. p d =p d;Relperm1 p d;Relperm2 . ............... 97 5.71 Multiphase - Relative Permeability Eect: The normalized permeability dierence in x-direction when fractures are perpendicular to horizontal wells. (log(k x )) d = log(k (x) ) d;Relperm1 log(k (x) ) d;Relperm2 . There is no fracture propagation in y-direction. ........................................... 97 5.72 Multiphase - Relative Permeability Eect: The normalized stress dis- tribution ( xx ) d dierence in x-direction at nine time steps, ( xx ) d = ( xx ) d;Relperm1 ( xx ) d;Relperm2 . .............................................. 98 xiv 5.73 Multiphase - Relative Permeability Eect: The normalized stress dis- tribution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;Relperm1 ( yy ) d;Relperm2 . ............................................... 98 5.74 Multiphase - Relative Permeability Eect: The normalized the deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;Relperm1 (S N ) d;Relperm2 . 99 5.75 Multiphase - Relative Permeability Eect for Relperm(2): the oil satura- tion (S o ) at nine time steps. ................................................... 99 5.76 Multiphase - Relative Permeability Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) Relperm(1) (S o ) Relperm(2) . ....................100 5.77 Multiphase - Relative Permeability Eect: the fractures propagation dis- tribution for the last time step (t d = 100). ...................................100 5.78 Multiphase - Relative Permeability Eect: Comparison of Lorenz Coe- cient (LC) between Relperm(2) and Relperm(2) for permeability (k x ) in x-direction. ......................................................................101 5.79 Multiphase - Fluid Viscosity Eect: (a) Normalized pressure (p d ) and fracture length (a k;d ) evolution depend on the location. ....................102 5.80 Multiphase - Fluid Viscosity Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ..................102 5.81 Multiphase - Fluid Viscosity Eect: The normalized pressure (p d ) dier- ence at nine time steps. p d =p d;=5 p d;=2 . .............................103 5.82 Multiphase - Fluid Viscosity Eect: The normalized permeability (logk x ) d dierence in x-direction when fractures are perpendicular to horizontal wells (log(k x )) d = log(k (x) ) d;=5 log(k (x) ) d;=2 . There is no fracture propagation in y-direction. .....................................................103 5.83 Multiphase - Fluid Viscosity Eect: The normalized stress distribution ( xx ) d dierence inx-direction at nine time steps, ( yy ) d = ( yy ) d;=5 ( yy ) d;=2 . .......................................................................104 5.84 Multiphase - Fluid Viscosity Eect: The normalized stress distribution ( yy ) d dierence iny-direction at nine time steps, ( yy ) d = ( yy ) d;=5 ( yy ) d;=2 . .......................................................................104 5.85 Multiphase - Fluid Viscosity Eect: The normalized deviatoric stress (S N ) d dierence at nine time steps, (S N ) d = (S N ) d;=5 (S N ) d;=2 . ....................................................................................105 5.86 Multiphase - Fluid Viscosity Eect for o = 2: the oil saturation (S o ) at nine time steps. .................................................................105 5.87 Multiphase - Fluid Viscosity Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) =2 (S o ) 0 =2 . .....................................106 5.88 Multiphase - Fluid Viscosity Eect: the fractures propagation distribution for the last time step (t d = 100). ..............................................106 5.89 Multiphase - Fluid Viscosity Eect: Comparison of Lorenz Coecient (LC) between o = 5 [cp] and o = 2 [cp] for permeability (k x ) in x- direction. ........................................................................107 5.90 Multiphase - Biot Modulus Eect: (a) pore pressure p, (b) permeability in x-direction k x , (c) water saturation S w and (d) watercut (WC) at producer. ........................................................................108 5.91 Multiphase - Biot's Modulus Eect: the pressure eld at 15 dierent time steps for M = 150 MPa. .......................................................108 xv 5.92 Multiphase - Biot Modulus Eect: the dierence in x-direction perme- ability when fractures are perpendicular to horizontal wells ( = 0). log(k x ) = log(k (x) ) M=150 log(k (x) ) M=100 . There is no fracture propa- gation in y-direction. ...........................................................109 5.93 Multiphase - Young's Modulus Eect: (a) Evolution in normalized pres- sure (p d ) and fracture length (a k;d ) depends on the location. ...............110 5.94 Multiphase - Young's Modulus Eect: the upper two gures illustrate water saturation (%) and permeability change at an element at the center of the domain between the injector and the producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ...................................................................110 5.95 Multiphase - Young's Modulus Eect: The dierence in normalized pres- sure at nine time steps. p d =p d;E1 p d;E2 . ................................111 5.96 Multiphase - Young's Modulus Eect: The normalized permeability (logk x ) d dierence in x-direction when fractures are perpendicular to horizontal wells (log(k x )) d = log(k x ) d;E1 log(k x ) d;E2 . There is no fracture prop- agation in y-direction. ..........................................................111 5.97 Multiphase - Young's Modulus Eect: Multiphase - Young's Modulus Eect: The normalized stress distribution ( xx ) d in x-direction at nine time steps, ( xx ) d = ( xx ) d;E1 ( xx ) d;E2 . .................................112 5.98 Multiphase - Young's Modulus Eect: The normalized stress distribution ( yy ) d in y-direction at nine time steps, ( yy ) d = ( yy ) d;E1 ( yy ) d;E2 . .112 5.99 Multiphase - Young's Modulus Eect: The normalized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;E1 (S N ) d;E2 . ..................113 5.100Multiphase - Young's Modulus Eect forE = 168 MPa: the oil saturation (S o ) at nine time steps. ........................................................113 5.101Multiphase - Young's Modulus Eect: the oil saturation dierence at nine time steps. S o = (S o ) E1 (S o ) E2 . ..........................................114 5.102Multiphase - Young's Modulus Eect: Comparison of Lorenz Coecient (LC) between the two cases of dierent Young's modulus, E1 and E2. The coecient is calculated using permeability in x-direction, k x . .........114 5.103Multiphase - Biot Coecient eect: (a) Normalized pressure (p d ) and fracture length (a k;d ) evolution depend on the location. ....................115 5.104Multiphase - Biot Coecient eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ..................116 5.105Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized pressure at nine time steps. p d =p d;=0 p d;=0:1 . .......116 5.106Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized permeability in x-direction when fractures are perpendicular to horizontal wells = 0. (log(k x )) d = log(k (x) ) d;=0:8 log(k (x) ) d;=0:1 . There is no fracture propagation in y-direction. .............................117 5.107Multiphase - Biot Coecient eect ( = 0:1 and = 0:8) : The normal- ized stress distribution ( xx ) d dierence inx-direction at nine time steps, ( xx ) d = ( xx ) d;o=0:8 ( xx ) d;o=0:1 . ......................................117 5.108Multiphase - Biot Coecient eect ( = 0:1 and = 0:8) : The normal- ized stress distribution ( yy ) d dierence iny-direction at nine time steps, ( yy ) d = ( yy ) d;o=0:8 ( yy ) d;o=0:1 . .......................................118 xvi 5.109Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized deviatoric stress at nine time steps, (S N ) d = (S N ) d;=0:8 (S N ) d;=0:1 . .....................................................................118 5.110Multiphase - Biot Coecient eect for = 0:1 : oil saturation (S o ) at nine time steps. .................................................................119 5.111Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): the oil saturation dierence at nine time steps, (S o ) = (S o ) o=0:8 (S o ) o=0:1 . .............119 5.112Multiphase - Biot Coecient eect for = 0:1 : fractures propagation distribution for the last time step (t d = 100). ................................120 5.113Multiphase - Biot Coecient eect: Comparison of Lorenz Coecient (LC) between o = 0:8 and o = 0:1 for permeability (k x ) in x-direction. .120 5.114Multiphase - Damage in Flow and Mechanics Properties Eect: (a) Pres- sure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. All three coupling techniques are shown dierent behav- iors. ..............................................................................121 5.115Multiphase - Damage in Flow and Mechanics Properties Eect: the up- per two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ...................................................................122 5.116Multiphase - Damage in Flow and Mechanics Properties Eect: The nor- malized pressure (p d ) dierence at nine time steps, p d =p d;F p d;F +M . ....................................................................................122 5.117Multiphase - Damage in Flow and Mechanics Properties Eect: The nor- malized permeability (logk x ) d growth in x-direction when fractures are perpendicular to horizontal wells (log(k x )) d = log(k (x) ) d;F log(k (x) ) d;F +M .There is no fracture propagation in y-direction. .....................................123 5.118Multiphase - Damage in Flow and Mechanics Properties Eect: The nor- malized stress distribution ( xx ) d dierence in x-direction at nine time steps, ( xx ) d = ( xx ) d;F ( xx ) d;F +M . .....................................123 5.119Multiphase - Damage in Flow and Mechanics Properties Eect: The nor- malized stress distribution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;F ( yy ) d;F +M . ......................................124 5.120Multiphase - Damage in Flow and Mechanics Properties Eect: The nor- malized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;F (S N ) d;F +M . .....................................................................124 5.121Multiphase - Damage in Flow and Mechanics Properties Eect: oil satu- ration (S o ) at nine time steps. .................................................125 5.122Multiphase - Damage in Flow and Mechanics Properties Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) F (S o ) F +M . .......125 5.123Multiphase - Damage in Flow and Mechanics Properties Eect: the frac- ture distribution at the last time step (t d = 100). ...........................126 5.124Multiphase - Damage in Flow and Mechanics Properties Eect: Compar- ison of Lorenz Coecient between two cases: Damage in ow property D(F) and Damage in ow and mechanics properties D(F+M). The coe- cient is calculated based on permeability in x-direction,k x , only. = 0 . ....................................................................................126 xvii 5.125Multiphase - Well Placement Eect : (a) Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. The three coupling methods show dierent behaviors. ...........................127 5.126Multiphase - Well Placement Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. ..................128 5.127Multiphase - Well Placement Eect ( = 0 and = 135): The dierence in normalized pressure dierence at nine time steps. p d = p d; =0 p d; =135 . .........................................................................128 5.128Multiphase - Well Placement Eect ( = 0 and = 90): The dierence in normalized pressure dierence at nine time steps. p d =p d; =0 p d; =90 . 129 5.129Multiphase - Well Placement Eect ( = 0 and = 135): The normalized permeability (logk y ) d in y-direction. ..........................................129 5.130Multiphase - Well Placement Eect for = 90: The normalized perme- ability (logk y ) d in y-direction when fractures are parallel to horizontal wells. There is no fracture propagation in x-direction. ......................130 5.131Multiphase - Well Placement Eect ( = 0 and = 135): The dier- ence in normalized stress in x-direction at nine time steps. ( xx ) d = ( xx ) d; =0 ( xx ) d; =135 . ......................................................130 5.132Multiphase - Well Placement Eect ( = 0 and = 90): The change in normalized stress inx-direction at nine time steps. ( xx ) d = ( xx ) d; =0 ( xx ) d; =90 . .....................................................................131 5.133Multiphase - Well Placement Eect ( = 0 and = 135): The dier- ence in normalized stress in y-direction at nine time steps. ( yy ) d = ( yy ) d; =0 ( yy ) d; =135 . ......................................................131 5.134Multiphase - Well Placement Eect ( = 0 and = 90): The dier- ence in normalized stress in y-direction at nine time steps. ( yy ) d = ( yy ) d; =0 ( yy ) d; =90 . .......................................................132 5.135Multiphase - Well Placement Eect ( = 0 and = 135): The dierence in normalized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d; =0 (S N ) d; =135 . .......................................................132 5.136Multiphase - Well Placement Eect ( = 0 and = 90): The dierence in normalized deviatoric stress ((S N ) d ) at nine time steps, (S N ) d = (S N ) d; =0 (S N ) d; =90 . ........................................................133 5.137Multiphase - Well Placement Eect for = 135: the oil saturation (S o ) at nine time steps. ..............................................................133 5.138Multiphase - Well Placement Eect for = 90: the oil saturation (S o ) at nine time steps. .................................................................134 5.139Multiphase - Well Placement Eect ( = 0 and = 135): the change in oil saturation at nine time steps, S o = (S o ) =0 (S o ) =135 . .............134 5.140Multiphase - Well Placement Eect ( = 0 and = 90): the change in oil saturation at nine time steps, S o = (S o ) =0 (S o ) =90 . ..............135 5.141Multiphase - Well Placement Eect for = 135: fracture distribution for the last time step (t d = 100). ..................................................135 5.142Multiphase - Well Placement Eect for = 90: fracture distribution for the last time step (t d = 100). ..................................................136 xviii 5.143Multiphase - Fracture Network Eect: Distribution of the fracture direc- tion for family (k=3) present in the domain. The fracture orientations are randomly distributed. ......................................................137 5.144Multiphase - Fracture Network Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) F 1 (S o ) F 2 . .........................................137 5.145Multiphase - Fracture Network Eect: The normalized pressure dier- ences (p d ) at nine time steps, p d =p d;F 1 p d;F 2 . ..........................138 5.146Multiphase - Fracture Network Eect: The normalized stress distribu- tion dierences in x-direction at nine time steps, ( xx ) d = ( xx ) d;F 1 ( xx ) d;F 2 . ........................................................................138 5.147Multiphase - Fracture Network Eect: The normalized stress distribu- tion dierences in y-direction at nine time steps, ( yy ) d = ( xx ) d;F 1 ( yy ) d;F 2 . ........................................................................139 5.148Multiphase - Fracture Network Eect: The normalized deviatoric stress (S N ) d dierence at nine time steps, (S N ) d = (S N ) d;F 1 (S N ) d;F 2 . ......139 5.149Multiphase - Fracture Network Eect (k=3): The normalized permeabil- ity (logk x ) d growth in x-direction. ............................................140 5.150Multiphase - Fracture Network Eect (k=1): The normalized permeabil- ity (logk y ) d growth in y-direction. ............................................140 5.151Multiphase - Fracture Network Eect (k=3): The normalized permeabil- ity (logk x ) d growth in x-direction. ............................................141 5.152Multiphase - Fracture Network Eect (k=3): The normalized permeabil- ity (logk y ) d growth in y-direction. ............................................141 5.153Multiphase - Fracture Network Eect (k=1): fractures propagation dis- tribution for the last time step. ...............................................142 5.154Multiphase - Fracture Network Eect (k=3): fractures propagation dis- tribution for the last time step. ...............................................142 5.155Multiphase - Fracture Network Eect (k=1 and k=3): Comparison of Lorenz Coecient (LC) between the two families for both permeabilities (k x and k y ). ....................................................................143 List of Tables 1.1 Dual-Porosity and Dual-Permeability theoretical models .................... 8 4.1 Simulation runs for the dierent time step range. ............................ 46 4.2 Simulation runs for dierent mesh renement cases. ........................ 47 5.1 Single-phase poromechanical properties. The subscript i for a property indicates the initial value of that time-dependent property................... 52 5.2 Multiphase poromechanical properties. ........................................ 87 xix Nomenclature a k fracture length in fracture family k m b critical fracture length m d k scalar damage variable for fracture family k - D damage tensor m 2 n k Unit normal vector of fractures in family k - m k number of cracks in family k per unit volume - porosity - k permeability tensor md c f uid compressibility psi 1 c m porous medium compressibility psi 1 K global stiness matrix Pa-m K s bulk modulus of the solid grain Pa M Biot modulus Pa N number of crack families - Biot coecient - density kg/m 3 viscosity cp p pressure psi strain - stress Pa t tip tensile stress in the crack-tip region Pa N total normal traction Pa S N deviatoric tension on the crack surface Pa C dr rank-4 drained elasticity tensor Pa fracture direction degree Poisson's ratio - E 0 Initial Young's modulus Pa ! Damage surface coecient - K I mode-I stress intensity factor (SIF) Pa p m K IC mode-I fracture toughness Pa p m xx Chapter 1 INTRODUCTION Signicant improvements in predicting complex reservoir performances by using ad- vanced methods and approaches have been adopted in the reservoir characterization techniques. However, the existence of natural fractures in low0permeability reservoir increases complication because they lead to complexity in geologic porous formations at multiple scales. Modeling fractures using varies methods is in uenced by fractures prop- erties, fracture distribution, and fracture propagation. With all of these aspects, fracture propagation will be indeed a target of interest in reservoir modeling and simulation. An applicable fracture modeling for any given fracture network can enhance the selection of well placement (orientation) and better estimation of permeability changes as well as scheduling production/injection rates, improving recovery factor. Therefore, essential studies for researchers is to apply various test models for the industry to evaluate the subsurface properties. An eective solution framework uses a sequential solution of three coupled problems: (1) mechanical deformation (2) uid ow, and (3) fracture propaga- tion. We use the nite element method to solve the quasi-static mechanical deformation problem, the nite volume method to solve the Darcy ow problem in the matrix, and a damage mechanics approach to model propagation of pre-existing spatially distributed fractures and the ensuing evolution in ow and mechanical properties of the reservoir. We use the xed stress scheme to couple the mechanical problem with the ow problem sequentially. The two problems are coupled with the fracture propagation problem via stresses, pressure, elastic moduli, and permeability. The literature review presented in chapter 1 provides a brief summary of the natural fracture. This chapter also discussed the existing approaches for geomechanical sim- ulation in fractured reservoirs. The damage mechanics methods (DMM) for fracture propagation is presented, including situations of the complex fracture network, where 1 2 discontinuous fracture representation can be complicated. The last part of the chapter will discuss the rock failure theory a fracture propagation mode. In chapter 2, we will discuss the mathematical modeling and governing equations used in our numerical simulation. Continuum mechanics ideas are used to represent the uid and solid skeleton as overlapping continua. The governing equations of uid ow and mechanics are obtained from mass and linear momentum balances. These equations are used in Biot's self-consistent theory of poroelasticity. The poroelasticity theory describes the relationship between the changes in total stress, the pore pressure, and the uid mass content in space and time. We use the continuum damage mechanics theory to model propagation of discrete fractures in the poroelastic domain. We will identify the mathematical terms and coecients in the governing equations that couple the three processes of uid ow, mechanical deformation, and fracture propagation. Chapter 3 discusses the sequential coupled simulator that we developed to solve the three sub-problems of ow (single-phase and multiphase), mechanical equilibrium, and frac- ture propagation sequentially. We discuss dierent types of sequential coupling methods and their advantages and disadvantages. We highlight that modeling and simulation of ow and deformation processes in fractured reservoirs using a computational mesh that explicitly represents individual fractures and conform to individual fracture geometry is computationally prohibitive. As a result, fractures are often treated as lower-dimensional objects within a higher dimensional matrix e.g. fractures as lines in a 2D domain. Even with such a simplied fracture representation, it is quite challenging to model the dy- namics of ow and fracture propagation in naturally fractured formations because of the complex geometry of the fracture network and the scale discrepancy between the fractures (low-dimensional objects) and the formation matrix they intersect. here are various types of coupling discussed in this chapter. Chapter 4 presents representative numerical simulations with inputs and outputs to verify the numerical formulation and its code implementation. Chapter 5 presents the implementation and validation of our numerical simulators where will be running dierent sensitivity analyses in which we will be discussing our results. The model approach is , a validated modeling, by using several synthetic, but realistic, test cases. We will illustrate the onset and evolution of permeability changes from uid injection and withdrawal. Using a two-way coupling method, we will use our single-phase ow simulators using dierent scenarios to estimate permeability changes. Similarly, however using one-way coupling, we will run the multiphase ow simulator. This will be examined through well rate, well type, fracture orientation, fracture distribution, fracture family and stress boundary stress. Also, we will be testing the uid ow eect in the multiphase ow simulator by comparing the water breakthrough and oil recovery. 3 1.1 Naturally Fractured Reservoirs The fractured low-permeability oil reservoir is one of the main resource bases for en- hanced reserves and productivity. A naturally fractured formation contains microcracks and fractures embedded in the host rock matrix across multiple length scales, which makes it dicult to obtain a petrophysical description of such reservoirs both in the eld and in the lab. Low permeability and low primary recovery from such reservoirs often favor application of horizontal drilling, hydraulic stimulation, and secondary re- covery technologies such as water ooding or gas injection. Petroleum reservoir engineers rely on numerical simulation tools to optimally assign well rates, accurately predict well performance, and eectively place wells in the reservoir subjected to such advanced technologies to oset the additional cost with additional recovery. However, fractured reservoirs are dicult to simulate numerically because of the nonlinearity of the govern- ing and constitutive equations and the presence of multiple length scales from a single fracture to the matrix to multiple interacting fractures. In order to analyze the fractures, we cannot ignore their geometries, spatial distribution, and propagation dynamics under production/injection-induced loading because they determine permeability and poroe- lastic properties of the reservoir. Therefore, the simulation model should be capable of modeling the changes in rock properties during production/injection if the goals is to accurately predict the oil recovery from such reservoirs. Another challenge is to nd the relation between operational decisions, e.g. well locations and ow rates, and subsurface stress. 1.1.1 Denition of a Fracture In general, a fracture is a surface of discontinuity in a brittle material. Joints, bedding planes, rock cleavage, foliation, shear zones, faults, and other contacts can be grouped under the denition of fracture [83]. We can dene any area with closely-spaced and highly interconnected discrete fractures as fracture zones. They could be very large zones and may contain both short (e.g. at centimeter scale) and long (e.g. at kilometer scale) fractures. For small discontinuities, fractures can be described as joints or breaks in rocks. A relatively long fracture with a signicant shear slip discontinuity is a fault. Dierent fractures have dierent curvatures, lengths and apertures. They propagate along dierent directions depending on the local rock composition and the state of stress around them. This results in dierent spatial distributions of the fracture number density, the fracture length, and the fracture orientation. Therefore, numerical models often use random space functions to randomly distribute multiple pre-existing fractures in the model domain to capture the ow and mechanical deformation response observed 4 in lab and eld data. A specic family or network of fractures, e.g. fractures trending northwest-southeast in a eld, can be distributed in a model with estimated values of the average fracture length and fracture spacing. However, under externally applied pressure or stress, many small fractures interact with each other and coalesce to create a larger fracture and/or a new fracture network. This can change the ow eld through the fractured region by altering the hydraulic conductivity. The three possible modes of fracturing are mode-I (tensile), mode-II (in-plane shear) and mode-III (anti-plane shear). In mode-I fractures, the walls of the fracture move orthogonally to the direction of the fracture plane. In mode-II and mode-III, the walls move in the plane of the fracture. During hydraulic fracturing, thermal contraction and diagenetic shrinkage, mode-I fracturing is dominant. In heterogeneous rocks or under complex stress states, a combination of these modes might occur together e.g. in the borehole breakout region. When these modes occur collectively, the fractures are called as mixed-mode fractures. 1.2 Theory of Poromechanics Poroelasticity is a continuum theory that is used to study and analyze mechanical re- sponse of uid-saturated rocks containing interconnected pores. The theory states that any porous material subjected to a change of stress, experiences matrix deformation and changes in pores size and shape. For a saturated rock, the presence of the uid alters the material stiness and causes ow between regions because of the dierence in uid pressure. The poromechanics theory was developed originally by Karl Terzaghi [87] to explain 1D consolidation of soil columns under quasi-static equilibrium. The work was extended to three dimensions by Maurice Biot [11]. Since then, multiple researchers have made signicant contributions to the theory, see e.g.[17]. The original theory of Terzaghi postulated that the deformation of soil results from the rolling and sliding of solid grains in the porous matrix instead of the expansion or contraction of the grains themselves. This rearrangement of the grains imparts a compressibility or bulk modulus to the porous skeleton that is dierent from the modulus of the grains or of the pore uid. For water-saturated soils considered by Terzaghi, the skeleton compressibility is larger than either the water compressibility or the solid grain compressibility. There- fore, he assumed incompressible grain and incompressible uid during his analysis. This assumption was relaxed by Biot, which resulted in the denition of the Biot coecient dened in terms of the skeleton modulus and the solid modulus. 5 1.2.1 Stress-Strain Relation and Poroelastic Properties of a Rock The relation between stress and strain in a poroelastic solid is usually dened by Hooke's Law, which is independent of time and loading history. It assumes a linear, proportional relation between eective stress and strain, i.e change in strain results from a change in eective stress. Using the general form of Hook's Law, any stress component near a point can be dened in terms of a linear function of strains around that point. For a 2D domain, this gives 0 ij =C ijkl kl (1.1) whereC ijkl is the elasticity tensor, 0 ij is the eective stress tensor, and kl is the strain tensor. For a homogeneous isotropic poroelastic body under a plane stress conguration, theC matrix is given by C ijkl = E 1 2 2 6 6 4 1 0 1 0 0 0 1 2 3 7 7 5 (1.2) where E is Young's modulus and is Poisson's ratio of the porous medium. Depend- ing on the ow boundary conditions (constant pressure drainage boundary or zero- ux undrained boundary), we can dened the drained Young's modulus E dr and Poisson's ratio dr or the undrained Young's modulusE u and Poisson's ratio u [94]. The drained bulk modulusK dr of a medium undergoing 3D deformation is related to E dr and dr as K dr =E dr =3(1 2 dr ). Consider a 2D porous domain with a porosity saturated with a uid and loaded with isotropic total stresses ( xx and yy ) under undrained conditions. This means no uid is allowed to leave the medium. The porosity represents the volume of the pores per the total volume of the material. The change of uid volume is aected by the increment of the pore pressure (p), the pressure in the pore uid is denoted by p and in the soil particles of magnitude [92], the change of the volume can be dened as: V f =C f pV (1.3) whereV f is the uid volume,V is the total volume of the element. The original volume of the pores is V . C f is the uid compressibility can be dened as: C f = 1 V f @V @p (1.4) 6 andK f = 1=C f is the bulk modulus of the uid. The volume change of the particles is, because the original volume of the particles is (1)V : V r = (1)C r pV (1.5) C r is the compressibility of the rock matrix material. The entire volume change of the porous medium includes the deformations of porous media due to the rearrangement of the particles. In the linearized system, it was assumed in the beginning that the change in volume can be dened as follows: V =C m ( p)V (1.6) where C m is the compressibility of the porous medium, sometimes called it bulk com- pressibility C b . The compression modulus (or bulk modulus) K b is the inverse of the compressibility C m of the porous medium (K b = 1=C m ). During the drained condition, whereC m = 1=K dr , the compressibility of the porous medium is obtained by measuring the volumetric strain because of changes in applied stress while holding pore pressure constant. The poroelastic expansion coecient describes how much the medium volume changes because of a pore pressure change while holding the applied stress constant. This denes the Biot modulus M as: 1 M =C f ( p)V (1.7) 1.3 Fracture Mechanics Fracture mechanics was developed initially to study the failure of materials (i.e. glass, metal, and engineering structures). Using similar fracture conditions, it was developed for rock fracture and other surface and subsurface failure such as in-situ stress identica- tion, earthquake movements, plate tectonics and uplift and erosion of crustal rocks [96]. The Mohr-coulomb or Hoek-Brown strength criterion describes phenomenological failure criteria and used as appropriate to determine rock failure strength and onset time. Be- cause of the complexity of the rock failure process, brittle and semi-brittle rocks mostly rely on the process of fracture initiation and propagation [30]. Most of the classical theories of continuum mechanics were developed using the experimental studies of the behavior of the metal in which dierent behaviors were observed such as linear and non-linear elastic, plastic and brittle failure. We used these constitutive models for subsurface rocks with the understating of the medium heterogeneity by considering the analogous behavior. 7 1.3.1 Grith Cracks Cracks propagate because of the applied load. Grith [35] used the thermodynamic theory to identify the criterion for fracture growth. His work concentrates on the energy required to create new surfaces of a certain area. The crack propagation criterion is dened such that crack propagates when the released elastic strain energy equals the required energy to create a new crack surface. The required tensile stress to create the new crack surface is given as follows: 0 = 2E s a 1=2 (1.8) where 2a is the fracture length, and s is the surface energy per unit area. Three dierent stress intensity factors can be dened for the three dierent modes of fracturing Mode I or opening mode, Mode II or sliding mode (In-plane shear), and Mode III tearing mode (Figure 1.1). The crack propagation condition can be dened in terms of the fractures toughness, K IC , which is the critical value of the stress intensity factor for the crack to grow. In 2D, there are two modes: (1) Mode I, which is in response to a normal load Figure 1.1: The three modes of crack propagation: (I) opening mode, (II) sliding mode, and (III) tearing mode [43] applied to the crack, and (2) Mode II, which corresponds to in-plane shear. 1.4 Fracture Modeling Modeling all types of fractures including connected and not connected pathways of nat- ural fractures is very complex. In the subsurface, these natural fractures create complex 8 networks and control the geomechanical and hydrological behaviors beside to the change of formation properties. For numerical modeling, it is necessary to identify the eect of natural fractures because of the heterogeneous stress elds and the complexity of the uid ow pathways. Our study used the damage mechanics method. It is a part of contin- uum mechanics that includes micro-scale eects into the continuum scale model through the damage variable [20]. Lemaitre [55] used the damage mechanics approach to model the deformation of metal under dierent loads. Then, Mazars [46] and Pijauder-Cabot continued the work to apply the damage mechanics techniques in a porous material application. Part of damage implementation is the simulation of hydraulic fracturing- induced permeability enhancement in shale reservoir addressed by [100]. In this section, dierent fracture models were discussed and highlight the advantages and disadvantages of each approach. 1.4.1 Dual Porosity and Dual Permeability Models Naturally fractured reservoirs are divided into two systems: matrix and fracture sys- tems. In the subsurface, dierent layers contain dierent formation distribution that contributes to uid transport. For the ow simulation modeling, formation properties (mainly porosity and permeability) are distributed in the model and each has a dierent in uence on the uid ow eld. The dual-porosity model [7] assumed that the matrix porosity has the largest volume of uid and the fractures transport uids to well-bore because they have extremely higher permeability. These fractures need to be identied and distribute their networks correctly in the reservoir simulator. Also, it is necessary to understand their behaviors because of reservoir changes. In general, there are four theoretical models listed in Table 1.1 for the naturally fractured reservoir as follow: Table 1.1: Dual-Porosity and Dual-Permeability theoretical models NFR Model Analytical Semi-analytical Dual porosity and Single per- meability De Swaan [85], Er- shaghi [42], Warren [95] Ershaghi [29], Bui [86] Dual porosity and Dual per- meability Hu [38] Van [4] Flow Models of Horizontal Well Production in a Single Medium Frick [89] Duan [97] Flow Models of Horizontal Well in NFR Du [24] Ozkan [67] 9 These models are very limited because of the explicit fracture assumptions. They also didn't account for the eect of mechanical properties changes. 1.4.2 Cohesive Zone Elements A very common fracture modeling approach in nite element analyses is the cohesive zone element approach to study fracture in a discontinuity system. The concept of a cohesive zone in fracture mechanics has been outlined earlier [7, 25] and are often referred to as Plastic Strip Yield Models. Figure 1.2 shows an application of the cohesive zone method. The cohesive zone method describes the creation of a fracture as the separation of the elements on the two sides of the fracture tip under the force of cohesive tractions. These region where this happens is called an extended crack tip or cohesive zone, which is occupied by cohesive elements. The fracture process of the cohesive zone approach Figure 1.2: Application of the cohesive zone element method in modeling fracture propagation. Cohesive elements are inserted along the boundaries of bulk elements lying on the two sides of the potential fracture [53]. works when a resisting force or traction acts on each crack surface as a function of the displacement jump across the surfaces. Hillerborg [36] applied a cohesive zone model in concrete to study crack growth in nite element analysis. Recently, the approach has been applied in hydraulic fracturing simulations. However, they need to initially dene the crack path in terms of locations and directions. The limitation of this application that it is mesh dependent fracture path solutions because the fractures can propagate only along element edges. 1.4.3 Extended Finite Element Method (XFEM) Cohesive zone approach was further extended to overcome the problem of mesh depen- dency by using the extended nite element method (XFEM). It allows for discontinuities to propagate directly through mesh elements along arbitrary directions as illustrated in Figure 1.3. This method was applied to model hydraulic fracturing of a poroelastic 10 medium [5] with one or a few individual fractures. XFEM can simulate crack propaga- tion using a direct integration technique by integrating XFEM method with a generalized algorithm of cutting triangulated domains [12]. Figure 1.3: A crack cutting the simulation mesh [12] The advantages of XFEM are as follows: It can model complex crack trends with multiple crack tips in the same element to model branching. It can model geometrically complex domains. It does not require re-meshing of the domain even when the discontinuities are large with respect to the domain size The applied algorithm uses a quadrature rule in order to create these branches where networking complexity is independent of the crack geometry. The micro-scale discontinuities exist and can aect permeability and other material behaviors and can be much smaller than the size of the nite elements. Although fracture is presented in both large and small scales in these models, the limi- tations are related to the implementation of the simplied assumptions [32]. 1.4.4 Discrete Fracture Network The discrete fracture network (DFN) uses an explicit representation for fracture geom- etry such as fracture direction, size, location, and shape. It also computes the relations 11 between each fracture and the fracture networks. These fracture networks as shown in Figure 1.4 are constructed using stochastic realization of the generated reservoir geol- ogy to map dierent types of rock fractures such as faults, joints, and bedding plans. Several advantages of DFN approach are maintaining natural fracture characteristics, representing a complex fracture network (including interaction, spacing, grouping, and intersection). The application can be suitable for 2D analysis more than a 3D prob- lem [75]. Figure 1.4: Discrete fracture networks in two and three dimensions [3]. Fracture- fracture intersections are zero-dimensional nodes in 2D and one-dimensional lines in 3D. 1.4.5 Lorenz Coecient A very useful representation of formation heterogeneity is the Lorenz Coecient tech- nique [68]. Heterogeneity in dierent properties such as the cumulative ow capacity, permeability and porosity can be computed using this technique by sorting the values from low to high magnitudes along a spatial direction. For example, the cumulative ow capacity over a series of stacked layers can be computed with the cumulative thickness in the vertical direction. The two values are plotted against each other as shown in Figure 1.5. The Lorenz Coecient L c is the measure of deviation of the curve from the diagonal straight line, which corresponds to a homogeneous distribution along that spatial direction. For heterogeneous formations such as a naturally fractured formation, the Lorenz Curve deviates from the diagonal. We can compute L c as twice the area between the two curves. To calculate th Lorenz coecient L c , we can use the equations below F c = P n elem i=1 k max;i max(k max ) (1.9) 12 Figure 1.5: Lorenz Plot [68]. S c = P n elem i=1 i max() (1.10) L c = 2 (L c ) 45 (L c ) property (1.11) where F c is the ow capacity and S c is the storage capacity. Chapter 2 MATHEMATICAL MODELING AND GOVERNING EQUATIONS In this section, we present the governing equations of coupled uid ow, mechanical deformation, and fracture propagation processes. The equations correspond to linear momentum balance of the porous medium, mass conservation of the uid components and the Grith failure condition for fracture growth. We use a classical continuum representation, where the uid and solid skeleton are viewed as overlapping continua. The ow-deformation equations are based on Biot's self-consistent theory of poroelastic- ity [10], and the continuum damage mechanics theory [55] is used to model the evolution in properties. The poroelasticity theory describes the relationship between the changes in total stress, the pore pressure, and the uid mass content in space and time. The continuum damage mechanics theory models the eect of discrete fracture propagation on hydraulic and elastic properties as a damage evolution process. We present our mathematical model for a two-dimensional domain (n dim = 2) representing a section of a reservoir containing micro-fractures. 2.1 Poroelastic Problem The model domain is assumed to be a 2D horizontal section of a reservoir containing large fractures and microcracks as shown in Figure 2.1. The section is a uid-saturated homogeneous rock embedded with randomly-placed idealized line microcracks. The domain is subjected to compressive loads representing the two horizontal stresses on its top and right boundaries, H and h , respectively. The stress ratio is dened as 13 14 SR = h = H . The bottom and left boundaries are xed in the normal directions. No- ow conditions are prescribed on all boundaries. The model is in hydrostatic and mechanical equilibrium, and the state of the rock before uid injection or production is regarded as the initial or undamaged state. Figure2.1: Left: 3D domain of a naturally fractured reservoir containing microcracks and larger fractures. Right: A two-dimensional conceptual model of a horizontal cross- section of the fractured reservoir showing multiple pre-existing microcracks. The model domain containing the randomly distributed, non-interacting microcracks is considered to be undamaged before the wells become active. 2.1.1 Single-phase Poroelasticity The following equation expresses the uid mass conservation during isothermal single- phase ow in a poroelastic medium [8]: dm dt +rw = f q w ; (2.1) where m = f (1 + v ) is the uid mass (per unit bulk volume), f is the uid density, w = f v is the uid mass ux,v is the Darcy velocity relative to the deforming skeleton, q w is the volume of water (subscript w) injected into the control volume per unit bulk volume, v is the volumetric strain dened as the trace of the strain tensor given below, and is the Eulerian or true porosity (ratio of pore volume to bulk volume at the current deformed conguration). evolves with pressure and stress. Initial porosity is 0 . Using Darcy's law in the absence of body forces (e.g. gravity), v = k rp; (2.2) where is the dynamic uid viscosity, k is the absolute permeability tensor andrp is the macroscopic pressure gradient. The quasi-static equilibrium of a porous medium 15 in absence of body forces (e.g. gravity) is governed by the linear momentum balance equation: r = 0; (2.3) where is the Cauchy total stress tensor containing normal total stresses and shear stresses. The eective stress describes the relationship between the total stress and the pore uid pressure, = 0 p1; (2.4) where 0 is the eective stress tensor, 0 1 is the Biot coecient, and 1 is the rank-2 identity tensor. We assume normal stresses are positive in tension and pressure is positive in compression. For time-dependent linear problems, where we are interested in solving for changes with respect to initial conditions, it is convenient to decompose quantities into their initial value, denoted with the subscript 0, and the change from the initial value. For example, the stress tensor can be decomposed as: = 0 +; (2.5) where 0 is the initial total stress and is the change in the total stress. The equilib- rium equation becomes: r 0 =r 0 +rp (2.6) The constitutive equations of isotropic poroelasticity in terms of increments in stress and pressure are [17]: = 0 p1 (2.7) m f;0 = v + 1 M p (2.8) whereM is the Biot modulus, 0 =C dr : is the increment in the eective stress tensor, C dr is the rank-4 drained isotropic elasticity tensor corresponding to the undamaged state, and is expressed in terms of Lame's parameter and the shear modulus G as C dr = ij kl +G( ik jl + jk il ) (2.9) where ij is the Kronecker delta matrix component. Young's modulus E and Poisson's ratio can be expressed in terms of and G. The linearized strain tensor is dened as := 1 2 (ru +r T u) (2.10) 16 in terms of the displacement vectoru. In Voigt notation, the stress and strain vectors under plane strain conditions can be written as: = [ xx ; yy ; xy ] T (2.11) = [ xx ; yy ; 2 xy ] T (2.12) which rendersC dr into a 3 3 rank-2 tensor. The 2D identity vector in Voigt notation becomes 1 = [1; 1; 0] T . The volumetric strain is v = tr(). The two coupling coecients of isotropic poroelasticity from Biot's theory are the Biot modulus M and the Biot coecient . They are scalars for an isotropic medium and dened as [17]: 1 M = 0 c f + 0 K s ; = 1 K dr K s (2.13) where c f = (1= f )d f =dp is the uid compressibility, K s is the bulk modulus of the solid grain, and K dr is the drained bulk modulus of the porous medium (independent of the pore uid). In 2D, K dr = 1 4 1 T C dr 1, which for plane strain yields K dr = +G. Therefore, v = 0 v p (2.14) 0 v =K dr v (2.15) where we assume linear isotropic elasticity for the mechanical response. To set the stage for volumetric total stress representation, we need to linearize the relation in Eq. (2.14) between volumetric total stress and volumetric strain v v;0 =K dr v (pp 0 ) (2.16) Thus, M is equal to [15]: M = K s 1 K dr Ks (1c f K s ) (2.17) The uid mass balance equation Eq. (2.1) needs to incorporate the eect of the me- chanical deformation. To apply the solid-to- uid coupling in the numerical solution, it is necessary to recognize the two contributions to the increment in uid mass content as shown in Eq. (2.8)): expansion of the pore space and increase in the uid pressure. Consider the model problem shown in Figure 2.2. We classify the cracks into crack families dened using the surface normal vector of the cracks. The subscript k is the crack family index (not permeability) and N is the number of crack families. A repre- sentative elementary volume (REV) containing a single elliptical crack with unit normal n k and subjected to total normal stresses xx and yy is shown. A local coordinate 17 Figure2.2: Stress eld in a REV around a crack of half-lengtha k induced by injection in a horizontal well of half-length L w in a domain under tectonic compressions h and H . is the fracture direction measured from the x-axis. t is the tension at the crack tip due to the normal projections of the eective stress and the deviatoric stress tensors. p is pore pressure induced by injection. The distance between the well and the crack is characterized with lengths r, r 1 , and r 2 , and angles , 1 , and 2 . The REV size is exaggerated with respect to the well size. system is dened with n k and the unit tangent vector along the crack long axis. The model boundary conditions are as follows. Zero normal displacements are prescribed on the left and bottom boundaries. Constant normal compressions are prescribed on the right and top boundaries equal to the minimum horizontal principal stress h and the maximum horizontal principal stress H . Zero uid ux is prescribed on all the bound- aries. Substituting Eq. (2.8) into Eq. (2.1) and assuming a slightly compressible uid, the uid mass balance equation in a homogeneous isotropic rock, in the pressure-stress formulation, becomes [45] 1 M + 2 K dr @p @t + K dr @ v @t rv =q w : (2.18) We can use the linear elastic fracture mechanics theory to approximately estimate the change in poroelastic stresses due to the presence of an injection well. Under the undrained deformation assumption (m = 0), valid for early times in low permeabil- ity rocks, the injection-induced poroelastic stresses xx , yy and xy at a point far 18 from the well can be given as [84]: xx + yy = 2p w 1 r p r 1 r 2 cos 1 + 2 2 = 2 v ; yy xx = 2p w rL 2 w (r 1 r 2 ) 3=2 sin sin 3 2 ( 1 + 2 ); xy =p w rL 2 w (r 1 r 2 ) 3=2 sin cos 3 2 ( 1 + 2 ): (2.19) wherep w is the injection-induced pressure, above the minimum principal stress, in the reservoir rock surrounding the wellbore. Eq. (2.18) can be re-written in terms of the volumetric strain as 1 M @p @t + @ v @t +rv =q w (2.20) Eq. (2.20) and Eq. (2.18) are two equivalent forms of the mass balance equation in single- phase poroelasticity. However, the two equations perform dierently when the ow and mechanics equations are solved sequentially and iteratively as part of a sequential solution scheme [44, 45, 51]. Eq. (2.18) is used in the xed stress sequential solution scheme where the@ v =@t term is frozen at its previous iteration value when this equation is solved forp. Eq. (2.20) is used in the xed strain sequential solution scheme where the @ v =@t term is frozen at its previous iteration value when this equation is solved for p. Even though the two schemes result in two-way coupling between ow and mechanics, the xed stress scheme is numerically more stable than the xed strain scheme and can handle strongly coupled problems in poroelasticity [51]. Here, we use the xed-stress scheme. 2.1.2 Multiphase Poroelasticity In multiphase ow system with n phase immiscible or partially miscible uid phases the equations of poroelasticity are dierent from those of single-phase and they show stronger nonlinearity due to saturation-dependent mobility and diusion behavior [18]. For an oil-water system, we denote the saturations by S w and S o and the phase pressures by p wat andp oil . From the denition of saturation, the summation of the phase saturations equals to 1: S w +S o = 1. The capillary pressure is dened as P c;ow (S w ) = p oil p wat . The mass conservation equation for multiphase system can be dened as: @ @t ((1 + v ) ` S ` ) +r ( ` v ` ) = ` q ` ; (2.21) where` is the phase index over all the phases (e.g. ` =w;o in an oil reservoir) andq ` is the volumetric source term of phase `. For multiphase ow, we can re-write the Darcy 19 equation, Eq. (2.2), as v ` = k r ` ` krp ` ; (2.22) where k r ` is the relative permeability of phase ` and ` is its dynamic viscosity. As per [17, 56{58], we know that coupled equations of multiphase ow and poroelasticity are nonlinear due to the nonlinearity of the capillary pressure relations e.g. P c vs. S w curve and relative permeability relations k r ` vs. S ` . Linearization of pressure and stress in the single-phase poroelasticity formulation has to be replaced with an incremental formulation that allows nonlinear evolution through time. As shown in Eq. (2.6), the eect of the total stress on the porous material composed of two parameters: (1) the eective stress which is responsible for deformation of the porous skeleton, and (2) the uid pressure which is responsible for the uid volume change: =C dr : n phase X =1 p 1 (2.23) where presents the Biot coecient for phase and the summation of all Biot co- ecients gives the Biot coecient of the saturated porous material, i.e. P = . Following [19, 56, 57], we assume =S . An equivalent pressure p E can be dened as the saturation-weighted average of phase pressures [57] such that the multiphase eective stress is = 0 p E 1; p E = X S p (2.24) 0 =C dr : (2.25) Some studies have proposed a dierent denition of p E to account for the energy at the interfaces between phases by including capillary pressures, P c;` for phases ` and , in calculation of p E [18, 45, 52, 66]. This denition imparts numerical stability to multiphase poroelasticity problems characterized by high capillarity ( 2 P c =K dr 1; P c > 0), e.g. deformation of soil columns and unconsolidated rocks. Here, we use the original denition because our focus is on consolidated oil reservoirs with negligible capillarity. The multiphase ow equation for phase ` in a deformable porous medium in the stress formulation is given as @ @t ` n phase X =1 M 1 ` + ` K dr 1 p + @ @t ` ` K dr 1 v +r ( ` v ` ) = ` q ` (2.26) where M 1 ` is one of the components of n phase n phase inverse Biot modulus tensor M 1 , which depends on phase saturationsS besides being dependent onK dr ,K s , and c ` . From the denition of c ` , we have ` = `;0 e c ` (p ` p `;0 ) . The formation volume factor 20 of a uid phase is introduced as B ` = `;0 = ` [6] such that B ` (p ` ) =B `;0 e [c ` (p ` p `;0 )] (2.27) 2.2 Fracture Dynamics and Propagation We use the classical Irwin-Grith theory [41] to derive the crack propagation criterion in terms of the stress intensity factor. Also, the continuum damage mechanics theory was used to update the petrophysical properties{permeability and stiness tensors{of a continuum model as functions of the damage tensor, which is dened using a damage variable. 2.2.1 Stress Intensity Factor (SIF) According to the linear elastic fracture mechanics theory [41, 54], the stress intensity factor around a crack of half-lengtha k and unit normaln k in an innite body subjected to the total stress and pore pressure p is K I (;p;a k ;n k ) = t p a k = r G c E 1 2 (2.28) where K I is mode-I stress intensity factor, t is the tip tensile stress in the crack-tip region, G c is the mechanical energy release rate (energy per unit crack surface area), and a k is the average fracture length in family k with unit normal n k . In the case of a nite size body containing multiple cracks, such as the physical model in Figure 2.1, the above equation is valid only under the assumption of discrete and non-interacting cracks. Here, we make that assumption because the models that account for crack in- teraction are mathematically more complex [26] and may not scale up for eld-scale simulation. An accurate determination of t is important because it controls the initia- tion and growth of new cracks, crack coalescence, and localization into macro-fractures, which lead to failure, strength reduction, seismicity, and permeability enhancement at macroscopic scales. The lack of knowledge of t around cracks and fractures is one of the major hurdles in understanding and predicting brittle processes in rocks. Based on some experimental observations, it has been proposed [81, 82] that the local stress near the crack tip can be expressed in terms of the eective normal traction and the normal component of the deviatoric stress on the crack surface as follows t = 0 N +f(a k )S N (2.29) 21 where 0 N = N +p is the eective normal traction, N =n k n k is the total normal traction, S N =n k Sn k is the deviatoric tension on the crack surface, and f(a k ) is a scalar function that accounts for the contribution of the deviatoric tension in inducing crack growth. Forn k = [ sin ; cos ] shown in Figure 2.2, S N = [( yy xx )=2] cos (2 ) xy sin (2 ); 0 N =S N + 0 v ; t = (1 +f)S N + 0 v : (2.30) Hence, t , and consequently K I , can be decomposed into an isotropic part and a devi- atoric part [23] similar to the decomposition of above. At early times, the injection- induced tractions can be evaluated analytically using Eq. (2.19). For example, S N = p w rL 2 w (r 1 r 2 ) 3=2 sin sin 3 1 + 3 2 4 2 ; (2.31) which shows us that S N can be dierent in the NE and NW quadrants of Figure 2.2 because of the dierent signs of . Within a quadrant, it will be dierent for cracks of dierent directions . For a crack not oriented along the principal stress directions e.g. 6= 0 ; 90 or for any crack under the injection-induced state of stress which is non-principal (Eq. (2.19)), the crack is under a mixed mode loading e.g. both mode-I and mode-II loading. Under these conditions, it is known that the crack grows following a curved path with the time-dependent crack direction (t) given by the maximum of G c ( ;a k ) corresponding to@G c =@ = 0, which changes witha k [27, 54]. In most models of brittle fracture the quasi-static evolution of a crack is governed by the energy release rate along the crack path [35, 65]. Here, we assume only mode-I failure to avoid the diculty associated with estimatingK II . Also, we do not model crack path to maintain the computational benets of our continuum damage mechanics approach. WhenK I starts to increase and reaches the equilibrium condition of mode-I crack prop- agation, K I =K IC , where K IC is the mode-I fracture toughness (a material property), the crack starts to grow under mode-I [16, 23, 41]. Propagation of a crack and associated evolution of the damage variable in a REV is governed by the propagation condition, which requires evaluation of the tip tensile stress t in the REV. t must be positive (tensile) in order for the fracture length to propagate, otherwise it will stabilize. This implies S N > 0 N =f for crack propagation. Once the crack length reaches a certain critical length b, crack coalescence begins and the material is considered to have failed because of self-sustaining mechanical instability. The function f(a k ) in Eq. (2.29) plays a role in crack coalescence because it can model local concentration of stresses around the crack as the crack length increases [80]. We also note thatf should decrease witha k 22 to model the drop in local tensile stress during crack propagation. f(a k ) has also been compared with the hardening-softening parameter in plastic deformation models [23, 82]. Here, we choose the following model for f [33, 100] f(a k ) = 8 < : !b=a k if (a k <b); ! if (a k b); (2.32) where ! is a material constant. A non-monotonic form of f(a k ) has also been pro- posed [81] to capture the stress relaxation at an early stage and the stress amplication for a k b because of the interaction between the growing crack and its neighbors. 2.2.2 A Model for Crack Length A fracture condition is required to determine the crack path and length accurately. A number of fracture conditions has been developed for a crack under critical mixed-mode loading [28, 39, 69, 70, 98]. Here, we consider the Grith condition. The propagation condition, which also represents a surface of the constant state of damage in the mean stress-deviatoric stress space [16], becomes p a k ( 0 N +f(a k )S N )K IC = 0 (2.33) which can be used to update the half-length a k . For early time growth when the cracks are below the critical length and do not interact i.e. a k <b, we obtain a k =!bR N + R I 2 R I q R 2 I 4!bR N (2.34) where R N =S N = 0 N and R I = (1= p )K IC = 0 N . For real positive values of a k , 0 N !b <S N < 1 4!b K 2 IC 0 N . Using Eq. (2.30), the lower bound condition can also be expressed as 0 v 1+!b <S N . In the subsurface, the cracks are generally under compression i.e. 0 N < 0 except during high pressure uid injection (e.g. hydraulic fracturing). A lithostatic state of stress (isotropic principal state of stress) such as the the one used below to initialize our Base Case yieldsSR = 1 andS N = 0, which implies R N = 0 and a k = [K IC =( p 0 N )] 2 as per Eq. (2.34). This can also be obtained from Eq. (2.28), Eq. (2.29) and the fracture equilibrium condition K I = K IC . Once a k is large enough for the cracks to start interacting among themselves through their local stress elds, unstable propagation and coalescence of cracks ensues culminating into a macroscopic fracture and structural failure of the material. We want to draw an analogy between the above model for computing the crack length a k during mode-I failure and a commonly used model for computing the fault slip s during mode-II (shear) failure 23 based on the Mohr-Coulomb condition. The Mohr-Coulomb failure condition is f , where =jn k N n k j is the shear traction magnitude on the fault/fracture plane, f = c f 0 N is the friction stress, c is cohesion, and f is the coecient of friction. After substitution, the shear failure equilibrium condition becomes + f 0 N c = 0; (2.35) Eq. (2.35) resembles Eq. (2.33) after the following substitutions in the former equation: with 0 N , f withf(a k ), 0 N withS N , and c withK IC =( p a k ). Therefore, as suggested earlier for plastic failure [23] and frictional failure [16], the Mohr-Coulomb failure surface is analogous to the mode-I failure surface given by Eq. (2.33). Dragon and Mroz [23] recognized the deviatoric stressS as a type of eective stress in crack development. As cohesion c or friction f changes, the Mohr-Coulomb failure surface moves in the stress space. Similarly, as a k changes in K IC =( p a k ) or in f(a k ), which are the respective counterparts of c and f , the mode-I failure surface moves. In fault modeling, f is usually modeled as a slip-dependent function e.g. the slip-weakening model [45], where f (s) = s [( s d )=s c ]s, s is the static friction coecient, d is the dynamic friction coecient, s c is the critical slip distance. This model of f (s) resembles the model of f(a k ) in Eq. (2.32), which further supports the analogy between these mode-I and mode- II failure models. A typical implementation of the Mohr-Coulomb failure model [1, 45] solves a contact problem to calculate the excess slip s required to accommodate the excess shear, = f . Since the excess slip must be associated with the displacement discontinuity across the crack surface, the contact problem is coupled to the elasticity problem in Eq. (2.3). The mode-I model proposed above for calculating a k (Eq. (2.34)) assumes that 0 N and S N do not change during crack propagation. This approximation is computationally much less expensive than solving the contact problem and is valid for microcracks that do not grow longer than the critical length and remain within the REV. 2.3 Continuum Damage Mechanics Fundamentals 2.3.1 Damage in Mechanical Properties We use the continuum damage mechanics theory [55] to dene the damage tensor in terms of the geometry and properties of propagating fractures. We assume an isotropic distribution of the microcracks in the initial undeformed conguration. In 2D, the 2 2 24 damage tensor is dened as [37, 80, 82] D := N X k=1 m k d k (n k n k ) (2.36) where m k = N k = denotes the number of cracks in family k per unit volume of the representative elementary volume (REV),d k is the dimensionless scalar damage variable. For the 2D case considered here, d k = (a 2 k a 2 0 )=a 2 0 (2.37) wherea 0 is the average initial length of the uniformly distributed fractures in the domain. Considering the initial state of the material containing a certain spatial density and dimensions (e.g. the half-lengths) of the cracks as the undamaged state, we can initialize d k = 0 at the beginning of the simulation. The damage variable evolves with the evolving state of stress through the crack length. Because of the growth of fractures and associated anisotropic damage, it is required to introduce the anisotropic elasticity tensor. Poroelastic properties can be updated as functions of the damage tensor as follows [60] C dr;ijkl (D) = ij lk +G( ik jl + jk il ) +A( ij D kl +D ij kl )+ (2.38) B( ik D jl + il D jk +D ik jl +D il jk ); 8i;j;k;l = 1; 2 (2.39) where the subscript k is not the fracture family index. Expanding the terms for plane strain yields C dr (D) = 2 6 6 4 + 2G + 2AD 11 + 4BD 11 +A(D 22 +D 11 ) AD 12 + 2BD 12 +A(D 22 +D 11 ) + 2G + 2AD 22 + 4BD 22 AD 12 + 2BD 12 AD 12 + 2BD 12 AD 12 + 2BD 12 G +B(D 22 +D 11 ) 3 7 7 5 (2.40) A and B are two constants dening the damage-induced modication of the strain en- ergy. The damage tensor componentD 11 indicates crack opening-induced damage along the x-direction on crack planes with unit normal in the x-direction. D 12 quanties sliding-induced damage along the y-direction on cracks with unit normal along x. This introduced anisotropy in the mechanical and ow properties of the damaged rock. Fol- lowing [15], we can write the components of the symmetric 2D Biot coecient tensor 25 (size 3 1 vector in Voigt notation in 2D) as functions of the damage tensor as follows 11 = 1 1 2K s 2( +G) + 2(A + 2B)D 11 +A(D 11 +D 22 ) 22 = 1 1 2K s 2( +G) + 2(A + 2B)D 22 +A(D 22 +D 11 ) 12 = 1 2K s 2(A + 2B)D 12 (2.41) Under anisotropic damage, the drained bulk modulus becomes K dr (D) = 1 4 (C dr;1111 +C dr;2222 +C dr;1122 +C dr;2211 ) = +G + (A +B)(D 11 +D 22 ) (2.42) where C dr;1111 , C dr;2222 and C dr;1122 =C dr;2211 are the components of plane strain elas- ticity tensor from Eq. (2.40). The Biot modulus, which is a scalar function, can be updated with the eects of damage by substituting K dr (D) in Eq. (2.17): M(D) = K s h 1 +G+(A+B)(D 11 +D 22 ) Ks i (1c f K s ) (2.43) We can re-write Equations (2.6){(2.8) as follows: r 0 =r 0 +r (p) =r 0 +rp +pr; (2.44) = 0 p; (2.45) m f;0 = : + 1 M p: (2.46) This has two implications compared to the isotropic equations. First, appearance ofr as a source term in Equation (2.44) suggests that spatial variation in the Biot coecient, which emerges because of spatially heterogeneous microcrack growth around a well, can contribute to deformation and enhance its heterogeneity. The simulation results discussed below validate this hypothesis of enhanced heterogeneity in the deformation eld. Second, appearance of 2 12 xy in the scalar product of the Biot coecient and the strain tensors in Equation (2.46) suggests a contribution to the change in uid mass because of shear strain. 2.3.2 Damage in Flow Properties Consider the REV dened earlier in terms of the porous matrix volume and the volume occupied by the microcracks, c . The Darcy ow velocityv is assumed to be given by the sum of a matrix ow component and a crack ow component, v =krp=, wherek =k 0 I +k c ,k 0 is the initial (undamaged state) isotropic permeability,k c is the 26 microcracking-induced permeability tensor (2 2 in 2D), and p is the global pressure gradient driving the mean ow through the domain. Here, it is assumed that the local pressure gradient driving the ow through a crack is the same as the global pressure gradient. This assumption is valid for randomly oriented cracks with relatively high permeabilities. k 0 does not change with damage and remains constant. We assume that k c (D) determines permeability changes caused by changes in the length, orientation, average aperture and number of microcracks. To determinek c (D), we express the crack ow part ofv as the average velocity through all cracks within [82] v = k 0 rp + 1 Z c v c d c ; (2.47) wherev c is the velocity in an individual crack, and integration over is equivalent to integration over c because cracks are only present in c part of . We assume the ow inside the crack to be laminar and parallel to the crack plane direction. For a given crack of orientationn k ,v c can be expressed as v c = e 2 k 12 (In k n k )rp; (2.48) where we used the assumption that the local pressure gradient is equal to the global pressure gradient, and there is no ow in the direction transverse to the crack plane. The positive scalar 0< 1 is introduced to account for the fact that not all parts of a crack may contribute to ow. The classic cubic law for ow rate through a crack created by two perfectly parallel surfaces, i.e. in-plane Poiseuille ow, is recovered with = 1. We assume a linear relation between the average aperture of microcracks and the local tensile stress such that e k ( t ;a k ) = (a k t )=E 0 ; where E 0 is Young's modulus of the undamaged material. Evaluation of the integral term in Eq. (2.47) requires estimation of the crack volume c , which depends on the number of cracks of all orientations within the REV that are contributing to ow. Let 0R k 1 be the crack connectivity coecient that represents the fraction of cracks contributing to ow. Above, we dened N k as the number of cracks oriented alongn k . Therefore, R k N k is the number of cracks oriented along n k and contributing to ow. The innitesimal crack volume for cracks oriented alongn k can be calculated as d c;k =R k N k e k a 2 k in 3D andd c;k =R k N k e k 2a k in 2D. This has to be averaged over all orientations 0 2 (all crack families) to obtain the average crack velocity. Therefore, we can rewrite the velocity in Eq. (2.47) as, v = k 0 rp + 2N k 1 2 Z 2 0 R k e k a k v c d (2.49) The crack permeability tensork c is regarded as a function of the crack orientationn k , average aperturee k , andN k . For ow through a single crack of uniform aperturee k and 27 unit normaln k , we obtain [82]: k c = 6 N k 1 2 Z R k a k e 3 k (In k n k )d : (2.50) Since the crack family index k varies with the crack direction , the integration over can be approximated with a summation over k, k c (D) = 6 N X k=1 m k R k a k e 3 k (In k n k ) (2.51) where m k = N k = . This model provides the components k xx , k yy , and k xy of the k tensor dened for the REV. The mass balance equation for coupled single-phase ow{poromechanics{damage me- chanics becomes @ @t M 1 (D) + T (D)C dr 1 (D)(D) p (2.52) + @ @t T (D)C dr 1 (D) +r k(D) rp =q w (2.53) We can also write this equation in the strain formulation as @ @t M 1 (D)p + @ @t ((D)) +r k(D) rp =q w (2.54) 2.3.3 Discretization of the Coupled Problem 2.3.3.1 Flow Problem We solve Eq. (2.52) using a nite volume method with a piecewise constant discretization of the pressure eld, two-point ux approximation of the Darcy ux term (second-order accurate) and the Backward Euler (rst-order accurate) time integration method [45]. The linear system of equations for the element-centered pressures can be written as follows A e t M 1 + T C dr 1 +T n+1 P n+1 + A e t T C dr 1 n+1 = A e t M 1 + T C dr 1 n P n + A e t T C dr 1 n +A e q n+1 w (2.55) where superscript n is the old time step at which we have the solution, and n + 1 is the new timestep at which a solution is sought. T is the transmissibility matrix composed of harmonic means of uid mobilities across element faces and, hence, depends on the 28 permeability tensor k(D). M is the matrix of element-centered Biot moduli in the mesh. The pressure vector, P , also includes unknown element pressures p w in the elements subjected to injectionq w . Area ofith element isA e;i = x i y i . For the strain formulation, the discretized system of equations is A e t M 1 +T n+1 P n+1 + A e t () n+1 = A e t M 1 P n + A e t () n +A e q n+1 w (2.56) We non-dimensionalize time in our simulations with the hydraulic diusion time scale: t ow c = x 2 S t k 0 ; t d = t t ow c (2.57) where x is the characteristic length given by the length of a typical mesh element, k 0 is the initial permeability value, and t d is the dimensionless time. S t = 1 M + 2 + 2G (2.58) S t is the storativity coecient dened in terms of initial or undamaged poroelastic property values. 2.3.3.2 Mechanics Problem We discretize the mechanics problem using the standard Galerkin nite element method with bilinear quadrilateral elements. Initial values of rock and uid properties such as Young's modulus, Poissons ratio, porosity, permeability, Biot coecient, and Biot modulus are assigned to each quadrilateral element and they are considered element- centered. U x and U y degrees of freedom are dened at the four nodes of each element. A local coordinate system in terms of is dened with origin at the center of the master element which is a square of size 2 in local coordinates. Linear shape functions are used in both directions of the element. The explicit representation of the four shape functions shown in Figure 2.3 are given below in terms of the local coordinates [49]: N s1 = 1 4 (1)(1) (2.59) N s2 = 1 4 (1 +)(1) (2.60) N s3 = 1 4 (1 +)(1 +) (2.61) N s4 = 1 4 (1)(1 +) (2.62) 29 Figure 2.3: Bi-linear nite element shape functions for the four nodes of a quadrilat- eral element [49] The strain-displacement matrix, B, of element i can be obtained from nodal strain- displacement matrices as follows [49]: B i = 1 det(J) [B 1 B 2 B 3 B 4 ] (2.63) whereB n matrix for noden is given in terms of partial derivatives of the shape functions, B n = 2 6 6 4 N sn;x 0 0 N sn;y N sn;y N sn;x 3 7 7 5 (2.64) where N sn;x =@N sn =@x. Therefore,B i becomes B i = 2 6 6 4 @Ns 1 @x 0 @Ns 2 @x 0 @Ns 3 @x 0 @Ns 4 @x 0 0 @Ns 1 @y 0 @Ns 2 @y 0 @Ns 3 @y 0 @Ns 4 @y @Ns 1 @x @Ns 1 @y @Ns 2 @x @Ns 2 @y @Ns 3 @x @Ns 3 @y @Ns 4 @x @Ns 4 @y 3 7 7 5 (2.65) 30 The Jacobian matrixJ for local-to-global coordinate transformation is dened as usual. Its determinant is det(J) = 1 8 h x 1 x 2 x 3 x 4 i 2 6 6 6 6 6 4 0 1 1 1 0 + 1 1 0 + 1 1 + 1 0 3 7 7 7 7 7 5 2 6 6 6 6 6 4 y 1 y 2 y 3 y 4 3 7 7 7 7 7 5 (2.66) The element stiness matrix for ith element in the local coordinate system is K i = Z 1 1 Z 1 1 B T C dr B dd (2.67) For plane stress, the elasticity tensor is C dr = E 1 2 2 6 6 4 1 v 0 1 0 0 0 1 2 3 7 7 5 (2.68) We follow the standard nite element recipe to obtain the weak form of the equilibrium equation: Z B T 0 d QP =F u (2.69) whereF u is the global nodal force vector that includes the eect of boundary loads and the coupling matrixQ is given by Q = Z B T 1N sp d (2.70) Since 0 =C dr : =C dr BU, we can write using Eq. (2.67): Z B T 0 d = Z B T C dr B d U =KU (2.71) We discretize the coupled equation using the standard nite element method with linear elements. The discretized linear system of equations that must be solved to calculate the nodal displacements is KU =F u +QP (2.72) whereK is the global stiness matrix containing Young's modulus and Poisson's ratio in the case of a linear elastic material considered here,U is the global nodal displacement vector, andP is the pore pressure vector. Given the time-dependency of the ow and fracture propagation problems, equation system (2.72) is solved at every time step. 31 The discretized system of equations is K n+1 U n+1 Q n+1 P n+1 =F u n+1 (2.73) where the load vectorF u n+1 at time leveln + 1 is known from the boundary conditions. The crack propagation process in uences the mechanics problem only through the change in the poromechanical properties, which is dierent compared to the approach based on discretized fractures e.g. cohesive zone, extended nite element, or interface element method. In the discretized fracture approach, a contact problem is solved on each crack surface at each timestep to obtain the traction and slip vectors on the surface. This increases the accuracy of the solution; however, at a signicant cost because the tractions are treated as unknowns similar to the nodal displacements. This approach is not suitable for modeling hundreds of fractures in a eld-scale simulator. It is also not required because discretizing individual fractures implies that the initial position, length, and direction of the fractures are known with certainty, which is rarely the case. Chapter 3 SEQUENTIALLY COUPLED SIMULATOR In this chapter, we discuss dierent solution methods to solve the coupled problem dened in the previous chapter. We describe in detail the sequential solution method that we used to implement our simulator. 3.1 Coupling Methods We present various coupling strategies to solve the coupled problem sequentially. The coupling between ow and mechanics problems is discussed rst because their solution provides us the primary unknowns, pressure and displacements. Next we discuss the damage mechanics problem to solve for fracture growth. In the last section of this chapter, we describe the coupling strategy we selected and its numerical implementation. 3.1.1 Fully Coupled or Simultaneous Solution Method Changes in stress and pore pressure in mechanically weak, stress-sensitive reservoirs re- quires two-way coupling between uid ow with mechanical deformation processes during simulation of subsurface injection and production. Fully coupled and iteratively coupled methods are the two main classes of the coupling methods [22, 79]. The fully-coupled method treats the medium as one system by solving ow and mechanics problems simul- taneously at each timestep [31, 47, 57, 71, 72, 93]. A fully coupled method is numerically unconditionally stable, however it is very expensive because of the time required to solve large ill-conditioned linear systems generated by the method. The formulation below 32 33 Figure3.1: Iteratively coupled methods for coupled ow and mechanics. (A) Drained, (B) Undrained (C) Fixed-strain and (D) Fixed-stress methods. p =p j+1 p j , m = m j+1 m j , _ v = (@ v =@t) j+1 (@ v =@t) j , and _ v = (@ v =@t) j+1 (@ v =@t) j . summarizes the solution method mathematically: " U n P n # $ ! " U n+1 P n+1 # ; where $ = ( r = 0 dm dt +rw = f q w (3.1) 3.1.2 Iteratively Coupled or Sequential Solution Methods To reduce the computational cost and to re-use existing simulator codes developed for uncoupled ow and mechanics problems separately in the computational mechanics com- munity, several sequential solution methods have been proposed recently [22, 44, 50, 51, 62, 64, 78, 88, 90, 91]. There are four sequential solution methods: drained, undrained, xed-strain, and xed-stress methods as shown in Figure 3.1. In the xed stress or xed strain sequential method, we solve the ow problem rst while freezing the stress or strain rate, and then we can solve the mechanics problem with the latest pressure and saturation solution. This procedure is repeated inside a sequential iterative loop until pressure and displacement solutions converge, which then becomes the solution at that time step. The converged solution is identical to the fully coupled solution within the tolerance of convergence. These sequential methods allow re-use of existing ow and mechanics solvers after minor modication to allow for the exchange of frozen quantities 34 from one solver to the other, while keeping the linear system size in each solver the same as in their respective uncoupled version. The computational cost is drastically smaller than a fully coupled method which allows coupled ow and geomechanics sim- ulation of real-world oilelds, aquifers, and geothermal reservoirs over multiple years of injection/production history [48, 99]. However, many such methods address only the coupling between the ow and mechanics problems. Dynamic changes in formation properties such as permeability and stiness due to fracture propagation associated with uid injection and production activities is usually not addressed in these methods. Change in permeability can change the ow eld, and change in Young's modulus can aect the deformation and stress elds. We need to account for this coupling because we want to simulate fracture propagation along with coupled ow and mechanics. 3.1.3 One-way Coupled Sequential Solution Methods The four sequential methods in Figure 3.1 are two-way coupled methods where both uid-to-solid and solid-to- uid couplings [94] are captured. For a weakly coupled sys- tem, for example gas production from a consolidated gas reservoir, solid-to- uid coupling is weaker than the uid-to-solid coupling because the gas compressibility is signicantly larger than that of the grains or the skeleton, and therefore the changes in uid pres- sure or mass due to deformation is negligible. The xed stress or xed strain method described earlier can be modied by xing the number of sequential iterations to one or by dropping the rate of stress/strain term in the mass balance equation. The result- ing sequential solution method is the one-way coupled method (Figure 3.2), which only accounts for uid-to-solid coupling. The communication between pressure and stress is considered insignicant. This means that an independent calculation of the ow eld can be applied. Therefore, the change of stress regime and displacements can be estimated as functions of position and time after determining the pressure from the ow solver. However, when the changes in stress aects the pore pressure behavior, the two-way coupling is signicant. For some problems in ow-induced seismicity, it may be impor- tant to only keep the solid-to- uid coupling, e.g. during undrained deformation of a low permeability basement rock under instantaneous loading, and neglect the uid-to-solid coupling. The resulting method has been called the decoupled method [76]. 35 Figure 3.2: OC: The one-way coupled poromechanics-damage mechanics sequential solution algorithm. Here, C 1 = (A e =t) M 1 + T C dr 1 andC 2 = 0. See text for denitions of other quantities. Time level information (n or n + 1) is shown only for pressure, displacement, stress, and the fracture length. Damage Mechanics Method (DMM) is used to update the poromechanical properties. 36 3.2 Proposed Methods Above, we described dierent types and degrees of coupling methods in the literature. Now, we are ready to present and discuss the solution methods proposed in this work. Although we have implemented and analyzed the results from both fully coupled and iteratively coupled methods, our objective remains to investigate and implement an it- eratively coupled scheme which is stable, convergent and computationally less expensive than the fully coupled approach. we will see next, and in later chapters too, that it is the stability that is more dicult to ensure for an iterative scheme than the economy of the computation. We extend the xed stress method, proposed earlier for isotropic deformation [14, 45, 51, 63], to model the anisotropic deformation behavior observed because of microcracking at the reservoir scale [2, 77] and the core scale [9]. We use the Irwin-Grith theory [34, 40] combined with the continuum damage mechanics (CDM) approach [21, 46, 55, 61, 80] to model fracture propagation and associated changes in mechanical and hydraulic prop- erties of the medium. Recently, the damage mechanics approach has been used to simulate hydraulic fracturing-induced permeability enhancement in geothermal reser- voirs [60, 73] and shale reservoirs [100]. The CDM approach allows one to capture the macroscopic eects{material anisotropy, inelasticity, pressure sensitivity, volumetric di- latation, and strength deterioration{induced by microscopic cracks into continuum scale models through constitutive equations based on a damage variable. The CDM approach is shown to be capable of handling complicated fracture networks with variable fracture volumes and lengths. 3.2.1 Extended Fixed Stress Method with Explicit Damage (CFM) We use the xed stress method to solve the coupled equations of mechanical deformation and ow sequentially. We solve the ow sub-problem rst then we solve the mechanical sub-problem inside a sequential iterative loop. The discretized system of equations for the ow problem is A e t M 1 + T C dr 1 +T n+1;j P n+1;j+1 = A e t M 1 + T C dr 1 n P n A e t ( T C dr 1 ) n+1;j ( T C dr 1 ) n +A e q n+1 w (3.2) where superscript j is the iteration counter of the sequential loop at the new timestep n+1. We discretize the rate of change of stress term in the mass balance equation, which 37 is the solid-to- uid coupling term, using the stress tensor at the previous sequential iteration, n+1;j . In other words, we keep the time rate of volumetric stress, @ v =@t, at its previous sequential iteration value. Hence, it appears on the right hand side of Eq. (3.2). The discretized system of equations for the mechanics problem is K n U n+1;j+1 =F u +Q n P n+1;j+1 (3.3) During the mechanics solve, we keep the pressure xed at its previous value from the current sequential iteration,P n+1;j+1 , which is known because the ow problem is solved before the mechanics problem. Figure 3.3 illustrates the solution algorithm. 3.2.2 Extended Fixed Strain Method The linear system of the extended xed strain method can be derived similar to the xed stress method above: A e t M 1 +T n+1;j P n+1;j+1 = A e t M 1 P n A e t () n+1;j () n +A e q n+1 w (3.4) The linear system for the mechanics problem is the same as for the extended xed stress method. 3.2.3 Extended Fixed Stress Method with Implicit Damage (CFMD) Another sequential solution algorithm that we implemented is presented in Figure 3.4 where damage mechanics is added inside the xed stress loop. Here, the pressure solver incorporates the eect of damage in ow properties (permeability in x and y) and the mechanics solver incorporates the eect of damage in mechanical properties (elastic moduli and Biot parameters). 3.2.4 One-way Coupled Method for Multiphase Geomechanics We use the one-way coupled method (Figure 3.2) to solve our coupled multiphase ow- poromechanics-damage mechanics problem. An implicit nite volume method is used to solve the multiphase ow problem using the MATLAB Reservoir Simulation Toolbox (MRST [59]). The pressures inside the P vector in Figure 3.2 refers to the equivalent pressure p E in Eq. (2.24). The multiphase pressures are used in the solution of the quasi-static mechanical deformation problem using the nite element method. 38 Figure 3.3: CFM: The coupled ow and mechanics sequential solution algorithm. Here, C 1 = (A e =t) M 1 + T C dr 1 and C 2 = (A e =t)( T C dr 1 ). See text for denitions of other quantities. Time level information (n or n + 1) is shown only for pressure, displacement, stress, and the fracture length. Damage Mechanics Method (DMM) is used to update the poromechanical properties. 39 Figure 3.4: CFMD: The coupled poromechanics-damage mechanics sequential solu- tion algorithm with implicit damage. 40 3.3 Implementation The numerical simulator is implemented in MATLAB. Our approach allows the use of a single discretized grid for the numerical solution of the sub-problems. On this grid, the simulator solves the discretized system of equations to determine the vectors of element pressuresP and nodal displacementsU. The element-averaged stresses are obtained by post-processing the displacement solution. After the pressure and displacement solutions have converged within the sequential loop, we accept them as the solution at the current time step and pass them on to the damage mechanics model of the fracture propagation problem. Here, we use the change in the state of stress computed from the ow-mechanics solution to model the changes in permeability and elastic properties via the damage tensor. The fracture length a k is updated by solving the damage-based propagation condition, Eq. (2.33). k,C dr ,, and M are updated using Eqs. (2.51), (2.40), (2.41), and (2.43), respectively. The following steps describe our work ow: 1. Create the 2D domain geometry and mesh. Distribute the rock and uid properties (initial values of porosity, permeability, density, compressibility, and uid viscosity) in the mesh. 2. Fractures are distributed in the mesh with their initial lengths, apertures, and orientations. See below for more details. 3. Dene wells and the simulation schedule. We consider both single-well (injector at the center) and double-well (injector and producer near boundaries) cases. For the injector-producer case, oil is produced at a prescribed rate from a horizontal well located near the left boundary of the domain and the injection well is located near the right boundary. 4. Apply boundary conditions: tectonic compression, xed displacement, and no- ow conditions. 5. Initialize the coupled system using a static step. We assume the reservoir is fully saturated with oil initially. We solve the system for one time step to create the initial state (pressures, saturations, stresses, and displacements) that satises equi- librium equations. 6. Loop through time steps. 41 3.3.1 Fracture Distribution in Space Initially, natural fractures are distributed uniformly in our model with given fracture density, angle, length, and aperture. We classify the cracks into crack families dened using the normal surface vectorn k of the cracks, where k is the crack family index. N families of the idealized line of microcracks are distributed. This initial state is consid- ered the undamaged state with the value of the damage variable as zero. Each single fracture is characterized using the fracture direction ( ), damage surface coecient (!), fracture unit vector (n k ), fracture toughness (K IC ), initial microcrack length (a i ) and connectivity coecients: t 1 = 1 and t 2 = 1. Changes in fracture length and aperture because of ow-induced deformation and rock failure is considered as damage to the hy- dromechanical properties (permeability and elastic moduli), which evolve in space and time as the damage variable evolves with the eective stress. 3.3.2 Fracture Evolution in time Each fracture is treated separately using the equilibrium condition K I K IC . In order to quantify the fracture growth, stress and pressure eld are quantied around the fracture. Fracture growth is incremental and quasi-static. Chapter 4 SIMULATOR VALIDATION In this chapter, we present representative numerical simulations to validate the math- ematical formulation and the computer implementation of our coupled simulator. We present results of mesh renement and time step renement studies to investigate the convergence properties of our formulation. 4.1 Uniaxial Consolidation Here, we discuss the validation of our simulator on Terzaghi's uniaxial consolidation problem in absence of fracture propagation. Figure 4.1 shows a discretized 2D physical domain and a representative 1D section relevant for simulation of the uniaxial consol- idation problem. The analytical solution to the mathematical problem can be taken from [13]. We simulated this problem with our numerical simulator in order to validate the uid-to-solid coupling. Figure 4.2 shows the solution at the last time step near the end of consolidation. The highest pressure drop is at the top where the uid is draining out of the constant pressure boundary and hence the highest displacement is at the top. The below Figures( 4.3 and 4.4) show a good agreement between analytical and our numerical solution at dierent time steps. The displacement prole across the domain at the last time step from the analytical and the simulated solutions are shown in Figure 4.4. There is a good agreement between the two solutions. 4.1.1 Eect of boundary condition We perform two simulations which are dierent only in terms of the mechanical boundary conditions applied on the left and right boundaries of the 2D domain, Figures 4.5 and 4.6. In the rst simulation, we x both horizontal and vertical displacements to zero 42 43 Figure 4.1: Physical model of Terzaghi's uniaxial consolidation problem. 0 10 0 20 40 60 80 100 Displacemenet of x (m) -1 -0.5 0 0.5 1 0 10 0 20 40 60 80 100 Displacemenet of y (m) -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 0 10 0 20 40 60 80 100 Pressure (psi) 1.62 1.64 1.66 1.68 1.7 1.72 1.74 10 5 Figure 4.2: Pressure, x-displacement, and y-displacement from our numerical simu- lation at the last time step. (pinned boundary). In the second simulation, we x the horizontal displacement to zero and vertical traction to zero (roller boundary). In this section, the numerical simulator ware evaluated with dierent boundary conditions as shown in Figures( 4.5 and 4.6), the number of grids both in x and y directions and pore pressure. Two cases will be discussed in chapter 5. 44 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Timestep 200 300 400 500 600 700 800 Pressure (psi) Last node pressure with 6 timesteps Orginal Code Terzaghi Figure 4.3: Comparison of analytical and simulated pressure evolution at the top element. 1 1.5 2 2.5 3 3.5 4 4.5 5 node nubmer -0.016 -0.014 -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 displacement y displacement for last timestep Orginal Code Terzaghi Figure 4.4: Comparison of analytical and simulated displacement prole across the domain at the last time step. 45 Figure 4.5: Experimental Solution: 2D displacement xing all boundaries Figure 4.6: Experimental Solution: 2D displacement xed bottom boundary and rolling left and right boundaries. 4.2 Model Sensitivity We evaluate the sensitivity of our simulation results to discretization errors in time and space. We conduct ve dierent simulations to study the eect of time step renement and eight dierent simulations to study the eect of mesh renement. Both sensitivities are done using our two-way coupled simulator for single-phase uid ow (CFM). For all simulations, each mesh element contains one fracture of initial fracture length of 0:1 m. To quantify the eect, we plot the evolution of fracture length, x-direction permeability, and pressure at an element located 25 m left of the injector. 46 Table 4.1: Simulation runs for the dierent time step range. Simulation Run t d 1 90E 4 2 22E 4 3 9:0E 4 4 6:0E 4 5 4:3E 4 6 2:8E 4 7 2:1E 4 8 1:9E 4 4.2.1 Time Step Renement Eect In this section, we evaluate the eect of rening the simulation time step on fracture propagation and permeability evolution. We create multiple simulation runs with dier- ent time steps and the time step is kept constant within a given simulation. The mesh is kept constant and N x =N y = 64 for all the runs. This gives a mesh element area of A e =2.441 m 2 . The time steps are dened as t d = t f;d n s (4.1) where t f;d is the nal dimensionless time, and n s is the number of time steps in the simulation. For a selected cell, we plot the fracture length, the log-permeability in the x-direction, and the pressure as shown in Figures (4.7, 4.8 and 4.9) for the seven simulations dened in Table 4.1. The solutions converge slowly as the time step is decreased. Analysis of the plots suggest that the time step size should be suciently small to capture fracture growth accurately. 4.2.2 Eect of Mesh Renement In the previous section, we used a mesh size ofN x N y = 64 64 and conducted a time step renement study. Here, we perform a mesh renement study by varying N x (=N y ) as shown in Table 4.2 while keeping the time step constant. The results are shown in Figures( 4.10, 4.11 and 4.12). There is one important dierence in this study compared to the time step renement study. There is always one fracture per element in all the renement cases. This means, for the same volume of rock, the number of fractures increase as the mesh is rened. More fractures means higher eective permeability and lower stiness in that volume. Therefore, the renement cases are not identical to each other in terms of properties and the renement study result will be aected by that. 47 Figure4.7: Time Step Renement Eect: Evolution of the normalized fracture length (a k ) d at a cell with dimensionless time (t d ) for dierent time step sizes. As the time step is decreased, the nonlinear behavior is captured better, and the fracture length solution slowly converges towards the true, unknown solution. Table 4.2: Simulation runs for dierent mesh renement cases. Simulation Run Area of an element A e [m 2 ] Grid size N x N y 1 30 32 32 2 9.8 64 64 3 2.4 100 100 4 1 120 120 5 0.7 150 150 48 Figure 4.8: Time Step Renement Eect: Evolution of the normalized permeability log(k x ) d at a cell with dimensionless time (t d ) for dierent time step sizes. Figure 4.9: Time Step Renement Eect: Evolution of the normalized pressure (p d ) at a cell with dimensionless time (t d ) for dierent time step sizes. 49 Figure 4.10: Mesh Renement Eect: Time evolution of the normalized fracture length (a k ) d for dierent mesh sizes. Figure4.11: Mesh Renement Eect: The normalized permeability log(k x ) d evolution comparison for dierent mesh sizes. 50 Figure 4.12: Mesh Renement Eect: The normalized pressure enhancement (p d ) for dierent mesh sizes. Chapter 5 REPRESENTATIVE NUMERICAL SIMULATIONS In this chapter, we present results from several simulations that are representative of production and injection in naturally fractured petroleum reservoirs. We conduct these simulations using the numerical simulator that we presented earlier. We evaluate the impact of geologic properties e.g. fracture orientation and number of fractures, engi- neering controls e.g. the well location and rate, and geophysical constraints e.g. the boundary stress ratio on changes in hydromechanical properties e.g. permeability and stiness. The level of modeling complexity changes from one simulation case to another. For each case, we discuss the spatial distribution and temporal evolution in properties to understand the impact of ow-induced stresses in the domain. We examine the impact of ow-mechanics-damage coupling on the simulation results by comparing the results from the one-way coupling method (Figure 3.2) with the results from two-way coupling methods (Figures(3.3 and 3.4). 5.1 Two-way Coupled Single-phase Flow-Geomechanics Sim- ulator We present results from multiple simulation cases to verify our computational framework and illustrate the eect of fracturing on deformation and ow. In each case, we analyze the spatial distribution and temporal evolution of pressure (p), stress ( xx and yy ), fracture length (a k ), deviatoric stress (S N ) and permeability (k) in order to understand the impact of ow-mechanics coupling on the damage in poromechanical properties. We evaluate the sensitivity of induced damage to the following model parameters: well ow 51 52 Mathematical symbol Numerical value i 0.1 k i 1 md f 1000 kg/m 3 1 cp p i 2000 psi i 0.21 E i 68 GPa i 0.8 a i 0.1 m ! 0.6 K IC 1.1 MPa p m b 1 m Horizontal compressions, h , H 8.48 MPa Table 5.1: Single-phase poromechanical properties. The subscript i for a property indicates the initial value of that time-dependent property. rate (q w ), fracture orientation ( ), number of fracture families (N), and the boundary stress ratio SR. We will discuss the two-way coupled simulator for single-phase uid ow in two parts. First, we will evaluate the coupling methods and the eect of damage in ow and mechanics properties. Second, we will go through an example of multiple injection rates and dierent fracture networks. 5.1.1 Physical Model The model domain is assumed to be a 2-D horizontal section of a reservoir containing the wellbore as shown in Figure 2.1 adding horizontal injector in the center of the domain. The section is a uid-saturated homogeneous rock embedded with randomly-placed ide- alized line microcracks. The domain is subjected to compressive loads representing the two horizontal stresses on its top and right boundaries, H and h , respectively. The stress ratio is dened as SR = h = H . The bottom and left boundaries are xed in the normal directions. No- ow conditions are prescribed on all boundaries. The model is in hydrostatic and mechanical equilibrium, and the state of the rock before the uid injection is regarded as the initial or undamaged state. The given domain properties already listed in Table 5.2. We consider that each element contains one randomly- oriented fracture. We assume the following values of important physical properties during the simulation. These values are representative of a tight formation. These val- ues imply a principal state of stress initially such that at t = 0, xx = yy , xy = 0, 0 v = 2:55 MPa, andS N = 0 (Eq. (2.30)). We present results in terms of the normalized pore pressure (p d = p=p i ), normalized stress ( d = =p i ), normalized deviatoric stress ((S N ) d =S N =p i ), normalized fracture length (a k;d =a k =a i ), and normalized permeabil- ity (k d =k=k i ). 53 5.1.2 Coupling and Sequential Iteration In uences In this section, we will test the coupling methods and the iteration loop eects on ow, mechanics and damage mechanics problems. Our comparison will use location (1) as an observation location located left of the horizontal injector as shown in Figure 5.1. Figure 5.1: Coupling Methods: Locations of the two fractures selected for detailed investigation. The two fractures are located in the left and right of the injector as shown above. 5.1.2.1 Coupling Methods We compared the three coupling methods described in Figures(3.2, 3.3 and 3.4). Using 30% of formation porosity and injection rate (q w ) of 5 bbl/day , we examined location (1) located left of the horizontal injector shown in Figure 5.1. Figure 5.2 illustrates that each method dier in the outcomes such as pressure, fracture length, deviatoric stress and local total stress. This is mainly because of iteration loop. The most realistic approach is given in Figure 3.4 where we update for ow and mechanical properties within the iteration loop. When ow and mechanical properties are updated inside the iteration loop, the calculation of pressure and stress are simultaneously computed to avoid inaccurate calculation in the next time-step. The estimated pressures in (CFM and CFMD) methods are in perfect match. However, stresses are slightly diered. At CFMD, the deviatoric stress (S N ) shows more tensile eect compared to CFM and OC because of more fractures propagation for CFMD method. For OC, one-way coupling, pressure and stress are not in communication by overestimating the stress. 54 Figure5.2: Coupling Methods: Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. All three coupling techniques are shown dierent behaviors 5.1.2.2 Adaptive Coupling Each sequential coupling method uses a dierent solution path to solve for pressure, stress and permeability variables. The change of pressure aects the stress calculation and ultimately changes the total permeability. Once the uid paths change, we need to go back and calculate the new pressure and stress. The iteration method will loop until these variables converge. The success of a method relies on understanding the formation sensitivity to pressure and stress changes. For example, the uid pressure might change the material volume, or the applied stress might aect the uid pressure calculations. Therefore, we compare the three coupling methods shown in Figure 5.3: (A) one-way coupling (OC), (B) two-way coupling (CFMD) and (C) adaptive coupling (OCFMD). OCFMD method presents the adaptive coupling where we use the two-way coupling method for the rst time step and then the simulator switches to use the one-way coupling method. During the rst time step, we used CFMD framework which allows ow and mechanical properties to be updated inside the iteration loop in which the calculations of pressure and stress are computed sequentially. This allows overcoming any delay in correcting the calculated values in the next time step. Afterward, OC is used. Consequently, OCFMD method is less expensive than the two-way coupled CFMD because the simulation run time of CFMD is longer. We applied the exact same conditions in all methods with an injection rate (q w ) of 5 bbl/day and used similar fractures directions and distributions to examine the fracture propagation. We illustrate 55 Figure 5.3: One-way Coupling (OC) and two-way Coupled Flow Mechanics Damage (CFMD) methods. CFMD uses the Converged decision box to iterate over the ow and mechanics solutions within a sequential loop. OC does not use that. the fracture behavior using four selected locations shown in Figure 5.4. Figure 5.5 and 5.6) illustrate that dierent methods have dierent outcomes. These discrepancies are related to the iteration loop in the two way coupling which allows for feedback from each problem to account for the changes in pressure and stress. During the iteration loop, the ow and mechanics problems are in communication until pressure and displacement are converged. For OCFMD, the iteration occurs only at the rst time step, then the model uses the one-way coupling method. The three methods don't have exact matching. However, the results of OCFMD are close to CFM fracture propagation values. Applying the adaptive coupling (OCFMD) improves the fracture propagation calculation compared to the one-way coupling method (OC). Therefore, OCFMD has the better determination of fracture propagation compared to the OC and it is a less expensive method compared to CFMD. Figure 5.7 shows the pore pressure, the fracture propagation length, the deviatoric stress and the local total stress of location(4). 5.1.2.3 Damage in Flow and Mechanical Properties During each time step, we update ow and mechanical properties using damage mechan- ics ( ow: k(D) and mechanics: C dr (D),(D) and M(D)). Most of the conventional reservoir simulators do not account for updating the mechanical properties changes. These changes signicantly in uenced fracture propagation, reservoir heterogeneity and 56 Figure 5.4: Adaptive Coupling: The four examined fractures locations for the single- phase simulator. also the oil recovery for multiphase simulations that will be discussed in the next section. Our approach can capture these changes. To describe the importance of the damage mechanics, we evaluated two scenarios: (1) Damage in ow eect (F ) and (2) Damage in ow and mechanics (F +M) eect. F models consider the eect of damage on per- meability only. F +M models consider the eect of damage on both permeability and mechanical properties (elasticity tensor). Both models use the CFM coupling method (Figure 3.3). The pressure dierence between the two scenarios (eect damage in ow and eect damage in ow and mechanics) demonstrated in Figure. 5.11. Referring to Eq. (2.52), the mechanical properties now became dynamic and changed with time. This allows pressure to increase more in the second scenario. The eect of the pressure and stress clearly resulted in an increase in permeability as shown in Figure 5.12. 5.1.3 Eect of Biot Coecient We evaluate the single-phase for two Biot coecients ( = 0:1 and = 0:8) using CFMD. From Eq. (2.52), the eect of pore pressure on deformation (linear elastic) is described through the Biot coecient () which describes the deformation in drained 57 Figure 5.5: Adaptive Coupling: The normalized fracture propagation of the three methods for location(1) and location(2). Figure 5.6: Adaptive Coupling: The normalized fracture propagation of the three methods for location(3) and location(4). conditions (i.e., the pore uid is free to leave the rock volume). Biot's theory of poroe- lasticity has two coupling coecients: the Biot modulus(M) and the Biot coecient (), which are scalars for an isotropic medium. They are related to rock and uid properties and are given in Eq. (2.13). Figure 5.19 shows pressure and fracture propagation for both Biot coecients ( = 0:1 and = 0:8). Clearly, the higher Biot coecient leads 58 Figure 5.7: Adaptive Coupling: (a) Normalized pressure p d , (b) Normalized fracture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . These plots show the comparison between the three methods at location (4). Figure 5.8: Single-phase Flow. Eect of Damage at Location(1): Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. 59 Figure 5.9: Single-phase Flow. Eect of Damage at Location(1): (a) Normalized stiness in x and y directions, C dr;ij =E i , (b) normalized Biot coecient in x and y directions, d , (c) normalized Biot modulus, M and (d) normalized permeability, log(kx;d). Figure5.10: Single-phase Flow. Eect of Damage. Snapshots of the permeability (k x ) dierence eld over 8 time steps, log(k x;d ) = log(k (x;d;F+M) ) log(k (x;d;F) ). F +M indicates coupled ow and mechanics. Higher permeability dierence (log(k x;d )) in north and south area of the domain. 60 Figure5.11: Single-phase Eect of Damage: Snapshots of the pressure dierence eld over 8 time steps, p d = p (d;F+M) p (d;F) . Pressure is higher with F +M in north and south of the domain compared to the F case. Figure 5.12: Single-phase Eect of Damage: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N =S (N;F+M) S (N;F) . More stress acts in north and south area. 61 Figure 5.13: Single-phase Eect of Damage: Comparison of Lorenz Coecient (LC) between damage in ow and damage in ow and mechanics properties for permeability (k x ) in x-direction. Figure 5.14: Single-phase Eect of Damage: Damage D yy evolution in y-direction. The damage D xx equals to zero. 62 Figure5.15: Single-phase Eect of Damage: Snapshots of the stress ( xx ) for damage in ow case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. Figure5.16: Single-phase Eect of Damage: Snapshots of the stress ( yy ) for damage in ow case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. 63 Figure5.17: Single-phase Eect of Damage: Snapshots of the stress ( xx ) for damage in ow and mechanics case eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. Figure 5.18: Single-phase Eect of Damage: Snapshots of the stress ( xx ) for dam- age in ow and mechanics eld over 8 time steps. The creation of anisotropy and heterogeneity in the stress eld increased with time. 64 Figure 5.19: Single-phase - Biot Coecient Eect: Pressure (p d ) and Fracture length (a k;d ) enhancement processes for dierent Biot coecients ( = 0:1 and = 0:8). to slightly higher pressure but strongly higher fracture propagation. Therefore, Biots coecient is a key parameter that quanties the contribution of pore pressure to the eective stress and dened in Eq. (2.4) and Eq. (2.52). Figure 5.20 illustrates the dif- ference of permeabilities between both cases and be described in Figure 5.12. Because of increase in permeability for = 0:8, the reservoir heterogeneity increases as shown in Figure 5.23. 5.1.4 Eect of Young's Modulus Young's Modulus E is a mechanical property which measures the ability of a material to withstand changes in length under tension or compression loads. That means in the second case (higher Youngs Modulus) the material is stier and resists to expansion or shrinkage. The drained bulk modulus is directly proportional to Young's modulus, K dr = E 2(12) . Because of the drained modulus term in the solid-to- uid coupling term in the pressure equation, the pore pressure increases causing a higher fracture propagation compared to the low Youngs Modulus case as shown in Figure 5.24. In Eq. (2.13), we can observe the eect of Young's Modulus. Therefore, the increase ofK dr results in an increase of pore pressure described in Figure 5.26. As a result, more uid diuses into the formation because of the increase of the total permeability. Because the fractures are oriented in thex-direction, the uid moves in this direction causing the 65 Figure 5.20: Single-phase - Biot Coecient Eect: Snapshots of the permeability (k x ) dierence eld over 8 time steps, log(k x;d ) = log(k (x;d;0:8 ) log(k (x;d;0:1) ).Higher permeability dierence (log(k x;d )) in north and south area of the domain for = 0:8 Figure5.21: Single-phase - Biot Coecient Eect: p d =p (d;0:1) p (d;0:8) . Snapshots of the pressure dierence eld over 8 time steps. Pressure diuses more of damage ow and mechanics case in north and south of the domain. More pressure accumulates in that case. 66 Figure 5.22: Single-phase - Biot Coecient Eect: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N =S (N;0:8) S (N;0:1) . More stress acts in north and south area. Figure 5.23: Single-phase - Biot Coecient Eect: Biot Coecient Eect: Com- parison of Lorenz Coecient (LC) between damage in ow and damage in ow and mechanics properties for permeability (k x ) in x-direction. 67 Figure 5.24: Single-phase - Young's Modulus Eect: Pressure (p d ) and Fracture length (a k;d ) enhancement processes for two Young's Modulus (E = 68 GPa and E = 136 GPa) . fracture to propagate and grow perpendicular to the horizontal injector. In this case, the eect of tensile load in the fracture resulting from the pressure increase is shown in Figure 5.25. In Figure 5.28, both cases show an increase of formation heterogeneity and the second case shows more heterogeneity. 5.1.5 Case I: Injection Well at the Center of the Domain (Base Case) We distribute the fractures in the domain with orientation distributed randomly as shown in Figure 5.29. This distribution corresponds to all fractures oriented in the general northwest-southeast direction of the domain and belonging to the same fracture familyk = 1. We inject water at 50 bbl/day for 88 days uniformly along the 100 ft length of the horizontal injector. Figure 5.30 shows the diusion of pressure over three months. At the beginning (t = 1 day), the pressure diuses uniformly around the injector because of an isotropic distribution of the poromechanical properties. At t = 0, the stress eld is homogeneous isotropic with xx = yy =8:48 MPa. As pressure increases shown in gure 5.30, the near wellbore region expands away from the well and applies compression on the far region. Figure 5.31 shows the distribution of stress eld at two time steps. The eect of in Eq. (2.31) leads to anisotropy in the stress and permeability elds 68 Figure 5.25: Single-phase - Young's Modulus Eect: Snapshots of the dierence in k x permeability eld over 8 time steps. The permeability dierence is dened as log(k x;d ) = log(k (x;d;E=136) )log(k (x;d;E=68) ). Higher permeability dierence in north and south area of the domain for E = 136 GPa. Figure 5.26: Single-phase - Young's Modulus Eect: Snapshots of the pressure dif- ference eld over 8 time steps, p d = p (d;E=68) p (d;E=136) . Pressure diuses more uniformly and higher for the case of E = 68 GPa. 69 Figure 5.27: Single-phase - Young's Modulus Eect: Snapshots of the deviatoric stress (S N ) dierence eld over 8 time steps, S N =S (N;2E) S (N;E) . More stress acts in north and south area. Figure 5.28: Single-phase - Young's Modulus Eect: Comparison of Lorenz Coe- cient (LC) between Youngs Modulus (E = 68 GPa andE = 136 GPa) for permeability (k x ) in x-direction. 70 Figure 5.29: Case I: Distribution of the initial fracture direction (t = 0) in the domain. All the fractures are oriented along the general northwest-southeast direction. All the fractures belong to one single family i.e. N = 1, and each mesh element contains one fracture i.e. m k = 1=A e . Figure 5.30: Case I: Snapshots of the pressure eld over approximately three months with a time step of t = 11 day between the snapshots. Pressure diuses uniformly for the rst three time steps such that the pressure contours are elliptic in shape with the well lying at the center along the north-south long axis of the ellipse. Later the contours start to orient along the NW-SE direction following the permeability enhancement in that direction because of injection-induced fracture growth. between NE-SW and NW-SE corridors. This leads to preferential ow along the NW- SE corridor and anisotropy in the pressure contours, which orient along the NW-SE corridor (as shown in Figure 5.30). Stress anisotropy ( xx 6= yy ) and rotation of the principal stress directions from the x and y axes ( xy 6=0) causes the magnitude of S N to increase from zero explained in Eq. (2.31). The actual increase depends on both the quadrant position and the fracture direction . For example, for a REV location in the NE corridor with =45 , 1 =60 , 2 =30 , and =135 ,S N =p w rL 2 w =(2(r 1 r 2 ) 3=2 ), which is positive. For a location in the NW corner, with similar values of the angles, S N becomes negative. The fractures start to propagate in regions with high S N . The heterogeneity in introduces heterogeneity in xx , yy , xy , S N , a k , k andC dr . The 71 Figure 5.31: Case I: stresses (MPa) at two dierent time steps showing the creation of anisotropy and heterogeneity in the stress eld, which is homogeneous isotropic at t = 0. Anisotropy and shear, resulting from the eect of in Eq. (2.19), leads to an increase of the deviatoric stress S and its crack-normal projection S N , which aects crack propagation through Eq. (2.34). The heterogeneity in the stress eld results from the heterogeneous distribution of , which makes S N and subsequently k and C dr heterogeneous. Figure 5.32: Case-I: the change in the permeability eld k max at the same 8 time steps as in Figure 5. The change in permeability increases in the beginning (changing color from blue to red along the NW-SE corridor) and decreases later (red to yellow) as the induced stresses relax because of pressure diusion and fracture growth processes. permeability growth eld is shown in Figure 5.32. Heterogeneity and anisotropy of the permeability and elasticity tensor elds enter the pressure solution through the solid- to- uid coupling term and the rock compressibility term as shown in Eq. (2.52). The pressure eld is relatively smooth because of the diusive nature of the mass balance equation. We select two locations in the domain: Location (1) in the NW corner and Location (2) in the NE corner (as shown in Figure 5.33) to conduct a more detailed 72 Figure 5.33: Case-I: Locations of the two fractures selected for detailed investigation. The two fractures are located in NW and NE corners of the domain as shown by the circles. The fracture lengths are dierent in space because the plots are at the end of Case I simulation. analysis. The volumetric compression v increases monotonically at both locations because of near-wellbore expansion under injection. However, deviatoric tension S N behaves dierently at the two locations because ofn < 0 in NE and > 0 in NW corners (Eq. (2.19)). Figure 5.34 shows that as the fracture grows at a location, stiness in the x direction decreases and the Biot coecient increases. This suggests a loss of poroelastic strength and an increase in the drained compressibility of the medium. 73 -8.44 -8.4 -8.36 < v [ P a ] #10 6 -2 0 2 S N [ P a ] #10 5 0 0.02 0.04 0.06 0.08 t d 1.2 1.3 1.4 1.5 C d r ; 1 1 1 1 = E o 0 0.02 0.04 0.06 0.08 t d 0.79 0.8 0.81 , 1 1 1.3 1.4 1.5 C d r ; 1 1 1 1 = E o 1 2 3 k m a x = k o Loc (1) Loc (2) Time (a) (b) (c) (d) Time Time Figure 5.34: Case I: (a) The deviatoric stress S N is positive (tensile) at Loc (2) and negative (compressive) at Loc (1) because of the dierence in signs of , 1 , and 2 in Eq. (2.31). At both locations,S N evolves towards zero as the stresses relax with fracture propagation. The volumetric stress v at both locations are negative (compressive) because of near-wellbore expansion. This is also evident from the negative sign of v in the rst equation of Eq. (2.19) for both positive and negative signs of . (b) Normalized stiness in the x direction, C dr;1111 =E 0 , decreases with time because of fracture growth, (c) 11 increases with time suggesting reduction in the drained bulk modulus, (d) normalized permeability k max =k 0 vs. normalized stiness C dr;1111 =E 0 curves show the two implications of fracturing: increase in permeability and decrease in stiness. 5.1.6 Case II: Eect of Variable Injection Rate In many oilelds, the injection well ow rate varies with time because of various opera- tional constraints e.g. availability of injection pumps, decrease in injectivity because of calcite scale build-up on the wellbore wall, and controlling the breakthrough of injected uid at producers. It is also a common practice to slowly ramp up the injection rate in a new injection well to avoid unintentional fracturing or activation of thief zones around the well. We want to evaluate the eect of variability in the injection rate on fracture propagation and permeability evolution. To that eect, we compare results from con- stant rate and variable rate injection simulations. The nal cumulative injected volume and the injection duration are the same between the two cases. Figure 5.35 shows the cumulative injection prole of the two cases. Figure 5.36 illustrates that, similar to 74 Figure 5.35: Case I and II: cumulative injection prole in the constant and variable rate injection cases. The variable rate well starts o at an injection rate (slope of the dotted curve above) lower than that of the constant rate well. Later it gradually increases and surpasses the constant rate. The total injected volume is the same in both cases. Case I, the Case II pressure eld develops anisotropy with time as a result of anisotropic fracture growth. However, the two pressure elds are very dierent in terms of the pressure values. This dierence results from the dierence in the injected mass at any given time and the dierence in fracture propagation with time. One observation is that the near-wellbore pressure is signicantly higher in Case II. This is related to a slower increase in fracture length and permeability because of a lower injection rate in the beginning of the variable rate case{compare Figure 5.37 with Figure 5.32. Constant Figure 5.36: Case II: evolution of pressure over three months with a time step of t = 11 day during variable rate injection. Compare with Case I above. rate injection has resulted in higher permeability values. Because of a slower fracture 75 growth in Case II, which is discussed below with Figure 5.38, water starts to accumulate around the wellbore, which leads to larger near-wellbore pressures. Figure 5.37: Case II: Changes in the permeability eld, k max , at dierent times with a constant time step of t = 11 day. Compare with Case I above. In Figure 5.38, we analyze the eect of rate variability on fracture propagation and damage at the two selected locations of Figure 5.33. By denition, fracture propagation 0 0.02 0.04 0.06 0.08 0.1 0.2 0.3 0.4 CR: Loc(1) CR: Loc(2) VR: Loc(1) VR: Loc(2) 0 0.02 0.04 0.06 0.08 0.2 0.4 0.6 0.8 0 0.02 0.04 0.06 0.08 2000 2500 3000 3500 4000 4500 0 0.02 0.04 0.06 0.08 -4 -2 0 2 4 10 5 Figure 5.38: Case II: fracture length (a k ) and permeability (k max ) enhancement processes are dependent on the location and the well injection rate. CR and VR denote constant and variable rate Cases 1 and 2, respectively. Loc (1) and Loc (2) denote the NW and NE corner locations selected earlier. is a direction-dependent process and depends on local stresses. At both locations, the fracture length a k for Case I is larger than that for Case II, which is mimicked by the 76 maximum permeabilityk max plot. An interesting observation is that the evolution in the location pressure does not mimic the evolution in a k or k max exactly. For example, a k for Case I Location (2) is lower thana k for Case II Location (1). A lower fracture length or permeability should result into a higher pressure. However, p for Case I Location (2) is lower than p for Case II Location (1). The reason is that the pressure equation is diusive and has a global character. Because of the permeability around wellbore is lower in Case II than in Case I, the pressure is higher in Case II at all locations. The S N plot in Figure 5.38 provides further insights into the coupling between fracture propagation and stress evolution processes. We observe a non-monotonic trend in S N in Case I, which is lost in Case II over the chosen duration of the simulation because of its slower injection. The eect of location i.e. on S N is twofold. First, < 0 at Location (2) makesS N positive (tensile) as per Eq. (2.19). Second, a smaller pressure p at Location (2) compared to Location (1) makes S N smaller in magnitude. This causes a slower and smaller growth in the fracture length a k at Location (2) (as shown in Figure 5.38). 5.1.7 Case III: Injecting at Constant rate for Dierent Times We studied the rate scheduling eects in the fracture propagation process. In real oil eld operation activities, the injection rates are not continues because of dierent rea- sons such as the water supply, the water treatment and the water handling facility limitations. For these reasons, we applied dierent water injection rates schedule shown in Figure 5.39 to evaluate the fracture prorogation. Dierent fracture locations were evaluated for the single horizontal injector as shown in Figure 5.4. During the injection period, we applied more compression in the domain and changed the local stress distri- bution that allow fractures to propagate. After we stopped the injection, we observed that the fracture continued to propagate until pressure stabilized. Clearly, this means that the fracture propagation is aected by injection rate behavior. Therefore, we moni- tored the normalized propagated length (a n;d ), dened asa n;d =a n+1 k;d a n k;d , as shown in Figure 5.40 to evaluate the behaviors of fracture propagation during and post-injection periods. Figure 5.41 illustrates the pressure behavior for dierent locations in the do- main. This behavior was also observed in the propagated length of the fracture (a n;d ). We also demonstrate the pressure and the stresses behaviors in Figure 5.42. In these plots, one can observe that the pressure and stress control the fracture propagation. For example, after the pressure reached to the peak point when we stopped injecting and the pressure started to drop, the slope of the fracture propagation began to change as shown in the fracture propagation plot. 77 Figure 5.39: Case III: the cumulative injected water volume (bbl) from the single horizontal well. Figure 5.40: Case III: The normalized fracture propagated length (a n;d ) for loca- tion(3) and location(4) with time. 78 Figure 5.41: Case III: the pressure behaviors at the injector and the observation location. Figure 5.42: Case III: (a) Normalized pressure p d , (b) Normalized fracture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . Two selected locations were investigated. Location (3) is located to the left of the horizontal injector and Location (4) in the right. 79 5.1.8 Eect of Fracture Orientation The fracture direction plays an important role in selecting the natural fractures that propagate during injection because it determines the magnitude of t at the fracture tip. This was mentioned in chapter 3 through -dependent quantitiesS N and 0 N (Eq. (2.29)) and also indicated through Figure 5.33. Therefore, we hypothesize that aects the evolution and distribution of poromechanical properties during injection. To test our hypothesis over a larger range of fracture directions, we conduct multiple simulations by varying at the two locations across the following range: = 90 ; 110 ; 135 ; 150 ; 180 , while keeping xed at all other locations. We analyze the simulation results in Fig- ure 5.43 for Location (1) and Figure 5.44 for Location (2). We observe that the growth in permeability is a function of 0 v and the deviatoric tension S N as expected from Eq. (2.30). Figure 5.44 shows that, at Location (2), the magnitude of S N is smaller compared to the values at Location (1). This explains the smaller permeability ob- served at Location (2) for all . 0 1 2 3 2 4 6 8 =180 o 150 o 135 o 110 o 90 o 0 1 2 3 -4 -2 0 10 5 0 1 2 3 0.9 0.95 1 1.05 1.1 -4 -2 0 10 5 2 4 6 8 Figure 5.43: Single-phase - Eect of Fracture Orientation: the eect of the fracture direction on evolution in the deviatoric tensionS N , the local stress ratio xx = yy , and the maximum permeability k max at Location (1) at early times. S N is non-monotonic with because of the sin 2 term in Eq. (2.30). Growth in k max is proportional to the growth inS N magnitude. k max vs. S N shows that as the magnitude of deviatoric stress increases, the permeability increases for all . 80 0 1 2 3 1 1.5 2 2.5 3 =180 o 150 o 135 o 110 o 90 o 0 1 2 3 -2 -1 0 1 2 10 5 0 1 2 3 0.9 0.95 1 1.05 1.1 -2 0 2 10 5 1 1.5 2 2.5 3 Figure 5.44: Single-phase - Eect of Fracture Orientation: the eect of the fracture direction on evolution in the deviatoric tensionS N , the local stress ratio xx = yy , and the maximum permeability k max at Location (2) at early times. k max is proportional to S N magnitude. Figure 5.45: Single-phase - Eect of Fracture Orientation ( = 0 ): Evolution of pressure over three months with a time step of t = 11 day. Compare with Case I above . 81 Figure 5.46: Single-phase - Eect of Fracture Orientation ( = 90 ): Evolution of pressure over three months with a time step of t = 11 day. Compare with Case I above . Figure 5.47: Single-phase - Eect of Fracture Orientation ( = 135 ): Evolution of pressure over three months with a time step of t = 11 day. Compare with Case I above . Also, we have evaluated the uniform fracture orientations for dierent fracture direction. We used same fracture direction in all REV and illustrated the pressure diusivity in the domain as shown in Figures(5.45, 5.46 and 5.47). We selected three fracture orientations : 0 , 90 and 45 . We used the sequential method described in Figure 3.4 for all scenarios. In Figure 5.48, we used location(A) to compare all fracture orientations. 82 Figure 5.48: Single-phase - Eect of Fracture Orientation - Location (A): Pressure (p) and Fracture length (a k ) enhancement processes are dependent on the location and the well injection rate. 5.1.9 Multiple Fracture Families In this case, we consider a random distribution of the fracture direction over the entire range from 0 to 180 as shown in Figure 5.49. This can be understood as a case with multiple fracture families with N 2. Eq. (2.51) predicts that the change in permeability will be larger in this case because of the increase in the total number of fractures, which creates additional connectivity. The injection rate is the same as in Case I. In Figure 5.50, we see the evolution of the permeability eld. Compared to Case I with the fracture family k = 1, we see that the permeability enhancement is azimuthally more uniform. This is expected from the eect of on S N . The result is that the Stimulated Reservoir Volume (SRV) is larger and spatially more uniform than in Case I. This highlights the importance of ascertaining the fracture direction in the reservoir region away from the injector. 83 2.5 3 3.5 4 4.5 5 5.5 6 F r e q u e n c y #10 -3 0 50 100 150 F r a c t u r e D i r e c t i o n A ( d e g ) 0 50 100 150 200 250 300 Nu m b e r of f r a c t u r e s Figure5.49: Distribution of the fracture direction for the two fracture families present in the domain. The fracture orientations are randomly distributed. Figure 5.50: Multiple Fracture Families for single-phase: Snapshots of the change in permeability eld for a fracture family with fracture directions distributed over the entire range of 0-180 . 84 5.1.10 Boundary Stress Anisotropy Boundary stresses can signicantly impact the onset, type and magnitude of failure in rocks subjected to hydraulic stimulation because they determine the solution of the equilibrium equation Eq. (2.73). We focus on the eect of anisotropy in the boundary stresses on permeability growth and poroelastic damage. We conduct four simulations with the following values of the boundary stress ratio, SR = 1, 1/2, 1/3, and 1/4 obtained by increasing the boundary stress in the y-direction, H . Figure 5.51 shows the eect of the stress anisotropy on permeability growth at the two selected locations in the domain. At Location (1) (> 0), as SR decreases from 1, i.e. the stress anisotropy increases, the permeability growth decreases and the deviatoric tension increases such that it become tensile at high anisotropy. At Location (2), as the stress anisotropy increases, both the permeability growth and the deviatoric tension increase. This is an interesting result that cannot be deduced from the permeability vs. deviatoric tension behavior at the isotropic case (SR = 1) shown in Figure 5.43 and 5.44. For isotropic case, the permeability grew proportional to the S N magnitude. 0.2 0.4 0.6 0.8 1 S R 1 2 3 4 5 6 7 k m a x =k o 0 5 10 15 S N [ P a ] #10 5 Location (1) t d = 0 : 04 0. 03 0. 02 0. 01 0.2 0.4 0.6 0.8 1 S R 1 2 3 4 5 6 7 k m ax =k o 0 5 10 15 S N [ P a ] #10 5 Location (2) k m a x S N k m ax S N Figure 5.51: Boundary Stress Anisotropy for single-phase: Evolution of the normal- ized permeability k max =k 0 and the deviatoric tension S N with the stress ratio (SR) at the two locations. The plot shows that as SR increases i.e. the boundary stress anisotropy decreases the permeability increases at Location (1) and decreases at Loca- tion (2). S N decreases with SR at both locations. 85 5.2 One-way Coupled Simulator for Multiphase Flow and Geomechanics 5.2.1 Comparison with Single-phase Simulator In this section, we used the one-way coupling simulator for the single-phase uid ow as our base case. In Figure 3.2, MRST uid ow solver replaced the single-phase ow problem and this coupling method is purposed in method(4). However, this new coupling method (method (4)) has to obtain a good match with method(1). The ow coupling equation was discussed earlier in Eq. (2.52) and the uid mass balance equation was addressed in Eq. (2.21) by keeping mechanical and damage mechanics problems constant. In Figure 5.52, we compared pore pressure (p), the fracture length (a k ), the deviatoric stress (S N ) and the horizontal stresses ( xx and yy ) for location (1). We obtained perfect match for the pore pressure. At the last time step, the pressure dierence between the two methods is around (2 psi). The fracture length, the deviatoric stress and the horizontal stresses also resulted in a reasonable match. The small dierence in pore pressure is mainly because MRST uses a dierent approach in updating formation and uid properties while using the damage mechanics to update for permeability growth. Figure 5.52: Compare one-way coupled simulators for single-phase and multiphase (MRST) codes at location one. 86 Figure 5.53: One-way coupling multiphase simulator : Pressure enhancement pro- cesses for 8 time steps. 5.2.2 Physical Model of an Injector-Producer Pair After bench-marking the multiphase model, we used the model to simulate the oil and the water system by injecting water into an oil reservoir in order to improve the oil recovery and maintain the reservoir pressure. As shown in Figure 5.54, we assumed a 2D model domain containing two wells horizontally placed near the left and right boundaries paralleled to the y-direction. The domain is assumed to be saturated with oil and embedded with randomly-placed idealized line microcracks. The horizontal sec- tions of these wells are around 20 m long. The domain is subjected to compressive loads representing the two horizontal stresses on its top and right boundaries, H and h , re- spectively. The bottom and left boundaries are xed in the normal directions. No- ow conditions are prescribed on all boundaries. The model is in hydrostatic and mechan- ical equilibrium, and the state of the rock before the uid injection is regarded as the initial or undamaged state Table 5.2 provides all poromechancial properties. Initially, we distribute the fractures in the domain with an orientation of = 0 (180 ). 87 Figure 5.54: Left: Two horizontal wells injector and producer placed in a natu- rally fractured reservoir containing microcracks and larger fractures. Right: A two- dimensional conceptual model of the two horizontal wells cross-section of the fractured reservoir showing the injector and producer intersecting multiple pre-existing microc- racks. The model domain containing the randomly distributed, non-interacting micro- cracks is considered to be undamaged before the well becomes active. Mathematical symbol Numerical value 0.1 k 0 10 md w 1000 kg/m 3 o 700 kg/m 3 w 1 cp o 5 cp p 0 2000 psi 0 0.21 E 0 68 GPa 0 0.8 a 0 0.1 m ! 0.6 K IC 1.1 MPa p m b 1 m q o 5 bbl/day q w 5 bbl/day Table 5.2: Multiphase poromechanical properties. 5.2.3 Eect of Damage in Flow properties In this part, we will use damage in ow properties scenario as the base case with oil production rate (q o ) of 1 bbl/day and water injection rate (q w ) of 5 bbl/day. All fractures are oriented with zero degree ( = 0 ). We will use damage in ow properties scheme to update the permeability (k). All other mechanical properties are held constant. Here, we will compare two models: Model(1) includes damage mechanics in which fractures propagate with time and Model(2) considers that fractures don't propagate with time. 88 Figure 5.55: Multiphase - Damage in Flow properties Eect: (a) Pressure (p d ) and Fracture length (a k;d ) evolution are dependent on the location. We will be able to demonstrate the importance of modelling natural fractures during injection and production activities. After water breakthrough (WBT), the oil well started to produce both oil and water presented as a total uid production rate (q f ). The simulation time for this case is 853 days (using 100 time steps). At t=0, all stress elds in both directions are homoge- neous, isotropic and equal. Although we are withdrawing from the domain, the pressure around the injector rapidly increases because of the high injection rate as shown in Fig- ure 5.57. As result, more compression is applied around the wellbore causing fracture to propagates. Therefore, the permeability started to increase in the x-direction because fractures are oriented in the same direction and parallel to the well placement. We notice that most of the permeability evolution occurs in North and South areas of the domain which is related to the deviatoric stress changes (Figure 5.61). This is mainly at early times. For comparison purposes, we will use the oil saturation map as shown in Fig- ure 5.62 to evaluate the cumulative oil production prole illustrated in Figure 5.56). The nal cumulative oil production for both cases are 0.52 Pore Volume (PV) for Model(1) and 0.5 PV for Model(2). The dierence of 0.02 PV is due to fracture propagation and the total permeability. Model(1) allows the water to move into the formation displacing more oil to the producer, which delays the water production, compared to Model(2). See Figure 5.62. 89 Figure 5.56: Multiphase - Damage in Flow properties Eect: the upper two gures illustrate water saturation (%) and permeability change (logk x ) d ) at the middle be- tween injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. Figure5.57: Multiphase - Damage in Flow properties Eect: The normalized pressure (p d ) enhancement at nine time steps. 90 Figure 5.58: Multiphase - Damage in Flow properties Eect: The normalized perme- ability (logk x ) d ) growth in x-direction when fractures are perpendicular to horizontal wells . There is no fracture propagation in y-direction. Figure 5.59: Multiphase - Damage in Flow properties Eect: The normalized stress distribution ( xx ) d in x-direction at nine time steps. 91 Figure 5.60: Multiphase - Damage in Flow properties Eect: The normalized stress distribution ( yy ) d in y-direction at nine time steps. Figure 5.61: Multiphase - Damage in Flow properties Eect: The normalized devia- toric stress (S N ) d at nine time steps. 92 Figure 5.62: Multiphase - Damage in Flow properties Eect: the oil saturation (S o ) at nine time steps. Figure 5.63: Multiphase - Damage in Flow properties Eect: the fractures propaga- tion distribution for the last time step (t d = 100). 93 Figure 5.64: Multiphase - Damage in Flow properties Eect: Comparison of Lorenz Coecient (LC) between Model(1) and Model(2)for permeability (k x ) in x-direction. 5.2.3.1 Relative Permeability Eect We assume two-phase ow with oil and water as the immiscible phases. The relative permeability curves are shown in Figure 5.65. We compared the results from two dierent relative permeability cases to evaluate the eect on fracture propagation process and uid movements. Because of the variation of phase velocities, the results showed that each case has dierent water breakthrough (WBT) illustrated in Figure 5.66. Both cases started to produce water at the same time (similar WBT); however, after a period of time, the water production proles started to deviate and this dierence increased over time. One can explain this behavior using Darcy's equation: v ` = kk r ` ` rp ` : We can observe that the higher relative permeability of water (k r w ) allows water to move faster. Using our approach and observation, we nd that the fracture propagation process also plays a role in water production from naturally fractured reservoir. The early fracture propagation allows more water to move toward the producer. At location (1), Figure 5.68 shows the comparison between the two cases. In the second case, we obtained earlier fracture propagation causing a higher water production as illustrated in Figure 5.66. Also, the dierence in water relative permeability curves (k r w ) in uenced the early water breakthrough and became more signicant when the water saturation is larger than 50% as shown in Figure 5.65. When the water, the less compressible uid, moves away from the injector displacing oil, it causes the pore pressure and deviatoric stress to increase (Figure 5.68). This increase allows fractures to 94 Figure 5.65: Multiphase - Relative Permeability Eect: Relative permeability curve for the second case. The residual oil saturation is 20% and the initial water saturation is 20%. Figure 5.66: Multiphase - Relative Permeability Eect: watercut at the producer for both cases. 95 Figure5.67: Multiphase - Relative Permeability Eect for Location(1): The examined fractures location for multiphase ow-geomechanics simulator. propagate. Consequently, the lower viscosity uid moves faster causing a higher water production. In conclusion, the recoverable oil is aected by the water movement and the fracture propagation process in naturally fractured reservoirs. Our observation showed that the process of fracture propagation in terms of the propagation time and the fracture length aects uid movements and consequently might in uence the amount of the produced oil in heterogeneous formations. The compression eects are illustrated in Figures(5.68 and 5.67). Although the pressure dierence is very small, permeability growth for the second case is higher as illustrated in Figure 5.71. This is because of the eects of deviatoric stress (S N ) as explained in Figure 5.67. As we mentioned before, the compressibility aects hydraulic diusion time scale and causes more fracture propagation. 5.2.3.2 Fluid Viscosity Eect Usually, the eect of formation uid properties on fracture propagation is neglected because of the Darcy ow conditions and stronger eects of other properties. We tested dierent oil viscosity ( o ) and compared the sweep eciency and the recoverable oil, assuming that w < o . Studying various oil viscosities allows us to examine the eect 96 Figure 5.68: Multiphase - Relative Permeability Eect: (a) Normalized pressure p d , (b) Normalized fracture length (a k ) d , (c) Normalized deviatoric stress (S N ) d and Normalized stress d . The plots illustrate the comparison between the two cases for the middle fracture. This fracture is located between the injector and producer. Figure 5.69: Multiphase - Relative Permeability Eect: the upper two gures illus- trate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumu- lative oil production (bbl) at the producer. 97 Figure 5.70: Multiphase - Relative Permeability Eect: The dierence in normalized pressure at nine time steps. p d =p d;Relperm1 p d;Relperm2 . Figure 5.71: Multiphase - Relative Permeability Eect: The normalized perme- ability dierence in x-direction when fractures are perpendicular to horizontal wells. (log(k x )) d = log(k (x) ) d;Relperm1 log(k (x) ) d;Relperm2 . There is no fracture propaga- tion in y-direction. 98 Figure 5.72: Multiphase - Relative Permeability Eect: The normalized stress distri- bution ( xx ) d dierence inx-direction at nine time steps, ( xx ) d = ( xx ) d;Relperm1 ( xx ) d;Relperm2 . Figure 5.73: Multiphase - Relative Permeability Eect: The normalized stress distri- bution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;Relperm1 ( yy ) d;Relperm2 . 99 Figure5.74: Multiphase - Relative Permeability Eect: The normalized the deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;Relperm1 (S N ) d;Relperm2 . Figure 5.75: Multiphase - Relative Permeability Eect for Relperm(2): the oil satu- ration (S o ) at nine time steps. 100 Figure 5.76: Multiphase - Relative Permeability Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) Relperm(1) (S o ) Relperm(2) . Figure 5.77: Multiphase - Relative Permeability Eect: the fractures propagation distribution for the last time step (t d = 100). 101 Figure 5.78: Multiphase - Relative Permeability Eect: Comparison of Lorenz Coef- cient (LC) between Relperm(2) and Relperm(2) for permeability (k x ) in x-direction. of the uid mobility on the fracture process. We compared o = 2 cp case with the base case ( o = 5 cp). The less viscous uid is more mobile than the more viscous uid. Also, oil is more compressible than water. For this reason, we observe lower pressures as per Eq. (2.2) and Eq. (2.52) and shown in Figure 5.81. Figure 5.87 shows that water displaced more oil in o = 2 compared to the base case. The water ngering phenomenon is also stronger at o = 5. For the fracture growth, the fracture starts to propagate early in the base case (Figure 5.79) because water (less compressible uid) arrives faster to the middle location. When we monitor the dierences in pressure (p d ) and deviatoric stress ((S N ) d ) over time, we observe that fracture is gradually propagating. This allowed more oil to be produced. We conclude that the fracture process plays an important role in the oil recovery from stress-sensitive, fractured reservoirs. 5.2.4 Eect of Mechanical Properties Many reservoir engineering studies solve only the ow problem and focus on modeling the spatial heterogeneity in permeability and porosity elds due to geological and petro- physical structure, which are considered stationary in time. However, the recent growth in advanced reservoir modeling has sparked an interest in the coupled ow-geomechanical models that attempts to model the evolution in rock mechanical properties as a result of stress change. This is especially required in formations containing pre-existing fractures, 102 Figure 5.79: Multiphase - Fluid Viscosity Eect: (a) Normalized pressure (p d ) and fracture length (a k;d ) evolution depend on the location. Figure 5.80: Multiphase - Fluid Viscosity Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. 103 Figure 5.81: Multiphase - Fluid Viscosity Eect: The normalized pressure (p d ) dif- ference at nine time steps. p d =p d;=5 p d;=2 . Figure 5.82: Multiphase - Fluid Viscosity Eect: The normalized permeability (logk x ) d dierence in x-direction when fractures are perpendicular to horizontal wells (log(k x )) d = log(k (x) ) d;=5 log(k (x) ) d;=2 . There is no fracture propagation in y- direction. 104 Figure 5.83: Multiphase - Fluid Viscosity Eect: The normalized stress distribution ( xx ) d dierence in x-direction at nine time steps, ( yy ) d = ( yy ) d;=5 ( yy ) d;=2 . Figure 5.84: Multiphase - Fluid Viscosity Eect: The normalized stress distribution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;=5 ( yy ) d;=2 . 105 Figure 5.85: Multiphase - Fluid Viscosity Eect: The normalized deviatoric stress (S N ) d dierence at nine time steps, (S N ) d = (S N ) d;=5 (S N ) d;=2 . Figure 5.86: Multiphase - Fluid Viscosity Eect for o = 2: the oil saturation (S o ) at nine time steps. 106 Figure 5.87: Multiphase - Fluid Viscosity Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) =2 (S o ) 0=2 . Figure 5.88: Multiphase - Fluid Viscosity Eect: the fractures propagation distribu- tion for the last time step (t d = 100). 107 Figure 5.89: Multiphase - Fluid Viscosity Eect: Comparison of Lorenz Coecient (LC) between o = 5 [cp] and o = 2 [cp] for permeability (k x ) in x-direction. either natural or hydraulically stimulated. The cumulative oil recovery and water break- through proles must strongly depend on the fracture evolution process, which means geomechanical modeling must be a part of the history matching exercise. Below we investigate the sensitivity of oil recovery and water breakthrough proles to variation in three mechanical properties: Biot modulus, Young's modulus and Biot coecient, which appear dierently in the ow and mechanics equations and therefore are expected to causes dierent sensitivities. We used the two Biot's modulus (M): 100 MPa for the base case and 150 MPa. We xed = 0:8. We see that the increase of Biot's modulus has delayed the water breakthrough by t d = 0:14 as shown in Figure 5.90. In Figures 5.90 and 5.91, we analyzed results at an element in the middle between injector and producer. We can estimate the eect of M on pressure using Eq. (2.52). Figure 5.92 shows the dierence between k x from the two cases. 5.2.4.1 Young's Modulus Eect In order to consider dierent types of sedimentary rocks in dierent parts of the world and at dierent depths, we vary the Young's modulusE of the porous medium and study its eect on fracture propagation and oil recovery. Dierent values of E cause dierent strain for the same value of the ow-induced stress near the wellbore, which then aects 108 Figure 5.90: Multiphase - Biot Modulus Eect: (a) pore pressure p, (b) permeability in x-direction k x , (c) water saturation S w and (d) watercut (WC) at producer. Figure 5.91: Multiphase - Biot's Modulus Eect: the pressure eld at 15 dierent time steps for M = 150 MPa. 109 Figure 5.92: Multiphase - Biot Modulus Eect: the dierence in x-direction per- meability when fractures are perpendicular to horizontal wells ( = 0). log(k x ) = log(k (x) ) M=150 log(k (x) ) M=100 . There is no fracture propagation in y-direction. the stress far from the wellbore and the fracture propagation process. Also, the fracture opening decreases as the stiness of the rock matrix increases, and consequently the fracture permeability decreases. We consider two values of E, uniformly distributed in space, in two simulations with the same initial and boundary conditions. We monitor ow, mechanical and fracture properties over 100 time steps. In multiphase ow systems, the water saturation distribution is aected directly by the permeability distribution because it controls the water-oil displacement process at a local scale. Figure 5.94 shows the evolution of permeability, where evolution curves at an element in the middle between the two wells are shown. At early times, the permeability growth and pressure evolution curves are identical for the two cases. At t d = 25, the permeability curves start to deviate. Figure 5.93 shows that the normalized deviatoric stress (S N ) d , which drives the fracture propagation process, behaves dierently. The faster arrival of less compressible uid (water) to the middle point results in less recoverable oil at the producer. 5.2.4.2 Biot Coecient Eect The second mechanical property for which we conduct a sensitivity study is the Biot coecient. Biot's theory of poroelasticity has two coupling coecients: the Biot mod- ulus M and the Biot coecient , which are scalars for an isotropic medium. They are 110 Figure 5.93: Multiphase - Young's Modulus Eect: (a) Evolution in normalized pressure (p d ) and fracture length (a k;d ) depends on the location. Figure 5.94: Multiphase - Young's Modulus Eect: the upper two gures illustrate water saturation (%) and permeability change at an element at the center of the domain between the injector and the producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. 111 Figure 5.95: Multiphase - Young's Modulus Eect: The dierence in normalized pressure at nine time steps. p d =p d;E1 p d;E2 . Figure 5.96: Multiphase - Young's Modulus Eect: The normalized permeability (logk x ) d dierence in x-direction when fractures are perpendicular to horizontal wells (log(k x )) d = log(k x ) d;E1 log(k x ) d;E2 . There is no fracture propagation iny-direction. 112 Figure 5.97: Multiphase - Young's Modulus Eect: Multiphase - Young's Modulus Eect: The normalized stress distribution ( xx ) d in x-direction at nine time steps, ( xx ) d = ( xx ) d;E1 ( xx ) d;E2 . Figure5.98: Multiphase - Young's Modulus Eect: The normalized stress distribution ( yy ) d in y-direction at nine time steps, ( yy ) d = ( yy ) d;E1 ( yy ) d;E2 . 113 Figure 5.99: Multiphase - Young's Modulus Eect: The normalized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;E1 (S N ) d;E2 . Figure 5.100: Multiphase - Young's Modulus Eect for E = 168 MPa: the oil satu- ration (S o ) at nine time steps. 114 Figure 5.101: Multiphase - Young's Modulus Eect: the oil saturation dierence at nine time steps. S o = (S o ) E1 (S o ) E2 . Figure 5.102: Multiphase - Young's Modulus Eect: Comparison of Lorenz Coef- cient (LC) between the two cases of dierent Young's modulus, E1 and E2. The coecient is calculated using permeability in x-direction, k x . 115 Figure 5.103: Multiphase - Biot Coecient eect: (a) Normalized pressure (p d ) and fracture length (a k;d ) evolution depend on the location. related to rock and uid properties as per Eq. (2.13). For multiphase ow system, we used = 0:1 and = 0:8 as the two sensitivity cases. Eq. (2.52) shows the eect of the Biot coecient on pore pressure evolution and spatial distribution. It is evident from the presence of in the system compressibility term multiplying the time derivative of pressure. It also determines the strength of solid-to- uid coupling term multiplying the time derivative of volumetric stress. Biot coecient aects total stress and the fracture length and therefore the permeability. So, the uid mobility term in Eq. (2.52) multiply- ing the pressure laplacian is also aected, although it is an indirect eect. Figure 5.103 compares the two cases. With a higher value of Biot coecient, the system compress- ibility goes up and we observe a lower pressure trend. The fractures grow longer, and the cumulative oil recovered at the producer is higher. 5.2.5 Eect of Damage in Flow and Mechanics Properties We continued the work to evaluate the eect of mechanical properties. In this case, we allow the model to update for the changes in the mechanical properties using the damage mechanics method. The changes in both ow and mechanical properties are 116 Figure 5.104: Multiphase - Biot Coecient eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. Figure 5.105: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized pressure at nine time steps. p d =p d;=0 p d;=0:1 . 117 Figure 5.106: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized permeability inx-direction when fractures are perpendicular to horizontal wells = 0. (log(k x )) d = log(k (x) ) d;=0:8 log(k (x) ) d;=0:1 . There is no fracture propagation in y-direction. Figure 5.107: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8) : The normal- ized stress distribution ( xx ) d dierence in x-direction at nine time steps, ( xx ) d = ( xx ) d;o=0:8 ( xx ) d;o=0:1 . 118 Figure 5.108: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8) : The normal- ized stress distribution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;o=0:8 ( yy ) d;o=0:1 . Figure 5.109: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): The dierence in normalized deviatoric stress at nine time steps, (S N ) d = (S N ) d;=0:8 (S N ) d;=0:1 . 119 Figure 5.110: Multiphase - Biot Coecient eect for = 0:1 : oil saturation (S o ) at nine time steps. Figure 5.111: Multiphase - Biot Coecient eect ( = 0:1 and = 0:8): the oil satura- tion dierence at nine time steps, (S o ) = (S o ) o=0:8 (S o ) o=0:1 . 120 Figure 5.112: Multiphase - Biot Coecient eect for = 0:1 : fractures propagation distribution for the last time step (t d = 100). Figure 5.113: Multiphase - Biot Coecient eect: Comparison of Lorenz Coecient (LC) between o = 0:8 and o = 0:1 for permeability (k x ) in x-direction. 121 compared with the previous case. Fracture length (a k ) is signicantly increased and this is because of the change of (C dr ; and M). Equations (Eq. (2.40), Eq. (2.41) and Eq. (2.43)) show how we can update for these mechanical properties. These changes can re ect the impact on recovery. We use the water breakthrough (WBT) prole as a comparison tool. Our observation identied that faster water breakthrough occurred for both cases because of the increase of permeability. However, there are slight delayed in the water breakthrough byt d = 0:02. More importantly, the formation heterogeneity increases in this case by incorporating the changes of the mechanics properties case as shown in gure (5.124). Figure 5.114: Multiphase - Damage in Flow and Mechanics Properties Eect: (a) Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. All three coupling techniques are shown dierent behaviors. 5.2.6 Eect of Well Placement In order to evaluate the eect of well position with respect to the fracture orientation on fracture propagation and eld recovery, we vary the fracture angle instead of changing the location of the two wells. This was discussed with one well in the single-phase simulator Section 5.1.8. Here, it is dierent because we have two wells (producer and injector). Our model aims to describe the eect using the water breakthrough (WBT) prole. We consider that both horizontal wells are placed parallel to each others. The simulations 122 Figure 5.115: Multiphase - Damage in Flow and Mechanics Properties Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. Figure 5.116: Multiphase - Damage in Flow and Mechanics Properties Eect: The normalized pressure (p d ) dierence at nine time steps, p d =p d;F p d;F+M . 123 Figure 5.117: Multiphase - Damage in Flow and Mechanics Properties Eect: The normalized permeability (logk x ) d growth inx-direction when fractures are perpendicu- lar to horizontal wells (log(k x )) d = log(k (x) ) d;F log(k (x) ) d;F+M .There is no fracture propagation in y-direction. Figure 5.118: Multiphase - Damage in Flow and Mechanics Properties Eect: The normalized stress distribution ( xx ) d dierence in x-direction at nine time steps, ( xx ) d = ( xx ) d;F ( xx ) d;F+M . 124 Figure 5.119: Multiphase - Damage in Flow and Mechanics Properties Eect: The normalized stress distribution ( yy ) d dierence in y-direction at nine time steps, ( yy ) d = ( yy ) d;F ( yy ) d;F+M . Figure 5.120: Multiphase - Damage in Flow and Mechanics Properties Eect: The normalized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d;F (S N ) d;F+M . 125 Figure 5.121: Multiphase - Damage in Flow and Mechanics Properties Eect: oil saturation (S o ) at nine time steps. Figure 5.122: Multiphase - Damage in Flow and Mechanics Properties Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) F (S o ) F+M . 126 Figure 5.123: Multiphase - Damage in Flow and Mechanics Properties Eect: the fracture distribution at the last time step (t d = 100). Figure 5.124: Multiphase - Damage in Flow and Mechanics Properties Eect: Com- parison of Lorenz Coecient between two cases: Damage in ow property D(F) and Damage in ow and mechanics properties D(F+M). The coecient is calculated based on permeability in x-direction, k x , only. = 0 . 127 are carried out for the water injection up to 100 days at three fracture orientations: = 0 (perpendicular to horizontal wells), = 90 (parallel to horizontal wells) and = 135 . Each case leads to a dierent fracture propagation behavior, which directly aects the permeability growth as dened in Eq. (2.51). Figure 5.126 describes the water saturation, the recoverable oil and the permeability growth for each fracture orientation. The fastest breakthrough is expected when the fractures are perpendicular to horizontal wells ( = 0 ), which is what we see in our result, Figure 5.126. Fluids ows in the direction of fracture illustrated in the oil saturation maps as shown in Figures 5.137 and 5.138. Figure 5.125: Multiphase - Well Placement Eect : (a) Pressure (p d ) and Fracture length (a k;d ) enhancement processes are dependent on the location. The three coupling methods show dierent behaviors. 5.2.7 Eect of Multiple Fracture Families In this case, we consider a random distribution of the fracture direction over the entire range from 0 to 180 as shown in Figure 5.49. This can be understood as a case with multiple fracture families withN 2. Eq. (2.51) predicts that the change in permeabil- ity will be larger in this case because of the increase in the total number of fractures, which creates additional connectivity. Dierent families of fracture distribution lead to dierent results in terms of pressure (as shown in Figure 5.145), uid sweep (as shown in 128 Figure 5.126: Multiphase - Well Placement Eect: the upper two gures illustrate water saturation (%) and permeability change (log(k x;d ) at the middle between injection and producer. The bottom two gures show water saturation (%) and cumulative oil production (bbl) at the producer. Figure 5.127: Multiphase - Well Placement Eect ( = 0 and = 135): The dier- ence in normalized pressure dierence at nine time steps. p d =p d; =0 p d; =135 . 129 Figure5.128: Multiphase - Well Placement Eect ( = 0 and = 90): The dierence in normalized pressure dierence at nine time steps. p d =p d; =0 p d; =90 . Figure 5.129: Multiphase - Well Placement Eect ( = 0 and = 135): The nor- malized permeability (logk y ) d in y-direction. 130 Figure 5.130: Multiphase - Well Placement Eect for = 90: The normalized per- meability (logk y ) d in y-direction when fractures are parallel to horizontal wells. There is no fracture propagation in x-direction. Figure 5.131: Multiphase - Well Placement Eect ( = 0 and = 135): The dif- ference in normalized stress in x-direction at nine time steps. ( xx ) d = ( xx ) d; =0 ( xx ) d; =135 . 131 Figure5.132: Multiphase - Well Placement Eect ( = 0 and = 90): The change in normalized stress inx-direction at nine time steps. ( xx ) d = ( xx ) d; =0 ( xx ) d; =90 . Figure 5.133: Multiphase - Well Placement Eect ( = 0 and = 135): The dif- ference in normalized stress in y-direction at nine time steps. ( yy ) d = ( yy ) d; =0 ( yy ) d; =135 . 132 Figure 5.134: Multiphase - Well Placement Eect ( = 0 and = 90): The dier- ence in normalized stress in y-direction at nine time steps. ( yy ) d = ( yy ) d; =0 ( yy ) d; =90 . Figure 5.135: Multiphase - Well Placement Eect ( = 0 and = 135): The dierence in normalized deviatoric stress (S N ) d at nine time steps, (S N ) d = (S N ) d; =0 (S N ) d; =135 . 133 Figure5.136: Multiphase - Well Placement Eect ( = 0 and = 90): The dierence in normalized deviatoric stress ((S N ) d ) at nine time steps, (S N ) d = (S N ) d; =0 (S N ) d; =90 . Figure 5.137: Multiphase - Well Placement Eect for = 135: the oil saturation (S o ) at nine time steps. 134 Figure 5.138: Multiphase - Well Placement Eect for = 90: the oil saturation (S o ) at nine time steps. Figure 5.139: Multiphase - Well Placement Eect ( = 0 and = 135): the change in oil saturation at nine time steps, S o = (S o ) =0 (S o ) =135 . 135 Figure 5.140: Multiphase - Well Placement Eect ( = 0 and = 90): the change in oil saturation at nine time steps, S o = (S o ) =0 (S o ) =90 . Figure 5.141: Multiphase - Well Placement Eect for = 135: fracture distribution for the last time step (t d = 100). 136 Figure 5.142: Multiphase - Well Placement Eect for = 90: fracture distribution for the last time step (t d = 100). Figure 5.144) and stress regime (as shown in Figures(5.146, 5.147 and 5.148)). We com- pared the two fracture networks by subtracting the two oil saturations for each time step illustrated in Figure 5.144. The gure shows that the oil distribution for each fracture distribution is dierent. The upper and lower parts of the domain have less oil movement for the base case whereas the middle has higher oil movement. These discrepancies are mainly linked back to the permeability enhancement and fracture propagation for each case. 137 Figure 5.143: Multiphase - Fracture Network Eect: Distribution of the fracture di- rection for family (k=3) present in the domain. The fracture orientations are randomly distributed. Figure 5.144: Multiphase - Fracture Network Eect: the oil saturation dierence at nine time steps, (S o ) = (S o ) F1 (S o ) F2 . 138 Figure 5.145: Multiphase - Fracture Network Eect: The normalized pressure dier- ences (p d ) at nine time steps, p d =p d;F1 p d;F2 . Figure 5.146: Multiphase - Fracture Network Eect: The normalized stress distribu- tion dierences in x-direction at nine time steps, ( xx ) d = ( xx ) d;F1 ( xx ) d;F2 . 139 Figure 5.147: Multiphase - Fracture Network Eect: The normalized stress distribu- tion dierences in y-direction at nine time steps, ( yy ) d = ( xx ) d;F1 ( yy ) d;F2 . Figure5.148: Multiphase - Fracture Network Eect: The normalized deviatoric stress (S N ) d dierence at nine time steps, (S N ) d = (S N ) d;F1 (S N ) d;F2 . 140 Figure 5.149: Multiphase - Fracture Network Eect (k=3): The normalized perme- ability (logk x ) d growth in x-direction. Figure 5.150: Multiphase - Fracture Network Eect (k=1): The normalized perme- ability (logk y ) d growth in y-direction. 141 Figure 5.151: Multiphase - Fracture Network Eect (k=3): The normalized perme- ability (logk x ) d growth in x-direction. Figure 5.152: Multiphase - Fracture Network Eect (k=3): The normalized perme- ability (logk y ) d growth in y-direction. 142 Figure 5.153: Multiphase - Fracture Network Eect (k=1): fractures propagation distribution for the last time step. Figure 5.154: Multiphase - Fracture Network Eect (k=3): fractures propagation distribution for the last time step. 143 Figure 5.155: Multiphase - Fracture Network Eect (k=1 and k=3): Comparison of Lorenz Coecient (LC) between the two families for both permeabilities (k x and k y ). Chapter 6 CONCLUSION AND FUTURE WORK 6.1 Conclusion We studied the problem of coupled single-phase/multiphase ow, poroelastic mechanics and fracture propagation (damage mechanics) in naturally fractured porous formations. A two-dimensional horizontal section of a reservoir containing pre-existing fracture line was successfully modeled using our MATLAB code. We used two coupling methods: the one-way coupling (OC) for multiphase ow and the two-way coupling (CFM and CFMD) for the single-phase ow. Our approach does not require the computational mesh to conform to existing or induced fractures. The continuum damage mechanics theory is used to quantify the changes of poromechanical properties and to identify the formation anisotropy by dening the damage tensor in terms of the geometry and properties of propagating fractures. We presented a novel numerical framework to model fracturing-induced changes in poromechanical propertiesthe elasticity tensor, the permeability tensor, the Biot co- ecient tensor, the Biot modulus, and the drained bulk modulusin a naturally fractured formation. The framework couples Biots poroelasticity theory, the Grith failure the- ory, and the continuum damage mechanics theory to model propagation of fractures with randomly distributed angles. The numerical coupling method extends the xed stress method to model the dynamic anisotropy in the properties. A sequential solution is used to solve the three coupled problems: (1) single-phase/multiphase ow, (2) mechanical deformation, and (3) fracture propagation. An implicit nite volume method is used to solve the single-phase ow and multiphase ow problems. Following a one-way coupled 144 145 scheme, multiphase pressures are used in the solution of the quasi-static mechanical de- formation problem using the nite element method. This is followed by the solution of the fracture propagation problem and associated evolution in rock properties using the damage mechanics method. We used a continuum damage mechanics approach to update the permeability and poroelastic properties of the rock. Because of the growth of fractures and associated anisotropic damage, it is necessary to introduce the anisotropic elasticity tensor using new fracture lengths and the dimensionless scalar damage variables. The change of damage tensor aects both ow and mechanical properties and can be used to update the formation permeability (k), the drained elasticity tensor (C dr ), the Biot coecient tensor (), the drained bulk modulus (K dr ), and the Biot modulus (M). One of the distinguished outcomes using the proposed model is the eect of coupling levels on the fracture propagation lengths. The results show that the two-way coupling method (CFM or CFMD) gives a larger fracture propagation length compared to the one-way coupling method (OC) because of the linkage among the three problems. CFM and CFMD allow for pressure, stress and other properties to get the feedback during the same time step. We also discussed the adaptive coupling (OCFMD) and showed enhanced results with less expensive simulation runs. The damage eect was evaluated to estimate the in uence of mechanical properties. When these properties are constant, lower values for fracture propagation were observed. As the fractures grow, formation stiness decreases and the Biot coecient () increases, which results in losing poroelastic strength and an increase in the drained compressibility (k dr )of the medium. Several numerical examples of water injection were compared to analytical solutions and the available experimental results. Results from multiple representative simulations are presented to validate the proposed framework and to investigate the eect of injec- tion/production rate variability, the initial distribution of fracture direction, and in-situ stress anisotropy on evolution in the properties. The role of deviatoric stress in property evolution is analyzed in detail. For multiple water injection rates, we observed that the fracture propagation is controlled by the pressure behaviors as explained from the crack aperture equation. For multiphase uid ow simulator, parameters of the uid viscosity and relative per- meability functions are assumed to be dierent for the matrix and fracture. Although uid mobilities are not directly correlated with fracture propagation, the propagation process is in uenced. Our model was able to estimate the water breakthrough (WBT) 146 for each scenario. The results showed that the fracturing process plays a signicant role in the water breakthrough where a high oil viscosity expedites the fracture process. This was described using pressure and deviatoric stress behaviors. The early fracture propagation with the high oil viscosity resulted from the compression stress and devi- atoric stress. However, with low oil viscosity, the fracture propagates because of the increase of pressure acting as a tensile load. We also examined the relative permeabil- ity functions eect against the water breakthrough. Additional to the eects on the dierence of water and oil velocities for the recoverable oil in certain saturations, the fracture propagation process also changes due to the change of increasing water veloc- ity toward the producers. The one-way coupling method of the multiphase uid ow simulator was also tested against dierent values mechanical properties. The results showed that the more consolidated rocks result in higher fracture propagation, whereas the higher Young's modulus (E) results in lower fracture propagation for the same stress and pressure conditions. It is important to accurately place horizontal wells in naturally fractured reservoirs. We evaluated dierent angles between fracture and well trajectory. We examined three fracture orientations: = 0 o , = 90 o and = 135 o assuming that each grid has a single fracture and a similar orientation. Obviously, fracture propagation length and process vary for each case in which the water moves in the direction for the fracture propagation. It is necessary to mention for a real oil reservoir that formations do not have a uniform fracture distribution. Realistically the classication of fractures in-terms of orientation, length and family (k) need to be modeled in our simulators. We applied dierent fracture families by randomly distributed fracture orientations. The dierent uid sweep was observed for each in which they can be tied back to the permeability growth and the fracture orientations distribution. The results show that this framework can be used to optimize the injection/production schedule rates, eectively place horizontal wells, monitor uid sweep and estimate the stimulated reservoir volume during induced fracturing of naturally fractured rocks. 6.2 Recommendation There are many interesting studies that can extend the work of this dissertation. These studies include (1) multiphase ow applications in petroleum industry to model three phases (oil, water and gas) in three-dimensional model under subsurface conditions in order to optimize production and recovery, and (2) applications in structural geology to estimate the direction of natural fractures and complex fracture network. 147 In this study, the cracks are assumed to be open under remote compressive stresses in the subsurface and they propagate because of the tensile eects, the volumetric normal traction and the deviatoric stress. Some experimental studies [60, 74] have been per- formed on open fractures under compressive loadings coupled with damage mechanics method; however, the results are limited to xed strain and single-phase ow. Fur- thermore, laboratory experiments in pre-existing cracks under compressive loadings are necessary to expand the applicability in reservoir scale. In addition, studies [82] have examined the mode required to fracture growth and represented in the mathematical formulation. There are many interesting studies that can be extended to the work in this dissertation. These studies include: 1. Two-way Coupling Methods in the Multiphase Flow Simulator. The two-way coupling allows for all three models to communicate. The model will cap- ture the formation changes considering for the feedback from each method in which stress changes aect fractures propagation and these propagation increases formation permeability that in uences uid ow in the same iteration loop. 2. Eect of Formation Wettablility. Under wettablility conditions when the countercurrent imbibition is signicant the cap- illary forces must be modeled with respect to matrix fracture interaction. Also, the dierences in the end points saturation must be modeled properly. 3. Three Phase Flow. Extending our coupled simulator to handle three phases (oil, water and gas systems) in a three-dimensional domain. 4. Induce Fracture in Gas Injection. One of the advanced technologies is evaluating CO2 injection in naturally fractured reservoirs. The existing convectional simulators do not account for the eect of uid compressibility in fracture creation and propagation. 5. Three-dimensional Mapping of Fracture Formation. Inverse modeling method to estimate the development of induced-fracture networks from permeability growth and injection/production activities. Understanding of distributed fracture network (DFN) is one of the main challenges for reservoir geology. Many re- searchers have presented the concept of production geology. The work focuses in de- veloping the fractures properties distribution considering the depositional background. Statistical techniques develop formation properties in the unknown areas by running 148 enormous realizations. This can be applied in estimating the DFN. Running dierent DFN realizations and monitoring injection and production rates, our coupled simulator provides the permeability growth due to activities. Applying the inverse modeling tech- niques, we can obtain the best estimation of DFN of the fractured formation. 6. Fracture Interaction Modeling Using Damage Mechanics Method Extending the crack model to include fracture interaction after fracture reaches the crit- ical length. Our work can be compared with experimental studies where the damage constants and fracture toughness can be updated due to fracture propagation. References [1] B. T. Aagaard, M. G. Knepley, and C. A. Williams. A domain decomposition approach to implementing fault slip in nite-element models of quasi-static and dynamic crustal deformation. J. Geophys. Res., 118:3059{3079, 2013. [2] O. Al-Harrasi, A. Al-Anboori, J-M. Kendall, and A. Wuestefeld. Seismic anisotropy in a hydrocarbon eld estimated from microseismic data. Geophys. Prosp., 59:227{243, 2010. doi: 10.1111/j.1365-2478.2010.00915.x. [3] Y. F. Alghalandis. ADFNE: Open source software for discrete fracture network engineering, two and three dimensional applications. Comput. Geosci., 102:1{11, 2017. [4] P.M. Boerrigter A.P.G. Van Heel and J.J. Van Dorp. Thermal and hydraulic matrix-fracture interaction in dual- permeability simulation. speree, 11(4):735{ 749, 2008. [5] N. Hosseinia A.R. Khoeia and T. Mohammadnejad. Numerical modeling of two- phase uid ow in deformable fractured porous media using the extended nite element method and an equivalent continuum model. Advances in Water Re- sources, 94:510{528, 2016. [6] K. Aziz and A. Settari. Petroleum Reservoir Simulation. Elsevier, London, 1979. [7] G.I. Barenblatt. The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics, 7:55{129, 1962. [8] J. Bear. Dynamics of Fluids in Porous Media. Wiley, New York, 1972. [9] P. Benson, A. Schubne, S. Vinciguerra, C. Trovato, P. Meredith, and R. P. Young. Modeling the permeability evolution of microcracked rocks from elastic wave ve- locity inversion at elevated isostatic pressure. J. Geophys. Res., 111:B04202, 2006. doi: 10.1029/2005JB003710. [10] M. A. Biot. General theory of three-dimensional consolidation. J.Appl.Phys., 12: 155{164, 1941. 149 References 150 [11] M.A. Biot. General theory of three-dimensional consolidation. J. Appl. Phys., 12: 64{155, 1941. [12] E. Sifakis J. Hellrung C. L. Richardson1, J. Hegemann and J. M. Teran. An xfem method for modelling geometrically elaborate crack propagation in brittle materials. IJNME, 2(2):1{41, 2009. [13] R. D. Carter and G. W. Tracy. An improved method for calculating water in ux. Petrol. Trans., AIME, 219:415{417, 1960. [14] N. Castelletto, J. A. White, and H. A. Tchelepi. Accuracy and convergence prop- erties of the xed-stress iterative solution of two-way coupled poromechanics. Int. J. Numer. Anal. Methods Geomech., 2015. doi: 10.1002/nag.2400. [15] A. H. D. Cheng. Material coecients of anisotropic poroelasticity. Int. J. Rock Mech. Min. Sci., 34:199{205, 1997. [16] L. S. Costin. A Microcrack model for the deformation and failure of brittle rock. J. Geophys. Res., 88:85{92, 1983. [17] O. Coussy. Mechanics of Porous Continua. John Wiley and Sons, Chichester, England, 1995. [18] O. Coussy. Poromechanics. John Wiley and Sons, Chichester, England, 2004. [19] O. Coussy, R. Eymard, and T. Lassabat ere. Constitutive modeling of unsaturated drying deformable materials. J. Eng. Mech., 124(6):658{557, 1998. [20] R. de Borst, M. Criseld, and C.Verhoosel. Damage mechanics. In: Nonlinear nite element analysis of solids and structures. John Wiley and Sons Ltd, Chich- ester, England, 2012. [21] R. de Borst, M. Criseld, and C. Verhoosel. Damage mechanics. In: Nonlinear nite element analysis of solids and structures. John Wiley and Sons Ltd, Chich- ester, England, 2012. [22] R. H. Dean, X. Gai, C. M. Stone, and S. E. Minko. A comparison of techniques for coupling porous ow and geomechanics. Soc. Pet. Eng. J., 11:132{140, 2006. [23] A. Dragon and Z. Mroz. A continuum model for plastic-brittle behaviour of rock and concrete. Int. J. Eng. Sci., 17:121{137, 1979. [24] K.F. Du and G. Stewart. Transient pressure response of horizontal wells in layered and naturally fractured reservoirs with dual-porosity behavior. In SPE Annual Technical Conference and Exhibition, number SPE 24682, Washington, 1992. References 151 [25] D. Dugdale. Yielding of steel sheets containing slits. J Mech Phys Solids, 8: 100{104, 1960. [26] M. T. Ebrahimi, D. Dini, D.S. Balint, A. P. Sutton, and S. Ozbayraktar. Dis- crete crack dynamics: A planar model of crack propagation and crack-inclusion interactions in brittle materials. Int. J. Solids Structures, 152:12{27, 2018. doi: 10.1016/j.ijsolstr.2018.02.036. [27] E. Erdogan and G. C. Sih. On the crack extension in plates under plane loading and transverse shear. J. Basic Eng., 85:519{525, 1963. [28] F. Erdogan and G. Sih. On the crack extension in plates under plane loading and transverse shear. Fluids Eng, 85:125 { 519, 1963. [29] A. Al-Ghamdiand I. Ershaghi. Pressure transient analysis of dually fractured reservoirs. SPEJ, 1(1):93{100, 1996. [30] K. L. Feigl, D. C. Agnew, Y. Bock, D. Dong, A. Donnellan, B. H. Hager, T. A. Herring, D. D. Jackson, T. H. Jordan, R. W. King, S. Larsen, K. M. Larson, M. H. Murray, Z. K. Shen, and F. H. Webb. Space geodetic measurement of crustal deformation in central and southern California, 1984{1992. J. Geophys. Res., 98:677{712, 1993. [31] Xiuli Gai. A coupled geomechanics and reservoir ow model on parallel computers. Ph.D. Dissertation, University of Texas at Austin, 2004. [32] E. Gallyamov, T. Garipov, D. Voskov, and P. Van den Hoek. Discrete fracture model for simulating water ooding processes under fracturing conditions. Int. J. Numer. Anal. Methods Geomech., 42:1445{1470, 2018. [33] A. Golshani, Y. Okui, M. Oda, and T. Takemura. A micromechanical model for brittle failure of rock and its relation to crack growth observed in triaxial compression tests of granite. Mech Mater, 38:287{303, 2006. [34] A. A. Grith. The phenomena of rupture and ow in solids. Phil. Trans. roy. soc. London Ser. A, 221:163{198, 1921. [35] A.A. Grith. The phenomena of rupture and ow in solids. Philos. Trans. R. Soc. Lond. Ser., A 22:163 { 198, 1921. [36] A. Hillerborg, M. Modeer, and P. Petersson. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and nite elements. Cem Concr Res, 6:773{382, 1976. References 152 [37] F. Homand-Etienne, D. Hoxha, and JF. Shao. A continuum damage constitutive law for brittle rocks. Comput Geotech, 22:135{151, 1998. [38] B.X. Hu and H. Huang. Stochastic analysis of reactive solute transport in het- erogeneous, fractured porous media: a dual-permeability approach. TPM, 48(1): 1{39, 2002. [39] M. Hussain, S. Pu, and J. Underwood. Strain energy release rate for a crack under combined mode i and mode ii. in: Fracture analysis. Math. Models Methods Appl., 1:2{28, 1973. [40] G. R. Irwin. Analysis of stresses and strains near the end of a crack traversing a plate. J. Appl. Mech., 29:361{364, 1957. [41] J. C. Jaeger and N. G. W. Cook. Fundamentals of Rock Mechanics. Chapman and Hall, London, 1979. [42] Y. Jalali and I. Ershaghi. Pressure transient analysis of heterogeneous naturally fractured reservoirs. In SPE California regional meeting, number SPE 16341-MS, California, 1987. [43] and R.W.Zimmerman J.C.Jaeger, N.G.W.Cook. Fundamentals of Rock Mechanics. Blackwell, Australia, 2007. [44] B. Jha and R. Juanes. A locally conservative nite element framework for the simulation of coupled ow and reservoir geomechanics. Acta Geotechnica, 2:139{ 153, 2007. [45] B. Jha and R. Juanes. Coupled multiphase ow and poromechanics: A computa- tional model of pore pressure eects on fault slip and earthquake triggering. Water Resour. Res., 50(5):3776{3808, 2014. doi: 10.1002/2013WR015175. [46] J.Mazars and G. Pijaudier-Cabot. From damage to fracture mechanics and con- versely: a combined approach. Int. J. Solids Structures, 33:3327{3342, 1996. [47] R. Juanes, E. J. Spiteri, F. M. Orr, Jr., and M. J. Blunt. Impact of relative per- meability hysteresis on geological CO 2 storage. Water Resour. Res., 42:W12418, 2006. [48] R. Juanes, B. Jha, B. H. Hager, J. H. Shaw, A. Plesch, L. Astiz, J. H. Dieterich, and C. Frohlich. Were the May 2012 Emilia-Romagna earthquakes induced? A coupled ow-geomechanics modeling assessment. Geophys. Res. Lett., 2016. doi: 10.1002/2016GL069284. References 153 [49] P. Kattan. MATALB Guide to Finite Element. Springer, Amman, Jordan, second edition, 2008. [50] J. Kim, H. A. Tchelepi, and R. Juanes. Stability and convergence of sequen- tial methods for coupled ow and geomechanics: Drained and undrained splits. Comput. Meth. Appl. Mech. Eng., 200:2094{2116, 2011. [51] J. Kim, H. A. Tchelepi, and R. Juanes. Stability and convergence of sequential methods for coupled ow and geomechanics: Fixed-stress and xed-strain splits. Comput. Meth. Appl. Mech. Eng., 200:1591{1606, 2011. [52] J. Kim, H. A. Tchelepi, and R. Juanes. Rigorous coupling of geomechanics and multiphase ow with strong capillarity. Soc. Pet. Eng. J., 18(6):1123{1139, 2013. [53] R. Kregtin. Cohesive zone models towards a robust implementation of irreversible behaviour. Philips Applied Technologies, 2005. [54] B. Lawn. Fracture of Brittle Solids. Cambridge University Press, Cambridge, 2nd edition, 1993. [55] J. Lemaitre. How to use damage mechanics. J Mech Phys Solids, 80:233{245, 1984. [56] R. W. Lewis and B. A. Schre er. The Finite Element Method in the Static and Dynamic deformation and Consolidation of Porous Media. Wiley, Chichester, England, second edition, 1998. [57] R. W. Lewis and Y. Sukirman. Finite element modelling of three-phase ow in deforming saturated oil reservoirs. Int. J. Numer. Anal. Methods Geomech., 17: 577{598, 1993. [58] R. W. Lewis and Y. Sukirman. Finite element modeling for simulating the sur- face subsidence above a compacting hydrocarbon reservoir. Int. J. Numer. Anal. Methods Geomech., 18:619{639, 1993. [59] K.-A. Lie, S. Krogstad, I. S. Ligaarden, J. R. Natvig, H. M. Nilsen, and B. Ska estad. Matlab Reservoir Simulation Toolbox, 2011. URL http://www. sintef.no/Projectweb/MRST. [60] Y. Lu, D. Elsworth, and L. Wang. Microcrack-based coupled damage and ow modeling of fracturing evolution in permeable brittle rock. Mech Mater, 49:226{ 244, 2013. [61] Y. F. Lu and J. F. Shao. Modeling of anisotropic damage in brittle rocks under compression dominated stresses. Int. J. Numer. Anal. Methods Geomech., 26: 945{961, 2002. doi: 10.1002/nag.230. References 154 [62] M. Mainguy and P. Longuemare. Coupling uid ow and rock mechanics: formu- lations of the partial coupling between reservoir and geomechanics simulators. Oil Gas Sci. Tech., 57:355{367, 2002. [63] A. Mikelic and M. F. Wheeler. Convergence of iterative coupling for coupled ow and geomechanics. Comput. Geosci., 17:455{461, 2013. [64] S. E. Minko, C. M. Stone, S. Bryant, M. Peszynska, and M. F. Wheeler. Coupled uid ow and geomechanical deformation modeling. J. Pet. Sci. Eng., 38:37{56, 2003. [65] M. Negri and C. Ortner. Quasi-static crack propagation by Griths criterion. Math. Models Methods Appl., 18:1895{1925, 2008. [66] M. Nuth and L. Laloui. Eective stress concept in unsaturated soils: Clarication and validation of a unied framework. Int. J. Numer. Anal. Methods Geomech., 32:771{801, 2008. [67] E. Ozkan and R. Raghavan. New solutions for well-test-analysis problems: part 2conputational considerations and applications. SPEJ, 6(3):369{378, 1991. [68] M. Love P. Fitch, S. Davies and T. Pritchard. The petrophysical link between reservoir quality and heterogeneity: Application of the lorenz coecient. In SP- WLA 54th Annual Logging Symposium, number SPWLA, Louisiana, 2013. [69] K. Palaniswamy. Crack Propagation Under General In-Plane Loading. Ph.D. Dissertation, California Institute of Technology, 1971. [70] K. Palaniswamy and W. Knauss. Propagation of a crack under general, in-plane tension. American Society for Testing and Materials, 8 (1):114 { 117, 1972. [71] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin nite element methods for poroelasticity I: the continuous in time case. Comput. Geosci., 11:131{144, 2007. [72] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin nite element methods for poroelasticity II: the discrete-in-time case. Comput. Geosci., 11:145{158, 2007. [73] J. Pogacnik, D. Elsworth, M. OSullivan, and J. OSullivan. A damage mechanics approach to the simulation of hydraulic fracturing/shearing around a geothermal injection well. Comput. Geosci., 29:338{351, 2016. [74] J. Pogacnik, D. Elsworth, M. OSullivan, and J. OSullivan. A damage mechanics approach to the simulation of hydraulic fracturing/shearing around a geothermal injection well. Comput. Geosci., 29:338{351, 2016. References 155 [75] J. Latham Q. Lei and C. Tsang. The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks. Comput. Geotech., 85:151{176, 2017. [76] E. Roelos. Fault stability changes induced beneath a reservoir with cyclic varia- tions in water level. J. Geophys. Res., 93:2107{2124, 1988. [77] M. Schoenberg and C. M. Sayers. Seismic anisotropy of fractured rock. Geophys., 60:204{211, 1995. [78] A. Settari and F. Mourits. A coupled reservoir and geomechanical simulation system. Soc. Pet. Eng. J., 3:219{226, 1998. [79] A. Settari and D. A. Walters. Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction. Soc. Pet. Eng. J., 6:334{342, 2001. [80] J. F. Shao. Poroelastic behaviour of brittle rock materials with anisotropic damage. Mech. Mater., 30:41{53, 1998. [81] J. F. Shao, Y. F. Lu, and D. Lydzba. Damage modeling of saturated rocks in drained and undrained conditions. J. Eng. Mech., 130:733{740, 2004. doi: 10. 1061/ASCE0733-9399. [82] J.F Shao, H. Zhou, and K. Chau. Coupling between anisotropic damage and permeability variation in brittle rocks. Int. J. Numer. Anal. Methods Geomech., 29:1231{1247, 2005. [83] B. B. S. Singhal and R. P. Gupta. Applied hydrogeology of fractured rocks. AHFR, ::13{33, 2010. [84] I. N. Sneddon. The distribution of stress in the neighborhood of a crack in an elastic solid. Proc. Royal Soc. London, 187:229{260, 1946. doi: 10.1098/rspa.1946.0077. [85] O. A. De Swaan. Analytical solutions for determining naturally fractured reservoir properties by well testing. SPEJ, 16(3):117{122, 1976. [86] D.D. Mamora T.D. Bui and W.J. Lee. Transient pressure analysis for partially penetrating wells in naturally fractured reservoirs. In The 2000 SPE rocky moun- tain regional/low permeability reservoirs symposium and exbibition, number SPE 60289-MS, Colorado, 2000. [87] K. Terzaghi. The shearing resistance of saturated soils and the angle between planes of shear,in proc. ICSMFE, 1:6{54, 1936. References 156 [88] L. K. Thomas, L. Y. Chin, R. G. Pierson, and J. E. Sylte. Coupled geomechanics and reservoir simulation. Soc. Pet. Eng. J., 8(4):350{358, 2003. [89] B. Schlager T.P. Frick, C.W. Brand and M.J. Economides. Horizontal well testing of isolated segments. SPEJ, 1(3):261{274, 1996. [90] D. Tran, A. Settari, and L. Nghiem. New iterative coupling between a reservoir simulator and a geomechanics module. Soc. Pet. Eng. J., 9(3):362{369, 2004. [91] D. Tran, L. Nghiem, and L. Buchanan. Improved iterative coupling of geomechan- ics with reservoir simulation. In SPE Reservoir Simulation Symposium, number SPE 93244, Houston, TX, January 2005. [92] J. P. Verdon, J.-M. Kendall, A. L. Stork, R. A. Chadwick, D. J. White, and R. C. Bissell. Comparison of geomechanical deformation induced by megatonne-scale co 2 storage at sleipner, weyburn, and in salah. Proc. Natl. Acad. Sci. U.S.A., 110: E2762{E2771, 2013. [93] J. Wan, L. J. Durlofsky, T .J. R. Hughes, and K. Aziz. Stabilized nite element methods for coupled geomechanics-reservoir ow simulations. In SPE Res. Simul. Sym. (SPE 79694), Houston, Feb. 2003. [94] H. F. Wang. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology. Princeton University Press, 2000. [95] J.E. Warren and P.J. Root. The behavior of naturally fractured reservoirs. SPEJ, 3(3):245{255, 1963. [96] C.C. Wu, S. Freiman, R. Rice, and J. Mecholsky. Microstructural aspects of crack propagation in ceramics. J. Mater Science., 13 (12):2659 { 2670, 1978. [97] W. Chen Y.G Duan. and Q.S. Li. A comprehensive well test model for horizon- tal wells with complex boundaries. In The Annual Technical Meeting, number PETSOC 98-21, Calgary, 1998. [98] Z. Zhang and J. Eckert. Unied tensile fracture criterion. Phys. Rev. Lett., 94:094 { 301, 2005. [99] X. Zhao and B. Jha. Role of well operations and multiphase geomechanics in con- trolling fault stability during CO 2 storage and enhanced oil recovery. J. Geophys. Res., 124, 2019. doi: 10.1029/2019JB017298. [100] J. Zhou, J. Shao, and W Xu. Coupled modeling of damage growth and permeability variation in brittle rocks. Mech. Res. Commun., 33:450{459, 2005.

Abstract (if available)

Abstract Injection of water or gas in a naturally fractured low permeability reservoir often leads to dynamic changes in permeability and elastic moduli of the reservoir due simultaneousagation of multiple pre-existing fractures under injection-induced stress. Wellbore injectivity, reservoir productivity, and spatial distribution of fluids within the reservoir depend on these properties. Numerical modeling of the evolution in the prop- erties require solution of the coupled physical processes of fluid flow, rock deformation, and fracture propagation, which is challenging when multiple fractures are propagating simultaneously due to nonlinearity of the governing equations and presence of multi- ple length scales. We present a numerical framework to model the dynamic evolution in permeability and elastic moduli based on a coupled formulation of multiphase flow, poroelasticity, and fracture-induced damage processes. We solve the fluid flow problem using the finite volume method and the quasi-static deformation problem using the fi- nite element method. We use the continuum damage mechanics method to propagate the fractures and evolve the properties defined at the representative elementary volume scale. We generalize the fixed stress method of solving the coupled flow-deformation problem to model fracturing-induced anisotropy in permeability, elastic moduli and in- duced stresses. The proposed framework can model complex fracture networks with randomly distributed fractures and does not require an explicit representation of frac- tures in the computational mesh. This reduces the computational burden associated with modeling fracture paths and flow through discrete fractures. As a result, fracture- specific multiphase flow functions are not required. We use the proposed framework to study the effects of variable injection/production rate, multiple fracture networks, boundary stress anisotropy, and multiphase flow properties on evolution in the proper- ties. Our work shows that the effect of fracture propagation can be included in existing, sequentially coupled flow-geomechanics models of fractured reservoirs relatively easily by implementing a damage-based solver module within the sequential loop.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

Reactivation of multiple faults in oilfields with injection and production

$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

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

PDF

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

$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

PDF

Modeling and simulation of complex recovery processes

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

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

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

PDF

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

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

PDF

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

PDF

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

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties

PDF

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

PDF

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

PDF

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

Asset Metadata

Creator Bubshait, Ahmed Khaled (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 12/11/2019

Defense Date 12/09/2019

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

Tag coupling,damage mechanics,fracturing,geomechanics,injection,Modeling,multiphase flow,OAI-PMH Harvest,permeability,poromechanics,propagation

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jha, Birendra (committee chair)

Creator Email bubshait@usc.edu

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

Unique identifier UC11673402

Identifier etd-BubshaitAh-8040.pdf (filename),usctheses-c89-249214 (legacy record id)

Legacy Identifier etd-BubshaitAh-8040.pdf

Dmrecord 249214

Document Type Dissertation

Rights Bubshait, Ahmed Khaled

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