UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulation of hydraulic fracture, stress shadow effects and induced seismicity in jointed rock Zangeneh, Neda 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2013_fall_zangeneh_neda.pdf [ 17.49MB ]
Metadata
JSON: 24-1.0074048.json
JSON-LD: 24-1.0074048-ld.json
RDF/XML (Pretty): 24-1.0074048-rdf.xml
RDF/JSON: 24-1.0074048-rdf.json
Turtle: 24-1.0074048-turtle.txt
N-Triples: 24-1.0074048-rdf-ntriples.txt
Original Record: 24-1.0074048-source.json
Full Text
24-1.0074048-fulltext.txt
Citation
24-1.0074048.ris

Full Text

NUMERICAL SIMULATION OF HYDRAULIC FRACTURE, STRESS SHADOW EFFECTS AND INDUCED SEISMICITY IN JOINTED ROCK by Neda Zangeneh   M.Sc., Memorial University of Newfoundland, Canada, 2003  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Geological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  August 2013  ? Neda Zangeneh, 2013   ii Abstract Hydraulic fracturing provides a means to optimize shale gas completions by enhancing the permeability of what is otherwise very tight rock. However, the coupled nature of the processes involved (e.g., thermo-hydro-mechanical-chemical), interlinked with geological variability and uncertainty, makes it extremely difficult to fully predict the spatial and temporal evolution of the hydrofrac and surrounding invaded zone. Numerical design tools have been developed to contend with this complexity, but these have largely focused on the mechanics of brittle fracture propagation at the expense of making simplifying assumptions of the host geology within which the hydraulic fracture is propagating, namely treating it as a linear elastic continuum. In contrast, the reservoir rock conditions are much more complex. Present are natural discontinuities, including bedding planes, joints, shears and faults superimposed by the in-situ stress field. The natural discontinuities under the applied in-situ stress have the potential to either enhance or diminish the effectiveness of the hydraulic fracturing treatment and subsequent hydrocarbon production.  Improved understanding of the interactions between the hydraulic fracture and natural fractures under the stress field would allow designers and operators to achieve more effective hydraulic fracturing stimulation treatments in unconventional reservoirs.  To better account for the presence of natural discontinuities in shale gas reservoirs, this thesis investigates the use of the 2-D commercial distinct-element code UDECTM (Itasca Consulting Group, 1999) to simulate the response of a jointed rock mass subjected to static loading and hydraulic injection. The numerical models are developed to illustrate some important concepts of hydraulic fracturing such as the effect of natural fractures in fracture connectivity, effects of stress shadowing in multiple horizontal well completion, and the   iii effect of fluid injection in induced seismicity, so they can be used to qualitatively evaluate the effects of the in-situ environment on the design and the consequences of the design on the in-situ environment.                 iv Preface This thesis was prepared as a manuscript based thesis.  Chapter 2 ?Investigation of the influence of natural fractures and in-situ stresses on hydrofrac propagation.? was co-authored by Neda Zangeneh, Erik Eberhardt and Marc Bustin. As the lead author, Neda Zangeneh built the numerical models, carried out the analyses, interpreted the results, performed the sensitivity analyses, and prepared the manuscript. The manuscript was reviewed by Erik Eberhardt and Marc Bustin.   Chapter 3 ?Investigation of the influence of stress shadows on multiple hydraulic fractures from adjacent horizontal wells.? was co-authored by Neda Zangeneh, Erik Eberhardt and Marc Bustin. As the lead author, Neda Zangeneh built the numerical models, carried out the analyses, processed the data, performed the sensitivity analyses, and prepared the manuscript. The manuscript was reviewed by Erik Eberhardt and Marc Bustin.  Chapter 4 ?A numerical investigation of fault triggered by hydraulic fracturing.? was co-authored by Neda Zangeneh, Erik Eberhardt, Marc Bustin and Amanda Bustin. The manuscript was reviewed and published in the book Effective and Sustainable Hydraulic Fracturing (edited by A.P. Bunger et al., InTech Publishers, Rijeka, ISBN 978-9535111375, Chapter 23, pp. 477-488). As the lead author, Neda Zangeneh built the numerical models, carried out the analyses, interpreted the results, and prepared the manuscript. The manuscript was reviewed by Erik Eberhardt, Marc Bustin, and Amanda Bustin.      v Table of Contents Abstract .................................................................................................................................... ii Preface ..................................................................................................................................... iv Table of Contents .................................................................................................................... v List of Tables .......................................................................................................................... ix List of Figures .......................................................................................................................... x Acknowledgements .............................................................................................................. xvi Dedication ............................................................................................................................ xvii Inspirations ......................................................................................................................... xviii Chapter  1: Introduction ........................................................................................................ 1 1.1 Problem statement ........................................................................................................ 1 1.2 Thesis objectives .......................................................................................................... 2 1.3 Thesis structure ............................................................................................................ 4 1.4 Literature review .......................................................................................................... 5 1.4.1 Overview ............................................................................................................... 5 1.4.2 Reservoir stimulation ............................................................................................ 6 1.4.3 Hydraulic fracturing .............................................................................................. 7 1.4.4 Modeling of hydraulic fracturing .......................................................................... 9 1.4.4.1 Early modeling .............................................................................................. 10 1.4.4.2 Programs specific to fracture modeling ........................................................ 11 1.4.4.2.1 FRACSIM ............................................................................................. 12 1.4.4.2.2 CVOR ................................................................................................... 12 1.4.4.2.3 WEM ..................................................................................................... 13   vi 1.4.4.2.4 StimPlan ................................................................................................ 13 1.4.4.2.5 FRANC3D/BES .................................................................................... 14 1.4.4.2.6 FEFFLAP .............................................................................................. 14 1.4.4.2.7 MFRAC-II ............................................................................................. 14 1.4.4.2.8 GOHFER ............................................................................................... 15 1.4.4.3 Summary of analysis approaches .................................................................. 15 1.4.4.4 Unified fracture design ................................................................................. 17 1.4.5 Evaluation and selection of numerical approach ................................................ 18 1.4.5.1 Numerical modeling approach ...................................................................... 20 1.4.5.1.1 Distinct element method (DEM) ........................................................... 20 1.4.5.1.2 Bonded particle method ........................................................................ 22 1.4.5.1.3 Finite element/Discrete element (FEM/DEM) ...................................... 24 1.4.5.2 Numerical modeling methodology ............................................................... 25 Chapter  2: Investigation of the Influence of Natural Fractures and In-situ Stresses on Hydrofrac Propagation using a Distinct-Element Approach ........................................... 35 2.1 Introduction ................................................................................................................ 35 2.2 Hydraulic fracture modeling and shale gas reservoirs ............................................... 36 2.3 Numerical modeling methodology ............................................................................ 38 2.4 Influence of natural fractures in hydraulic fractures .................................................. 40 2.5 UDEC simulations ..................................................................................................... 41 2.5.1 Invaded and dilated zones ................................................................................... 41 2.5.2 Influence of natural fractures on initiation of hydraulic fractures ...................... 43   vii 2.5.3 Effect of differential stress and angle of approach on interaction between hydraulic fractures and natural fractures......................................................................... 44 2.5.3.1 Crossing interaction ...................................................................................... 45 2.5.3.2 Offset interaction .......................................................................................... 46 2.5.3.3 Arresting interaction ..................................................................................... 47 2.5.3.4 Bounding relationships as a function of approach angle and stress state ..... 48 2.6 Discussion and conclusions ....................................................................................... 49 Chapter  3: Investigation of the Influence of Stress Shadows on Multiple Hydraulic Fractures from Adjacent Horizontal Wells ........................................................................ 76 3.1 Introduction ................................................................................................................ 76 3.2 Influence of stress shadows arising from multiple hydrofracs .................................. 77 3.3 Numerical modeling methodology ............................................................................ 80 3.4 Model set-up and simulation scenarios ...................................................................... 81 3.4.1 Stress shadowing simulation ............................................................................... 83 3.4.1.1 Conventional hydrofrac (single well) ........................................................... 83 3.4.1.2 Two alternating injection wells (zipperfrac) ................................................. 85 3.4.1.3 Two simultaneous injection wells (simulfrac) .............................................. 86 3.4.2 Scenario comparison ........................................................................................... 86 3.5 Influence of in-situ stress and operational factors ..................................................... 88 3.5.1 Influence of reservoir depth ................................................................................ 89 3.5.2 Influence of in-situ stress ratio ............................................................................ 90 3.5.3 Influence of wellbore spacing ............................................................................. 91 3.5.4 Influence of injection rate ................................................................................... 93   viii 3.6 Discussion and conclusions ....................................................................................... 93 Chapter  4: A numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing............................................................................................................................................... 112 4.1 Introduction .............................................................................................................. 112 4.2 Effect of fluid injection on induced seismicity ........................................................ 113 4.3 Numerical modeling methodology .......................................................................... 115 4.4 Simulation setup and parameters ............................................................................. 116 4.5 Post-injection seismicity simulations ....................................................................... 117 4.6 Fluid pressure build-up during injection and diffusion after shut-in ....................... 117 4.7 Shear slip of critically stressed fault ........................................................................ 118 4.8 Seismic moment and moment magnitude calculations ............................................ 118 4.9 Discussion and conclusions ..................................................................................... 120 Chapter  5: Thesis Discussion and Conclusions ............................................................... 127 5.1 Chapter summaries ................................................................................................... 128 5.1.1 Influence of natural fractures and in-situ stresses on hydrofrac propagation ... 128 5.1.2 Influence of stress shadow on fracture propagation in multi-horizontal wells . 129 5.1.3 Influence of hydraulic fracturing on fault slip and post-injection seismicity ... 129 5.2 Key conclusions and scientific contributions .......................................................... 130 5.3 Further research ....................................................................................................... 132 Bibliography ........................................................................................................................ 134 Appendix .............................................................................................................................. 147     ix List of Tables Table 1.1    Comparison of common geotechnical numerical codes ...................................... 30 Table 2.1    Properties assigned to the modeled planes of weakness and natural fractures .... 52 Table 2.2    Assigned stresses and angle of approach to investigate hydraulic fracture behavior with respect to natural fracture ................................................................................ 52 Table 3.1    Strength and deformation properties assigned to the modeled planes of weakness................................................................................................................................................. 96 Table 4.1    Properties assigned to the modeled planes of weakness and fault ..................... 122                  x List of Figures Figure 1.1   Hydraulic fracturing in a vertical wellbore. The fracture propagation direction is perpendicular to the minimum horizontal stress ..................................................................... 31 Figure 1.2  Conventional (vertical) and unconventional (horizontal) hydraulic fracturing wells ........................................................................................................................................ 32  Figure 1.3  a) PKN fracture model, and b) KGD fracture model .......................................... 33 Figure 1.4  Conceptual model of a rock as used in the Distinct Element Method (DEM) ..... 33 Figure 1.5  Conceptual model of a rock as used in the Bounded Particle Model (BPM) ......  33 Figure 1.6  Discontinuum model: distinct element discretization .......................................... 34 Figure 1.7 Conceptualization (plan view) of the invaded and dilated zones (after Dusseault and McLennan, 2011) ............................................................................................................. 34 Figure 2.1 a) Northeast BC Devonian shale outcrop b) Image log data in a horizontal well through the Maskwa/Otter Park formation, and core log showing: c) partially open fractured) closed fractures in Muskwa/Otter Park formation. After Reine and Dunphy, 2011............... 53 Figure 2.2  Discontinuum model: distinct element discretization using Voronoi tessellation 54 Figure 2.3 Rock mass model with background incipient fractures only, generated using Voronoi tessellation (20m x 20m) .........................................................................................  54 Figure 2.4 Increased fracture aperture showing simulated hydraulic fracture initiation pattern in reservoir rock comprised of incipient fractures only .......................................................... 55 Figure 2.5 Conceptualization (plan view) of the invaded and dilated zones (after Dusseault and McLennan, 2011) ............................................................................................................  55 Figure 2.6   Simulated pore pressure distribution around the borehole .................................. 56   xi Figure 2.7 Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red ............................................................................................................................ 56 Figure 2.8   Displacement vectors around the pressurized borehole ...................................... 57 Figure 2.9 Rock mass model including background incipient fractures superimposed by a persistent fault ......................................................................................................................... 57 Figure 2.10   Displacement vectors around the pressurized borehole and along the fault ..... 58 Figure 2.11 Corresponding shear displacement along the fault in response to the effective stress change ........................................................................................................................... 58 Figure 2.12 Distinct element rock mass model with the distributed natural fractures. Zoomed in window provided above ...................................................................................................... 59 Figure 2.13 Simulated pore pressure distribution showing influence of natural fractures on development of the invaded and dilated zones. Zoomed in window provided above ............ 60 Figure 2.14 Modeled apertures showing interactions that develop between the induced hydrofrac and natural fractures, including the development of invaded and dilated zones. Zoomed in window provided above ....................................................................................... 61 Figure 2.15 Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red .............................................................................................................. 62 Figure 2.16 Diagram of hydraulic fracture intersecting natural fracture ................................ 62 Figure 2.17 Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 60 degrees angle of approach. .............................. 63   xii Figure 2.18 Modeled response of the induced hydrofrac (i.e., incipient fractures that have opened) as it intersects a natural fracture at a 60 degree approach angle under high differential stresses. In this case the hydrofrac crosses the natural fracture. Zoomed in window provided above .......................................................................................................... 64 Figure 2.19  Simulated pore pressure distribution for scenario assuming a natural fracture at a 60 degree approach angle under high differential stresses. Zoomed in window provided above ....................................................................................................................................... 65 Figure 2.20 Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above ................................................... 66 Figure 2.21 Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 30 degree angle of approach ................................. 67 Figure 2.22  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 30 degrees approach angle under low differential stresses. In this case the hydrofrac path is offset by the natural fracture. Zoomed in window provided above ........................................ 68 Figure 2.23 Simulated pore pressure distribution for scenario assuming a natural fracture at a 30 degree approach angle under low differential stresses. Zoomed in window provided above. ...................................................................................................................................... 69 Figure 2.24 Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above ................................................... 70 Figure 2.25  Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 45 degree angle of approach ................................. 71   xiii Figure 2.26 Modeled response of the induced hydrofrac as it intersects a natural fracture at a 45 degree approach angle under intermediate differential stresses. In this case the hydrofrac is stopped by the natural fracture. Zoomed in window provided above ................................. 72 Figure 2.27 Simulated pore pressure distribution for scenario assuming a natural fracture at a 45 degrees approach angle under intermediate differential stresses. Zoomed in window provided above. ....................................................................................................................... 73 Figure 2.28 Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above ................................................... 74 Figure 2.29   Correlation between the angle of approach and differential stress based on UDEC based numerical experiments ...................................................................................... 75 Figure 3.1  Discontinuum model: distinct element discretization .......................................... 97 Figure 3.2 Rock mass model assuming two orthogonal sets of incipient fracture planes (i.e., planes of weakness), parallel and perpendicular to bedding. Shown are the locations of the two horizontal wellbores and 8 monitoring points referred to in subsequent figures ............. 97 Figure 3.3   Three dimensional configuration of the two horizontal wellbores ...................... 98 Figure 3.4 Conventional hydrofracing scenario (single well): Fluid pressure histories for different monitoring points between the two wellbores, with injection occurring from Wellbore A .............................................................................................................................. 99 Figure 3.5 Conventional hydrofracing scenario (single well): Stress histories for different monitoring points between the two wellbores, with injection occurring from Wellbore A . 100 Figure 3.6  Stress distribution, resulting from the pressurization of wellbore A for a period of 90 minutes ............................................................................................................................. 100   xiv Figure 3.7 Zipperfrac scenario between two alternating injection boreholes: Fluid pressure histories for different points between the two wellbores. At 90 minutes, the injection at wellbore B is shut-in and injection started at wellbore A ..................................................... 101 Figure 3.8  Zipperfrac scenario between two alternating injection boreholes: Stress histories for different points between the two wellbores. Point 1 is closest to the initial injection (wellbore B), with point 8 being closest to the subsequent injection (wellbore A) .............. 102 Figure 3.9 Stress distributions arising from the zipperfrac scenario, after injection into: (a) wellbore B for a period of 90 minutes, and (b) then into wellbore A, for a period of 90 minutes .................................................................................................................................. 103 Figure 3.10 Simulfrac scenario between two simultaneously injecting hydrofrac wellbores: Fluid pressure histories for different points between the two wellbores. Point 1 is closest to injection well B and point 8 closest to injection well A ....................................................... 104 Figure 3.11 Simulfrac scenario between two simultaneously injecting hydrofrac wellbores: Stress histories for different points between the two wellbores. Point 1 is closest to injection well B and point 8 closest to injection well A ...................................................................... 105 Figure 3.12 Stress distributions around simultaneously pressurized wellbores during pressurization (simulfrac scenario): (a) after 20 minutes of injection, (b) after 40 minutes of injection, (c) after 60 minutes of injection, (d) after 90 minutes of injection ....................... 106 Figure 3.13   Hydraulic fracture comparison between the three modeled scenarios ............ 107 Figure 3.14 Maximum induced hydraulic fracture aperture comparison between the three modeled scenarios ................................................................................................................. 107 Figure 3.15   Hydraulic fracture length comparison between the three scenarios ................ 108   xv Figure 3.16 Normalized hydraulic fracture length comparison between the three scenarios for different reservoir depths ...................................................................................................... 108 Figure 3.17 Normalized hydraulic fracture length comparisons between different completion scenarios for different horizontal to vertical in-situ stress ratios .......................................... 109 Figure 3.18  Geometry of an elliptical opening in a biaxial stress field ............................... 109 Figure 3.19  Hydraulic fracture comparison for different wellbore spacing ........................ 110 Figure 3.20 Normalized hydraulic fracture length for simulfrac and zipperfrac scenarios with different wellbore spacing..................................................................................................... 110 Figure 3.21 Normalized hydraulic fracture length for simulfrac and zipperfrac scenarios with different injection rate ........................................................................................................... 111 Figure 4.1  Rock mass model with two sets of weakness planes and a fault ........................ 123 Figure 4.2  Pore pressure distribution at the time of shut-in ................................................. 123 Figure 4.3  Pore pressure distribution 2 ? hours after shut-in .............................................. 124 Figure 4.4  Fault shear displacement 2 ? hours after shut-in ............................................... 124 Figure 4.5 Saturation of various magnitude scales: Mw (moment magnitude), ML (Richter local magnitude), Ms (surface wave magnitude), mb (short-period body wave magnitude), mB (long-period body wave magnitude), and MJMA (Japanese Metrological Agency magnitude).After Idriss (1985) and Kramer (1996) .............................................................  125 Figure 4.6 Relationship among various scaling parameters for earthquakes. The larger the earthquake, the larger the fault and amount of slip, depending on the stress drop in a particular earthquake. Observational data indicate that earthquake stress drops range between 0.1 and 10 MPa (modified after Zoback and Gorelick, 2012) .............................................  126    xvi   Acknowledgements  I express my gratitude to the faculty, staff and my fellow students at the University of  British Columbia (UBC), who encouraged me to complete my work in this research. I owe particular thanks to Professors Erik Eberhardt and Marc Bustin, whose critical questions taught me to question any aspect of this research more deeply. I thank Professor Roger Beckie for providing coherent answers to my questions, and supporting me in different aspects of this thesis. I also thank Professor Oldrich Hungr for his encouraging and supporting statements throughout this study. The financial support of Natural Sciences and Engineering Research Council of Canada (NSERC), Trican Well Service, Faculty of Graduate Studies (FoGS), and Geoscience BC facilitated the completion of this research. Heartfelt gratitude and appreciation owed to my parents Fereydoun and Mehri, who have supported and encouraged me in all my pursuits. My sisters Noushin and Negar have always been a source of inspiration, encouragement and persistence. My brothers, Hamed and Hesam have been a reason for achievement as being their role model. My parents-in-laws Reza and Rouhi have been another source of encouragement throughout this study.  I would like to thank my loving, supportive, encouraging, and patient husband Dr. Ali Azizian, whose faithful support made this Ph.D. research possible. Also I would like to thank the love of my life, my joyful daughter Dana, for holding on to a dream to have a mother who no longer has to study days and nights, weekdays and weekends.   xvii Dedication    I dedicate this dissertation to my beloved Father,  Fereydoun Zangeneh who taught me: how to live, how to love, and how to learn ?                  xviii Inspirations  ?Don't be satisfied with stories, how things have gone with others. Unfold your own myth.? Molana Jalal eddin Rumi  ?No one saves us but ourselves. No one can and no one may. We ourselves must walk the path.? Gautama Buddha  ?A man is but the product of his thoughts. What he thinks, he becomes.? Mahatma Gandhi  ?I am thankful to all those who said no. It's because of them I did it myself.? Albert Einstein  "Genius is one percent inspiration, ninety-nine percent perspiration." Thomas Edison  ?Everything should be made as simple as possible, but not simpler.? Albert Einstein  ?Science may set limits to knowledge, but should not set limits to imagination.? Sir Bertrand Russell   ?To do this job, I have worked infinitely hard; whoever works just as hard would reach at the same result.? Johann Sebastian Bach   1 Chapter  1: Introduction 1.1 Problem statement Hydraulic fracturing is a process for maximizing the extraction of oil and natural gas resources by fracturing the reservoir rock by applying hydraulic pressure to enhance its permeability. Designing a treatment to achieve the desired results is intimately connected with different factors such as geology, in-situ stress, pumping pressure and fluid flow. The coupled nature of the processes involved in hydraulic fracturing (e.g., thermo-hydro-mechanical-chemical), interlinked with geological variability and uncertainty, makes it extremely difficult to firmly establish the different variables required for designing a hydrofrac. To better understand the behavior of rock during a hydraulic fracturing treatment, numerous studies have been undertaken and several software programs developed to predict reservoir rock behavior during the hydraulic injection. Most of these are designed in a way that treats the rock mass as an elastic continuum. A different technique that treats the geology in a more realistic way, is needed to account for not only the natural fractures, but also the fractures and geological features that are not clearly visible. With the oil and gas industry now moving more towards fracturing shale with ultra-low permeability, i.e. unconventional reservoirs, the success of hydraulic fracturing stimulation will depend on developing a better understanding of hydraulic fracture growth in such unconventional reservoirs. The multiple hydraulic fractures generated from neighboring horizontal wells can provide a better completion option for stimulating unconventional and low permeability oil and gas reservoirs (e.g, Yost and Overbey, 1989). While the industry has   2 made significant advancement in developing mechanical systems for carrying out these completions, progress in the design of multiple hydrofracs and overall production optimization requires a more thorough mechanistic understanding of their influence.  It has also been well established that the fluid injections during hydraulic fracturing also induce seismicity (e.g. Hollister and Weimer, 1968). The occurrence of induced seismicity represents an important consideration in the risk profile of an injection well. Estimation of the corresponding maximum seismic magnitude that is likely to occur in response to the hydraulic fracture treatment is very important. Improved understanding of this condition will allow designers and operators to control the amount of injection during a hydraulic fracturing treatment to minimize fluid pressure diffusion and subsequent slip of a neighboring fault. 1.2 Thesis objectives The objective of this research is to further the understanding of hydraulic fracturing by addressing limitations in current empirical and continuum-based numerical design by explicitly accounting for the controlling influence of natural discontinuities present in the rock on hydraulic fracture propagation and stimulation response. Thus, the research aims to examine how reservoir geology impacts the design of hydraulic fracturing stimulations and influences production efficiency of natural gas reservoirs. To do so, several thesis objectives are defined as listed below.  i. To approach the modeling and analysis of hydraulic fracturing from a discontinuum framework, in which the controlling influence of geological structures is explicitly accounted for in the model. Rock masses are comprised of both an   3 intact rock component and the discontinuities that transect it. The discontinuities include faults, joints and bedding planes. The presence of the discontinuities generally controls the behavior of the rock mass under static or dynamic loading (Goodman, 1976).  ii. To investigate the influence of natural fractures on the modeling of hydraulic fracture propagation. The analysis is carried out for several geometrical variants, related to hypothetical geological scenarios simulating a naturally fractured reservoir, and their effect on hydraulic fracture size, orientation and path. Also examined are the propagation limiting effects of leak-off to the natural fractures and enhancing effects of reactivating natural fractures instead of using energy to fracture the intact rock.  iii. To investigate the effect of different multiple horizontal well completion scenarios on the propagation of hydraulic fractures with respect to altered in-situ stress fields during stimulations. The effects of multiple hydraulic fractures on stress shadowing are also studied as a function of the reservoir depth, in-situ stress ratio, wellbore spacing, and injection rate.  iv. To investigate the relationship between hydraulic fracturing and the response of a naturally fractured rock mass to transient fluid pressures. The relationship between injection pressure, fluid diffusion and maximum magnitude of an induced seismic event is analyzed. The simulation is then used to interpret the associated seismic events and their relationship to the injections and shut-in pressures, and to estimate the maximum magnitude of the induced seismic event.    4 1.3 Thesis structure This thesis consists of an introduction chapter, three research chapters, and a conclusion. Chapter 1 introduces the problem and explains the thesis objectives and the literature review carried out. The next three chapters elaborate on the procedures and findings of the research conducted to achieve the stated objectives, each using 2-D distinct-element modeling techniques. Chapter 2 evaluates the use of 2-D distinct element techniques, including Voronoi tessellation procedures, to account for both shear and dilation along existing joints and tensile rupture of intact rock as a result of effective stress change following the fluid injection. The analysis is carried out for several geometrical variants, related to hypothetical geological scenarios simulating a naturally fractured reservoir, and their effect on hydraulic fracture size, orientation and path. Also the propagation limiting effects of leak-off to the natural fractures and enhancing effects of reactivating natural fractures instead of using energy to fracture the intact rock are examined. Chapter 3 describes a series of numerical experiments investigating the influence of stress shadow on multiple hydraulic fracture treatments. The utility of the 2-D distinct element approach is evaluated for modeling shear slip and aperture dilation along existing rock joints and tensile rupture of intact rock as a result of changes to the effective stress field following fluid injection. The analysis is carried out for different completion scenarios to investigate their effect on the propagation of hydraulic fractures. The changes in fluid pressure and corresponding effective stress changes around each wellbore during different completion techniques are examined. The effects of multiple hydraulic fractures on stress   5 shadowing are also studied as a function of the reservoir depth, in-situ stress ratio, wellbore spacing, and injection rate. Chapter 4 investigates the relationship between hydraulic fracturing (i.e. fluid injection) and the response of a naturally fractured rock mass, including fault slip, to transient fluid pressures through a series of 2-D distinct element simulations. The conceptual reservoir model consists of a critically stressed fault plane and the surrounding rock mass containing planes of weakness, for which a hydraulic fracture is numerically simulated and the response modeled using a transient, coupled hydro-mechanical solution. The results demonstrate the influence of fluid diffusion generated by the fracturing fluid after shut-in on the triggering of fault slip. The simulation is then used to interpret the associated seismic events and their relationship to the injections and shut-in pressures, and to estimate the maximum magnitude of the induced seismic event. 1.4 Literature review 1.4.1 Overview Hydraulic fracturing (or hydrofrac) is a method to improve well productivity in tight reservoir rocks. Hydraulic fracturing is a procedure for maximizing the extraction of underground resources, such as oil and natural gas, by fracturing the resource rock using hydraulic pressure. Designing a treatment to reach the best possible results is involved with: i) the reservoir geological conditions and the properties of rock mass which control the geometry of fracture; ii) the orientation and magnitudes of the principal stresses which control the direction and length of hydrofrac propagation, iii) the mechanics of fluid which controls flow of fluid in the fractures and the placement of proppant in the fractures; and iv)   6 the frac fluid chemistry which governs the functioning of the treatment fluid used to perform the fracturing treatment (Economides and Nolte, 2000). The coupled nature of the processes involved in hydraulic fracturing (e.g., thermo-hydro-mechanical-chemical), interlinked with geological variability and uncertainty, makes it extremely difficult to firmly establish the different variables required for designing a hydrofrac. In this thesis numerical models are elaborated to illustrate some important concepts of hydraulic fracturing, so they can be used to qualitatively evaluate the effects of the in-situ environment on the design and the effects of the design on the in-situ environment. It is desirable to study several introductory papers that were published between the late 1950s and early 1970s, which establish the foundation of hydraulic fracture modeling, as well as the different hydraulic fracture modeling software currently available to evaluate their strengths, weaknesses and ability to simulate hydraulic fracture propagation. 1.4.2 Reservoir stimulation The two main activities in the petroleum industry to extract oil and natural gas from a reservoir are artificial lift and reservoir stimulation (Economides and Nolte, 2000). Artificial lift is a method to increase the flow of oil or natural gas from a production well.  The flow of oil or natural gas is increased by using mechanical devices such as pump or velocity string inside the well or by injecting gas into the liquid down the production well to decrease the weight of the hydrostatic column. The produced fluid typically involves a mixture of oil, water and/or gas (Economides and Nolte, 2000). The unconventional reservoirs have a rock matrix with very small inter-particle pore networks having very poor fluid-flow characteristics (generally, less than 1 millidarcy   7 permeability). A great amount of oil or gas can be trapped in these rocks. To increase the permeability of these reservoirs, stimulated fracture networks should be added to the natural fractures so the reservoir achieves an economic recovery of the hydrocarbon. Rock types in this category include shales and methane-bearing coal. The term shale refers to any rock consisting of a fine graines with extremely small pores storing hydrocarbon. The unconventional reservoirs in some regions like the Rocky Mountains are the desired places for reservoir stimulation activity (Cramer, 2008).  The reservoir stimulation is performed to increase the well productivity by increasing permeability around the wellbore. The stimulation can be done by injecting acids or by creating fractures to increase the conductivity of the formation. The first of these are referred to as ?matrix stimulation? and the latter as ?hydraulic fracturing? (Economides and Nolte, 2000). These two ways can be used to improve the natural connection between the wellbore and the reservoir. 1.4.3 Hydraulic fracturing Hydraulic fracturing is done by injecting fluids into a reservoir formation. As the flow rate continues, the fluid pressure increases. When the pressure of the fluid exceeds the formation rock maximum stress, it causes the formation rock to fracture. (Economides and Martin, 2007). One of the early applications of hydraulic fracturing was for estimating in-situ stresses in rock in which the localized rock mass responses to injection were modeled to resolve the in-situ stress components. In those models, the geological environment was   8 simplified by considering the rock mass as an equivalent continuum and therefore applying continuum mechanics. Hubbert and Willis (1957), Ito and Hayashi (1993), and Haimson and Cornet (2003) are some examples of the main studies done to determine in-situ stress orientation and magnitude using hydraulic fracturing.  Hydraulic fracturing is the process of creating a fracture initiated by applying hydraulic loading using a pressurized fluid (Adachi et al. 2008). When the pumping fluid does not escape into the surrounding rock as fast as it is injected into the formation, pressure rises and at some point it exceeds the rock strength and as a result a fracture initiates. The initiated fracture propagates in the direction of the major principal stress (i.e. it opens in the minor principal stress direction in which the work needed for cracking the rock is minimum). In the early days, most wells were vertical and because the reservoirs were at great depths, the maximum stress was vertical gravity loading. As a result, the initial fracture initiated in the vertical direction and was planar and perpendicular to the minimum horizontal stress (Figure 1.1). As the fracture grows, new formation area becomes exposed to the injected fluid and the rate of fluid leaking off into the formation gradually increases. However, if the pumping rate is sustained in a way that it becomes higher than the fluid-loss rate, then the hydraulic fracture will continue to propagate. Once pumping stops, the fracture gradually closes as the high-pressure fluid in the fracture leaks off, and the newly-created surface area in the formation will be less effective for production. To keeping the fractures open, a propping agent like sand or a high-strength granular substitute for sand is usually added to the hydraulic fracturing fluid.   9 The 3-D coupled hydromechanical processes involved are very complex and hence, the description of the hydraulic fracturing process is highly idealized. Therefore it is important to be able to increase the effectiveness of the hydraulic fracture treatment and generate more productive wellbores. Hydraulic fracturing can be used to: i) bypass wellbore damage to return a well to its initial condition for having a better productivity; ii) make the formation more conductive by fracturing and increase its productivity; or iii) modify fluid flow in the formation. These three goals can overlap and affect the neighboring wells with respect to placement of other wells and the number of wells to be drilled. The increase in natural gas production following stimulation is an important objective of a successful hydraulic fracture design. To achieve higher production, hydraulic fracturing techniques started from vertical production wells.  Treatments then evolved to horizontal wells for providing greater well contact with a larger volume of reservoir rock than vertical wells. Evolving towards horizontal fracturing well over vertical well was due to the limited vertical thickness of the reservoir compared to its lateral extent (Figure 1.2). 1.4.4 Modeling of hydraulic fracturing The main purpose of reservoir stimulation is to enhance the reservoir permeability, assisting the more efficient hydrocarbons delivery, and also to increase ultimate economic recovery.  An evaluation of the economics of a hydraulic fracturing treatment includes: increase in production, treatment cost, and risks associated with mechanical problems. For these tasks, computer models can be used to model the flow of fluids through the porous media.    10 1.4.4.1 Early modeling Several introductory papers were published between the late 1950s and early 1970s that developed the foundation for hydraulic fracture modeling by introducing two of the most common models used in early fracture treatment design, the PKN (Perkins-Kern-Nordgren) and KGD (Khristianovic-Geertsma-de Klerk) models.  The difference between the two models is on how they convert a 3-D fracture mechanics problem into a 2-D plane strain problem.  Figure 1.3 illustrates these models. For a fracture treatment design, the height of the fracture is considered to be fixed or completely confined in a given layer. One of the two assumptions is made in each model with regards to the ratio between the length and height of the fracture. The PKN geometry (Figure 1.3a) is used when the fracture height is much smaller than the fracture length (Perkins and Kern, 1961) and the KGD geometry (Figure 1.3b) is used where the fracture height is much greater than the fracture length. In this case, the fracture width in the vertical direction changes more slowly than in the horizontal direction (Geertsma and de Klerk, 1969). Therefore all horizontal cross sections act equivalent and they are identical. Thus, the PKN model has a constant height with an elliptical vertical cross-section. The KGD model has a constant height with constant width along the height. In both models, the rock is assumed to be a continuous, homogeneous, isotropic linear elastic solid. If the parameters, such as in-situ stress, Young?s modulus, formation permeability and total leak-off coefficient are entered correctly, the model will give realistic estimates of generated fracture length and width subject to using correct value of fracture   11 height in the 2-D model.   Fluid flow in these models is assumed to occur only along the length of the fracture and the width and shape of the fracture are related to the fluid pressure distribution in the fracture.  These 2-D models offer a simplified but computationally efficient means to model hydraulic fracture propagation. However, the 2-D models are not able to simulate both lateral and vertical propagation. With advances in computation capabilities, an evolution was made towards pseudo-3-D models (P3D). In P3D models, the assumption of constant and uniform height was not considered (Settari and Cleary, 1986; Morales, 1989). The height in P3D models is a function of position along the fracture and time. In P3D models, the fracture length is assumed much greater than the height. Adding the vertical fluid flow element is the main difference between the P3D and the 2-D PKN and KGD models.  P3D models have been favored by many fracture design engineers (Economides and Nolte, 2000). P3D models are better than 2-D models for most situations. They compute the height, width and length of the fracture using the data for the pay zone. P3D models have also some limitations. For example they cannot handle fractures of arbitrary shape and orientation and they are not able to do fracture network simulation. In that case, full 3-D models are required.  1.4.4.2 Programs specific to fracture modeling For a given reservoir, a set of values for various treatment parameters like fluid injection rate and time, viscosity of fracturing fluid, and proppant concentration must be chosen for creating a hydraulic fracture geometry that meets the design objective (Rahman et al., 2001). Therefore, improvements are made in field equipment and technology as well as software for   12 modeling, monitoring and interpreting hydrofrac generation. Industry uses several of the codes widely to design and optimize fracture treatments. These codes are described in the following sections. 1.4.4.2.1 FRACSIM FRACSIM (Ventura, 1985) is a single-phase, 2-D reservoir model. It is used to simulate flow of gas through a porous medium for vertically fractured gas wells. The rock mass is treated as continuum and the required input parameters are the reservoir pressure, formation permeability, and fracture half-length. The program can be used to investigate various single-well problems such as non-Darcy flow, damage around the fracture, wellbore storage, and fracture permeability reduction over time. It can also be used to model different fracture lengths for a specific well and subsequently use the output results to FRACOP (Ventura, 1985) to estimate the associated cost. The estimation of cost allows the comparison of costs for each modeled fracture length to select an optimum fracture length based on reservoir permeability.  1.4.4.2.2 CVOR Consortium for Virtual Operations Research, CVOR (Mohaghegh et al., 1999) is a hybrid software toolbox. The model is used to design and optimization of hydraulic fracture treatment consisting of neural network and genetic algorithm routines in Virtual Fracturing Optimizer program. In cases where the detailed reservoir data are limited, unavailable or expensive, the use of CVOR is favorable. The virtual intelligence techniques can be applied using only the basic well information, production history, and results of previous hydraulic fracturing job treatments. Alternatively, in cases where no hydraulic fracture treatments have   13 yet been performed, reservoir data such as permeability profile data can be used instead of production history and historical data from prior treatments. The rock mass in this method is treated as an equivalent continuum. This method uses neural networks that establish relationships between the reservoir conditions and treatment parameters instead of mathematical formulations. The model computes the fracture geometry by processing historical treatment and well productivity data. Accordingly, this methodology is not a substitution for physics-based approaches (Rahman et al., 2003). 1.4.4.2.3 WEM Well Evaluation Model (Reynolds et al., 2004) is a program for simulating the fractured wellbore performance by means of a fracture size optimization plot after a successfully placed fracture treatment. This plot shows the predicted flow rate versus fracture half-length for different formation permeabilities. It also models the well performance versus fracture conductivity. The rock is modeled here as an elastic continuum surrounding the propagating hydraulic fracture. 1.4.4.2.4 StimPlan StimPlan (Reynolds et al., 2004) is a pseudo-3-D numerical modeling code that treats the rock mass as a continuum and utilizes an implicit finite difference formulation to determine the optimum fracture dimensions, sizing, and sand schedule for hydro-fracture design. The program provides conductivity contour plots along the fracture length.     14 1.4.4.2.5 FRANC3D/BES The FRANC3D/BES software system, which stands for FRacture ANalysis Code for 3-D problems using a Boundary Element Solver, is a program developed to simulate crack initiation, reorientation and coalescence (Dong and de Pater, 2002). The software treats rock mass as a continuum and it can model multiple crack propagation in 3-D. It can also model the slow reorientation of a crack to the preferred crack plane as a result of a decrease of the differential stress around the crack or an increase of pressure within the crack. 1.4.4.2.6 FEFFLAP The Finite Element Fracture and Flow Analysis Program, FEFFLAP (Shaffer et al., 1984), is another numerical modeling program that models fluid-driven discrete fracture propagation in jointed media. It was developed in support of the Unconventional Gas Recovery Program at the Lawrence Livermore National Laboratory. FEFFLAP is a 2-D finite element program. It combines linear solid elements with nonlinear joint elements, coupling discrete fracture propagation with fluid flow in the fracture. 1.4.4.2.7 MFRAC-II This hydraulic fracturing program is a pseudo-3D simulator that accounts for the coupled hydro-mechanical parameters that affect the propagation of fracture and proppant transport. The fracture propagation in the rock is simulated by coupled hydro-mechanical equations. Integral methods are utilized to solve the partial differential equations (Meyer, 1989). The program can treat the rock as a continuum and recently added discontinuity by means of adding Discrete Fracture Network (DFN) (Meyer, 2011).   15 1.4.4.2.8 GOHFER GOHFER is a finite difference code for calculating fluid flow within a crack using an integral equation for fracture width. The program considers the rock as a continuum and is able to consider any number of randomly arranged rock formations. The model uses the linear-elastic equation for the displacement of a semi-infinite half space under a distributed load. In this model, pressure distribution in the fracture is calculated using a 2-D finite difference formulation. A coupled solution to fluid pressure and fracture width are utilized in the program (Barree, 1983). 1.4.4.3 Summary of analysis approaches The main draw-backs of the hydraulic fracturing design and optimization software reviewed can be summarized as follows (e.g. Rahman et al., 2003): i. In most models, parametric sensitivity analysis is used to reach the optimum design, which is a lengthy process and cannot always achieve the globally optimum design. ii. In most models, the net present value (NPV) is used to evaluate achieving an optimum design. There are other potential design objectives in addition to NPV, such as maximum production, and minimum treatment costs which need to be considered in optimality of a design. iii. In most models, for finding the optimum design, a few parameters are varied while most of the parameters kept fixed. In most cases, a non-optimum fracture geometry and reduced reservoir productivity become a result of such an optimum design.   16 iv. In most models, mathematics and physics of design models are not presented comprehensively. Most of the models are software/simulators with partial flexibility with respect to considering other important factors in reservoir stimulation. v. And most importantly, all of these methods treat the rock mass as an elastic continuum. A different technique, that treats the geology in a more realistic way, is needed to account for not only the natural fractures, but the fractures and geological features that are not clearly visible. The Net Present Value (NPV) is usually used for hydraulic fracturing treatment design by the industry (Meng and Brown, 1987; Balen et al., 1988; Economides et al., 1994; Aggour and Economides, 1998). Improvement in NPV is done by parametric sensitivity analyses. The height of the fracture is known (often equal to the pay zone height) and for each of the treatment parameters a set of discrete values are assumed. The NPV is then calculated for different treatment parameters combination and a known fracture length, (Balen et al., 1988). This procedure is repeated for different fracture lengths. Then the NPV versus fracture length curves are produced. Using these curves, the values of stimulation parameters and the fracture length are optimized based on a maximum value of NPV. This procedure can be used to optimize the hydrofrac design but it is not possible to examine all potential design scenarios and to satisfy the operational limitations and fracture growth control requirements. Also, there are other design objectives that are needed to be considered for an optimum hydraulic fracturing design. Those include maximizing production and minimizing treatment cost. Thus, the optimization of hydraulic fracturing   17 design can be done by finding a ?best possible? set of values for treatment and other parameters while maximizing the revenue/production. 1.4.4.4 Unified fracture design Economides and Demarchos (2008) state that almost all commercial ?hydraulic fracture models? are simulators that use injection variables (fluids and proppant) to predict fracture geometry, whereas ?fracture design models? first determine the optimum fracture geometry to maximize production and then establish the best way to achieve the desired geometry. Most industry models do not provide for optimum fracture performance suitable for production engineering applications. The best proppant schedule is selected by trial and error which means ?maximum conductivity? is difficult to achieve. Economides and Demarchos (2008) proposed the Unified Fracture Design (UFD) theory to ?pinpoint the optimum fracture dimensions? and then use P3D/3D model to design the necessary injection. This approach accounts for constraints, real reservoir configurations, and heterogeneity. The word ?unified? has been selected to mean that the design process is applicable to all petroleum reservoirs: low permeability to high permeability, hard rock to soft rock.  In this approach, treatment sizes are unified as they are characterized by the dimensionless Proppant Number. The proppant number determines the theoretically optimum fracture dimensions for achieving the maximum productivity index.  For a chosen proppant mass, the fracture size is physically optimized and based on the selected fracture dimensions, the revenue is maximized for certain economic criteria (e.g. NPV). In this method, the mass of the proppant used in a treatment can be used as the main design   18 optimization variable for hydraulic fracturing treatment. The NPV can be calculated for different proppant masses. The NPV will be maximized for one value of the proppant mass, within range of proppant mass determined by the UFD model (Marongiu-Porcu et al., 2008).  After using a model to predict fracture propagation and performance, in is necessary to perform a trial and error process to determine the optimum treatment design from a set of calculated physical designs for achieving an economic optimization. 1.4.5 Evaluation and selection of numerical approach Numerical modeling includes a range of methods which are different based on the way they represent the problem domain and solve the governing equations. Key methods include boundary element, finite element and finite difference. These methods can be grouped as methods that treat the problem domain (i.e. rock mass) as a continuum while the distinct and discrete element methods treat the problem domain as a discontinuum (Hoek et al., 1991). Each of these methods is capable of solving complex hydrogeological and geomechanical interactions and has been widely applied to a range of civil and mining geomechanics problems, but has not been fully employed for hydraulic fracture modeling. A detailed evaluation was undertaken to explore which methods would be most suitable for the research objectives defined. A number of studies was reviewed in which a fracture mechanics approach was used to model tensile fracture growth under different conditions such as stress situation, reservoir and fluid properties and pumping rates (e.g., Detournay and Cheng, 1991). Boundary element and finite element numerical models have been used by Shah et al. (1997), Warpinski et al. (1994) and Vandamme and Curran (1989) to simulate 3-D hydrofractures. One of the important findings of these studies is that the   19 traditional analytical solution and numerical models do not correctly predict the hydraulic fracture determined by microseismic records. These models only assume tensile fracture growth whereas a considerable amount of shear-type seismic events is usually observed in microseismic records (Al-Busaidi et al., 2005). Therefore, a number of researchers have recently been investigating the use of discontinuum methods like the distinct element method. With these techniques, the continuum is divided into distinct blocks or particles. The fluid can flow between those distinct blocks. This method allows a better understanding of hydraulic fracture growth in rocks which may be occupied with sets of pre-existing fractures. These studies show that the degree of refinement of the numerical representation of the flow network depends on the adopted mechanical discretization; the pattern and size of the discontinuities and assemblage of blocks in a hydraulic fracture simulation (Choi, 2003).  Table 1.1 summarizes a comparison between several commercially available state-of-the-art numerical modeling codes based on the different methods described above. Also included are their respective advantages and disadvantages with respect to the research objectives. Based on this comparison, the boundary element, finite element and finite difference codes are not suitable as the continuum representation of the problem domain they require explicitly neglects the influence of discontinuities. Only the discontinuum codes are capable of simulating fracture propagation in a jointed rock mass and are chosen here as the primary methods to be used. These are described in more detail in the following section.     20 1.4.5.1 Numerical modeling approach Numerical methods, through a variety of constitutive material models, can be used to qualitatively and quantitatively evaluate the effects of geology on a hydrofrac design and the consequences of the hydrofrac on modifying the geology. To meet the research objectives defined, it is desirable to be able to include in the model a large number of random or oriented discontinuities to more accurately represent the planes of weakness that exist within tight reservoir rock masses. It is hypothesized that these, combined with the superimposed stress field, can have a significant effect on the propagation direction and characteristics of a hydrofrac. For this discontinuum approaches are ideally suited. These can be divided into two classes: distinct element and bonded particle methods.  The advantages of these methods over continuum-based methods like finite element are listed by Cundall and Hart (1992) as:  ? they allow finite displacements as well as rotations of discrete bodies, including slip and detachment; ? they recognize new contacts between bodies automatically during calculations.  Without the first characteristic, some important mechanisms in a discontinuous medium cannot be molded by the program; without the second, a limited numbers of bodies, which their interactions are known in advance, can only be modeled. 1.4.5.1.1 Distinct element method (DEM) The Distinct Element Method (DEM) is a Lagrangian numerical technique in which the problem domain is divided through by discontinuities of variable orientation, spacing and   21 persistence, which may undergo large deformations in shear and opening in response to a change to the modeled stresses. The blocks, although deformable, remain intact. A detailed discussion of Distinct Element Method is presented by Cundall and Hart (1992). Figure 1.4 shows an idealized DEM discretization of a rock mass with two orthogonal sets of discontinuities. The DEM also includes coupling to computational fluid dynamics formulations. One fundamental advantage of the DEM is that the rock mass pre-existing joints can be integrated into the model and allowed to undergo large deformations (Bobet et al., 2009). The commercial code developed based on this method is UDECTM (Universal Distinct Element Code; Itasca Consulting Group, 1999), a 2-D code that simulates the response of jointed rock masses under either static or dynamic loading. Additional control and customization are available in UDEC through a powerful built-in programming language FISH. UDEC is capable of modeling the behavior of jointed rock masses. The important controlling factors include yielding of weak rock in tension and slip along pre-existing discontinuities. Progressive failure associated with crack propagation can be modeled by breaking of pre-existing contacts between pre-defined blocks. The blocks are deformable but they remain intact during the crack propagation simulation.  Key for simulating hydraulic fracturing, UDEC has the capability to model fluid flow through the fractures. A fully coupled hydro-mechanical analysis can be performed, in which fracture conductivity is dependent on mechanical deformation and, conversely, joint water   22 pressures affect the mechanical computations. The blocks in this assemblage are treated as being impermeable. 3DEC is the three dimensional version of UDEC program which is capable of modeling of hydraulic fracturing in three dimensional space. In one such recent example, McLennan et al. (2010) modeled and assessed complex fracture growth and predicted production associated with the fracture propagation through a generated fracture using 3DEC. A naturally fractured reservoir was built to model complex, non-planar, hydraulic fracture propagation and the results showed the complex fracture network evolution during stimulation. Hamidi and Mortazavi (2012) used 3DEC for dimensionally small models to study the effects of different fracture fluid properties and fluid rate, in-situ stress state, and rock mass properties on hydraulic fracture propagation.  Modeling three dimensional representation of hydraulic fracturing using 3DEC is computationally intensive and requires a very long time running. Thus UDEC represents a medium which can more efficiently capture the influence of geology without being as much computationally heavy. 1.4.5.1.2 Bonded particle method Bonded particle methods (Potyondy and Cundall, 2004) initiated from the Discrete Element Method application to a discontinuous mass represented by two dimensional discs or three dimensional spheres. The key idea of the bounded particle method is that the rock mass can be approximated by a bounded mass of cemented particles (Figure 1.5). The particles are assumed to be rigid. The interactions between the particles are through their contacts in a way that deformation is created at the particle contacts or by relative displacements between   23 particles (Bobet, 2009). Tensile and/or shear cracks can be simulated between particles where the tensile or shear strength of the contact is reached. The Particle Flow Code (PFC; Itasca Consulting Group, 1999) is a commercially available 2-D and 3-D bonded particle code. It can be used in analysis, testing, and research in different fields for presenting the interaction of many discrete objects that undergo large-strain and/or fracturing. This allows the modeling of hydrofracturing to be approached by assuming that the rock is made up of individual particles with specific stiffness assigned to them bonded together with bonds with specific strength. Under the applied load, the bonds between particles can break and develop fracture. Gil et al. (2010) modeled a hydraulic fracturing using the PFC code to show that the model is capable of capturing the true physics of injection into low permeability formation. Damjanac et al. (2010) modeled a hydraulic fracturing with a field-driven Discrete Fracture Network (DFN) superimposed on an intact rock particle model as a new approach to hydraulic fracturing modeling using PFC.  In PFC, the mechanical properties (stiffness and strength) of the particles and the bonds are selected by numerical experiments to match the overall mechanical behavior of a given rock. It is often difficult to choose material properties for PFC model in a way that the behavior of the resulting synthetic material, consisting the particles, shows similar physical behavior of chosen material. In PFC code, which synthesizes macroscale material behavior (referred to the scale of the model) from the microscale components interactions (referred to the scale of the particles), the microscopic elements properties are unknown. Therefore, first the response of the physical material must be determined, then using calibration process, the appropriate micro-properties can be chosen to be used in the model. In this process, the   24 responses of the synthetic material used in a numerical model are directly compared with the responses of the material in the lab or in the field.  The main drawback of this code is that the mechanical properties of the particles and the bonds are selected by numerical experiments to be used as the inputs for a numerical experiment. The rock mass is conceptually presented as a microscale granular assemblage. Therefore, in bigger models which have significantly more elements, the computation time increases. This considerable computational time is noted in Gil et al. (2010), where they changed the model size to a smaller one when the model reached quasi-steady state, in order to minimize simulation time. 1.4.5.1.3 Finite element/Discrete element (FEM/DEM) Hybrid Finite Element/Discrete Element (FEM/DEM) code utilizes both finite elements and discrete elements, for simulating fracture propagation, strain, and block kinematics (Munjiza, 2004). ELFEN (Rockfield Software, 2009) is an advanced Finite Element/Discrete Element commercial code. This technique includes brittle fracture capabilities enabling the modeling of a continuum passing to a discontinuum in response to changes in stress and strains. The code provides linear elastic, critical state, and plastic constitutive models for different geological materials. The fracture propagation is governed by coupled rock and fluid mechanics. In a recent study, Rogers et al. (2010) used ELFEN to do a series of 2-D simulations with embedded Discrete Fracture Network (DFN). They performed a number of simulations on generic fracture models by differing properties such as fracture length, aperture, and intensity, and studied the impact of those changes on the resultant micro-seismic pattern. The disadvantage of the code is that it is computationally intensive and requires more significant time for performing a coupled hydromechanical analysis.   25 1.4.5.2 Numerical modeling methodology In UDEC, the problem domain is converted into a system of interacting blocks through the generation of one or more joint sets to representing a system of natural fractures based on the in-situ geological structure. Key for simulating hydrofracturing, UDEC has the capability to model fluid flow through the fractures. UDEC performs a fully coupled hydro-mechanical analysis, in which mechanical deformation of joint apertures changes fracture conductivity, and joint water pressures affects the mechanical computations of joint aperture (Figure 1.6). The blocks in this assemblage are treated as being impermeable. The cubic law for flow within a planar fracture is used, where the flow rate is given by:                                                               lPkaq ?? 3                                                               (1.1) where, k is a joint conductivity factor, a is the contact hydraulic aperture, ?P is the pressure difference between the two adjacent domains, and l is the length assigned to the contact between the domains. Since the UDEC formulation is restricted to the modeling of fracture flow, it should be noted that leak-off along the fractures diffusing into the rock matrix is assumed to be zero (only leak-off into other fractures is considered). Furthermore, the fracture flow is idealized by means of a parallel plate model and cubic law, which disregards tortuosity. When a joint contact is broken, the fluid flows into the joint with a finite permeability. The hydro-mechanical coupling is completed by associating the aperture to the mechanical joint stiffness and acting effective normal stresses. The rock mass modeled in this study is represented by either two persistent orthogonal planes of weakness oriented parallel and perpendicular to bedding or by using the Voronoi tessellation generator to represent the rock features. These serve as incipient planes   26 along which hydrofrac propagation is restricted (as previously noted, the blocks are otherwise non-divisible) and are referred to as ?incipient fractures? as defined by Dusseault and McLennan (2011). The concept of incipient fractures is utilized to emphasize that it is not only the clearly visible fractures identified in televiewer, geophysical and core logs that are opened or sheared. Since, in the process of hydraulic fracturing, significant tensile stresses are generated, opening incipient fractures that can then serve as a path for hydraulic fracture propagation. Dusseault and McLennan (2011) describe two zones that develop around the hydrofrac injection point: an invaded zone and a dilated zone (Figure 1.7). The invaded zone is the region of propped fractures and the dilated zone is the region where the fractures are permanently opened and propped by shear displacement, or wedging and block rotation.  As mentioned earlier, one fundamental advantage of UDEC is that pre-existing joints in the rock mass can be incorporated directly and are allowed to undergo large deformations in shear or opening. The fracture mechanics techniques concentrate on tensile fracture mode, which is appropriate for hydrofrac, but misses equally important shear dilation in invaded and dilated zones.  Therefore, these zones can be seen in the UDEC model which accounts for both shear and dilation along existing joints and tensile rupture of rock in response to fluid injection and subsequent changes to the effective stress field. Generally, there are other advantages and some disadvantages in using UDEC for modeling hydraulic fracturing. The most important advantage of UDEC over commercial hydraulic fracture models is that UDEC treats the rock mass as a discontinuum as opposed to other commercial models which adopt oversimplifying assumptions of the geological environment, treating the rock mass as an equivalent continuum. The other advantage of   27 choosing UDEC over other commonly used tools for hydrofrac analysis is that UDEC allows large strain behavior as oppose to commercial models which are limited to small strain behavior. Also, UDEC is capable of coupled stress-flow modeling (i.e. coupled hydromechanical modeling) whereas the commercial models are not performing a coupled mechanical and hydraulic analysis.  There are also some disadvantages for choosing UDEC over other numerical approaches like 3DEC and ELFEN. The disadvantage of using UDEC versus 3DEC is that the three-dimensional nature of the process, represented by interconnectivity of both natural and hydraulically induced fractures in the 3D space is left out in a 2D model. It can be argued that even though the models are two-dimensional, the approach proves the concept of using a Discrete Element Method (DEM) as a modeling tool for hydraulic fracturing in naturally fractured rock mass. Also, the DEM is capable of modeling the behavior of fractured rock mass during hydraulic fracturing, which is the main focus of this thesis. ELFEN is also an advanced code for modeling hydraulic fracturing with having the option of simulation with embedded DFN, but it is requires significantly more running time. Thus UDEC represents a practical medium that is capable of adequately capturing the influence of geology, but does so without being too computationally heavy.   30 Table 1.1.  Comparison of common geotechnical numerical codes. Numerical Codes Coupled Fluid-Mechanical Analysis Fluid Flow Modeling Constitutive Models (representation of behaviour)  Natural Fractures (representation of geology) Induced Fractures Finite-Element/Finite Difference Capable Fluid flow is modeled through a continuum porous medium.                                         Full range of constitutive models, including plasticity, strain-dependent strength degradation and time dependency. Low stiffness assigned to an interface; only a small number can be included.  Simulated by associating yielding of elements as representing fracture damage of the modeled continuum. Distinct Element Capable Fluid flow is modeled through fracture flow between a system of impermeable blocks. Full range of constitutive models, including plasticity, strain-dependent strength degradation and time dependency. Modeled as contacts between blocks, which can open, slip and close. Can efficiently model a large number of randomly oriented joints. Simulated through the splitting and opening of contacts between blocks. Bonded Particle  Capable Fluid flow is modeled through a coupled fluid-particle simulation using the Navier-Stokes equation for a fluid interacting with a solid phase. Based on micro-mechanical contact properties, which are calibrated by simulating laboratory tests. Modeled as zero tensile strength contacts between particles. Requires advanced techniques to represent a joint network. Uses breakage of bonds between particles to represent induced fracturing.   31   Figure 1.1. Hydraulic fracturing in a vertical wellbore. The fracture propagation direction is perpendicular to the minimum horizontal stress.       32 Figure 1.2. Conventional (vertical) and unconventional (horizontal) hydraulic fracturing wells.    33 a) b) Figure 1.3. a) PKN fracture model, and b) KGD fracture model.   Figure 1.4. Conceptual model of a rock as used in the Distinct Element Method (DEM).   Figure 1.5. Conceptual model of a rock as used in the Bounded Particle Model (BPM).   34    Figure 1.6. Discontinuum model: distinct element discretization.   Figure 1.7. Conceptualization (plan view) of the invaded and dilated zones (after Dusseault and McLennan, 2011).         35 Chapter  2: Investigation of the Influence of Natural Fractures and In-situ Stresses on Hydrofrac Propagation using a Distinct-Element Approach 2.1 Introduction Hydraulic fracturing has been one of the primary means for improving well productivity in the oil and gas industry, over the last 60 years. Hydraulic fracturing design was originally applied assuming homogeneous and stiff reservoir rocks that behave as a linear elastic continuum. Significant advances have since been made in hydraulic fracturing theory and practice, based largely on linear elastic fracture mechanics; however, these too rely on a continuum treatment of the rock. Current standard hydrofrac simulations are capable of modeling the complex hydro-mechanical interactions involved in fracture propagation, projecting the hydrofrac as bi-wing, perfectly planar, symmetrical fractures. But again, this is done at the expense of making simplifying assumptions regarding the host geology within which the hydraulic fracture is propagating. In-situ, the reservoir rock conditions are complex considering the natural fractures. Natural discontinuities are present, including bedding planes, joints, and faults. The  influence of these pre-existing fractures on the propagation of the induced hydraulic fracture is often neglected in conventional hydrofrac models. Yet, one of the key elements of a successful hydrofrac job is the effective connectivity of the hydraulic fracture, which would interact with and incorporate the natural discontinuities in the rock mass. Further complexity is the superimposition of the in-situ field, which includes its own complexity and heterogeneity. Together, these have the potential to either enhance or diminish the effectiveness of the hydrofrac treatment and subsequent hydrocarbon production.    36 This study describes a series of numerical experiments investigating the influence of natural fractures on the modeling of hydraulic fracture propagation. The primary objective is to evaluate the use of 2-D distinct element techniques that account for both shear and dilation along existing joints and tensile rupture of intact rock in response to fluid injection and subsequent changes to the effective stress field. The analysis is carried out for several geometrical variants, related to hypothetical geological scenarios simulating a naturally fractured shale gas reservoir, and their effect on hydraulic fracture size, orientation and path. Also examined is the propagation limiting effects of leak-off to the natural fracture system and as a result of the leak-off, enhancing effects of reactivating natural fractures instead of creating new fractures the intact rock. 2.2 Hydraulic fracture modeling and shale gas reservoirs  The main purpose of reservoir stimulation is to enhance the reservoir permeability, assisting more efficient hydrocarbons delivery, and also to increase ultimate economic recovery. An evaluation of the economics of a hydraulic fracturing treatment includes: increase in production, treatment cost, and risks associated with mechanical problems. For these tasks, computer models can be used to predict the propagation of hydraulic fractures through the porous media. In the two introductory models upon which early hydraulic fractures were designed, the PKN (Perkins and Kern, 1961) and KGD (Geertsma and de Klerk, 1969), the rock mass is assumed to be a continuous, homogeneous, isotropic, linear elastic solid. Many subsequent programs specific to fracture modeling have likewise repeated this assumption, treating the rock mass as an elastic continuum (e.g. Barree 1983, Shaffer et al. 1984, Ventura 1985,   37 Meyer 1989, Mohaghegh et al. 1999, Dong and Pater 2002, Reynolds et al. 2004). In actuality, however, the geology of the reservoir rocks, especially those related to unconventional resources (e.g. shale gas), involve heterogeneous, anisotropic rock masses containing numerous natural discontinuities and other defects that serve as planes of weakness.   Geological discontinuities such as joints, bedding planes, shear fractures and faults are found commonly in gas shale reservoirs (e.g. Figure 2.1). Field experiments in gas shales suggest a significant interaction between natural fractures and induced hydraulic fractures. Thus the influence of natural fractures on the final performance of a hydraulic fracture treatment is significant and needs to be considered in its design.  To meet the objective of considering the effect of natural fractures on hydrofrac performance, it is desirable to include in the model a large number of oriented discontinuities to more accurately represent the planes of weakness that exist within the reservoir rock masses. These include weak bedding planes and persistent joints forming along bedding that typically represent the dominant discontinuity set. Other non-persistent natural fracture sets would also be present.  The changes that occur during burial (diagenesis) can create a non-persistent set of cross joints that are sub-normal to bedding. There could also be present a set of fractures at acute angle that are fold-related as a result of tectonic forces (Nelson, 1985). Together, these discontinuity sets form a network of pre-existing fractures of varying persistence, spacing, dip angle and direction relative to the orientation of the in-situ stress field and injection wellbore.    38 2.3 Numerical modeling methodology The 2-D commercial code UDEC (Universal Distinct Element Code; Itasca Consulting Group, 1999) is used here to simulate the response of a jointed rock mass subject to static loading and hydraulic injection. The Distinct Element Method (DEM) is a Lagrangian numerical technique in which the problem domain is divided through by discontinuities of variable orientation, spacing and persistence, which may undergo large deformations in shear and opening in response to a change to the modeled stresses. UDEC can simulate the progressive failure associated with crack initiation and propagation by the breaking of pre-existing contacts between the pre-defined joint bounded blocks. The blocks are deformable but they remain intact. A detailed discussion of the Distinct Element Method is presented by Cundall and Hart (1992). A limiting factor in the modeling of stress-induced fracturing using UDEC is that all potential fracture pathways must be predefined. To counter this limitation and provide added degrees of freedom for fracture propagation, a Voronoi tessellation scheme can be used to generate randomly sized polygonal blocks (Figure 2.2). These Voronoi contacts can then be assigned strengths equivalent to weak intact rock, and are referred to here as ?incipient fractures?. These are used in this study to provide a background network of potential fracture pathways to model the propagation path of a hydraulically driven fracture. Superimposed on this can be discontinuities of a set dip angle, spacing and persistence, assigned zero cohesion and tensile strength values, representing natural bedding and cross joints in the rock mass. This integrated system was developed and tested here to enable the modeling of the interactions between the natural fracture network and induced fracture as directed by the in-situ stress field.    39 Key for simulating hydraulic fracturing, UDEC is capable of performing a fully coupled hydro-mechanical analysis to model fluid flow through a network of fractures. Mechanical deformation of joint apertures changes conductivity and, conversely, the connectivity changes the joint water pressure, which affects the mechanical computations of joint aperture (Figure 2.2). The blocks in this assemblage are treated as being impermeable. The cubic law for flow within a planar fracture is used, where the flow rate is given by:                                                        lPkaq ?? 3                                                                           (3.1) where, k is a joint conductivity factor, a is the contact hydraulic aperture, P? is the pressure difference between the two adjacent domains, and l is the length assigned to the contact between the domains. Since the UDEC formulation is restricted to the modeling of fracture flow, it should be noted that leak-off along the fractures diffusing into the rock matrix is assumed to be zero (only leak-off into incipient fractures is considered). Furthermore, the fracture flow is idealized by means of a parallel plate model and cubic law, which disregards tortuosity. When a joint contact is broken, the fluid flows into the joint with a finite permeability, k:                                                            ?121?k                                                                                (3.2) where, k is the joint conductivity factor and ?   is fluid dynamic viscosity. The hydro-mechanical coupling is completed by associating the hydraulic aperture to the mechanical joint stiffness and acting effective normal stresses (i.e., mechanical aperture).    40 2.4 Influence of natural fractures in hydraulic fractures Hubbert and Willis (1957) showed that hydraulic fractures in an isotropic medium always propagate perpendicular to the orientation of the minimum principal stress. In other words, the fracture propagates in the direction for which opening is most easily accommodated, propagating parallel to the maximum stress so that opening is against the minimum principal stress. However, in a non-isotropic medium, hydrofrac mine-back experiments by Warren and Smith (1985) have shown that pre-existing fractures and faults can have an influence on fracture propagation, although the overall trajectory is still controlled by the orientation of the minimum principal stress, as also noted by Zoback (2007).  Blanton (1982, 1986) showed the interaction of hydraulic fractures with natural fracture systems under different angles of approach and applied triaxial stress fields in a series of laboratory experiments on Devonian shale as well as blocks of pre-fractured hydrostone. These experiments showed three classes of induced hydraulic fracture behavior whereby the propagating fracture: i) crosses the pre-existing fracture; ii) is diverted by the pre-existing fracture; or iii) is arrested by the pre-existing fracture. The behavior was found to be dependent on the background stress conditions and the angle of approach to the pre-existing fracture. Warpinski and Teufel (1987) presented results investigating the influence of geological discontinuities on the propagation of hydraulic fractures observed in mine-back experiments. They likewise confirmed that the geological discontinuities influence the overall geometry and effectiveness of the hydraulic fracture.  They also showed that when a hydraulic fracture intersects a natural fracture or other discontinuity at some angle with   41 respect to the hydraulic fracture (or the direction of maximum stress), the fracture will propagate across the natural fracture, or activate and dilate the natural fracture in shear. Olson et al. (2012) examined and showed the effects of cemented (or healed) natural fractures, under different angles of approach to the hydraulic fracture direction. These tests were done in a series of laboratory experiments on cast hydrostone containing embedded planar glass discontinuities as proxies for cemented natural fractures. These experiments showed that the induced hydraulic fracture can bypass the natural fracture by propagating around it, be arrested by it or diverting along it, or a combination of bypass and diversion.   Dershowitz et al. (2010) noted the significance of the interaction of natural fractures with the propagating fracture based on Discrete Fracture Network (DFN) realizations. They studied the influence of natural fractures on hydraulic fracture propagation in the cases where the distributed natural fractures cross the path along which the hydrofrac propagation is to occur and where other natural fractures are in proximity to the hydrofrac propagation path. 2.5 UDEC simulations 2.5.1 Invaded and dilated zones The rock mass modeled in this study includes both natural fractures and incipient fractures as defined by Dusseault and McLennan (2011). Construction of the UDEC models began with the generation of the incipient fractures using the Voronoi tessellation generator (Figure 2.3). The concept of incipient fractures is utilized to emphasize that it is not only clearly visible fractures identified in televiewer, geophysical and core logs that are opened or sheared, but that incipient fractures can also serve as potential pathways for hydraulic fracture propagation.   42 Figure 2.4 shows the result of numerical simulation of fluid injection in the model shown in Figure 2.3. The results here show the initiation and propagation of the hydraulic fracture, with fracture development indicated by increasing fracture aperture during fluid injection into the wellbore. Fluid flow is accommodated along the joint network by assuming a hydraulic conductivity based on the aperture of the joints and the viscosity of the permeating fluid. Strength and deformation properties for the incipient fractures are provided in Table 2.1. The stress field is superimposed across the problem domain to represent the in-situ stresses at depth as shown in Figure 2.3.  In the process of hydrofracing, the permeability of the modeled rock mass is enhanced. Dusseault and McLennan (2011) describe two zones that develop around the wellbore: an invaded zone and a dilated zone (Figure 2.5). The invaded zone is the region of propped fractures and the dilated zone is the region where the natural fractures are permanently opened by wedging, block rotation, or dilation by shear displacement (Dusseault and McLennan, 2011). These zones are seen in the UDEC model results in plots of the increased fracture aperture (Figure 2.4), simulated pore pressure distribution (Figure 2.6) and fractures that shear and dilate in response to the corresponding change in effective stresses (Figure 2.7). Figure 2.4 shows that in addition to the main fracture pathway, in which the tensile incipient fracture strength has been exceeded and invaded by the injected fluid, there is a zone in which the injected fluid dilates neighboring incipient fractures. This agrees with the conceptualization of the dilated zone developed by Dusseault and McLennan (2011). The simulated hydraulic fracture in Figure 2.4 and pore pressure distribution in Figure 2.6 show the dilated zone has enhanced flow properties with fluid flow in the incipient fractures. By   43 applying the fluid pressure in the wellbore, the initial model response to fluid injection is tensile hydraulic fracturing, followed by a combination of tensile opening and shear dilation in the neighboring incipient fractures (Figure 2.7). Figure 2.8 shows the displacement vectors corresponding to the pressurized wellbore and induced hydraulic fracture. The displacement vectors indicate that the main mode of deformation is fracture opening with the incipient joints along the hydraulic fractures failing in tension (Figure 2.8).   A different response is observed when a persistent discontinuity, for example a fault, is added to the rock mass in near proximity to the wellbore. Figure 2.9 shows the model geometry for this scenario (properties for the fault are given in Table 2.1). The wellbore is pressurized, and as a result, a hydraulic fracture is generated, which intersects the fault and induces slip along the fault (Figure 2.10). Figure 2.11 shows the shear displacements along the fault, indicating that the main mode of deformation is shear dilation. 2.5.2 Influence of natural fractures on initiation of hydraulic fractures As demonstrated in the previous section, the complexity of the hydraulic fracture and surrounding dilated zone is dependent on the presence of incipient and natural fractures populating the shale formation. These fractures are ubiquitously distributed throughout the rock mass. To investigate this closer, UDEC was utilized here to simulate the interactions of hydrofrac with the natural fractures which are in the hydrofrac propagation path and the natural fractures which are in the proximity to the hydrofrac propagation path, described by Dershowitz et al. (2010).  The UDEC experiment involves a set of non-persistent natural fractures distributed in the rock mass. The simulation of hydraulic fracture propagation is conducted based on a rock   44 mass model where the distributed natural fractures cross the path along which the hydrofrac propagation is to occur, as dictated by the orientation of the principal stresses (Figure 2.12). The induced fracture propagates in the direction perpendicular to the minimum in-situ stress, and coalesces with one of the natural fractures as it propagates. Although the overall trajectory of the induced hydraulic fracture is controlled by the in-situ stresses, the natural fracture in the proximity of the propagation path is seen to have an influence on the induced fracture. This effect is shown in the simulated pore pressure distribution (Figure 2.13). This figure indicates that the pore pressure response is asymmetric, with higher pressures developing to the right of the wellbore where the hydraulic fracture path directly intersects a natural fracture. This natural fracture facilitates the advancement of the hydraulic fracture (Figure 2.14), diverting the injection fluid due to its lower resistance to aperture opening (i.e. higher hydraulic conductivity). Further interrogation of the pore pressure distribution (Figure 2.13) indicates that the presence of other natural fractures in proximity to the hydrofrac propagation path serve to contribute to the development of the invaded and dilated zones (Figure 2.14). Therefore, in addition to their influence on the hydrofrac propagation path, the presence of the natural fractures also influences the dilated zone (Figure 2.15).  2.5.3 Effect of differential stress and angle of approach on interaction between hydraulic fractures and natural fractures  Blanton (1982, 1986) showed that the interaction of a hydraulic fracture with a natural fracture under different angles of approach (i.e., the angle between the propagating hydraulic fracture and natural fracture it?s approaching; shown in Figure 2.16), may result in one of   45 three responses. Firstly, if the injection pressures are sufficient, the hydraulic fracture may propagate across the natural fracture without changing direction, with the hydrofrac remaining essentially planar. The second scenario is where the hydraulic fracture does not cross the natural fracture but fluid pressures are sufficient to offset and initiate a continuation of the hydrofrac from the tip of the natural fracture. The third scenario is where the hydrofrac may stop propagating (arrest) after it approaches and intersects the natural fracture; in this case, the fluid pressures drop below those needed for the hydrofrac to continue due to leak-off into the natural fracture. Whether the hydraulic fracture crosses or doesn?t cross the natural fracture, depends on the in-situ stress conditions relative to the orientation and properties of the natural fracture. A series of UDEC experiments was carried out to investigate the types of interactions that can be expected between induced and natural fractures under various angles of approach and states of stress. Table 2.2 summarizes the differential stress and angle of approach under which each case has been modeled. 2.5.3.1 Crossing interaction An example where the UDEC simulation produced a hydraulic fracture that crossed a pre-existing natural fracture is shown in Figure 2.17. In this case, the hydraulic fracture approached the natural fracture with a high incidence angle of 60 degrees under high differential in-situ stresses (?1 ? ?3 > 11 MPa). The results presented in Figures 2.18 and 2.19, show the induced fractures (via incipient fractures that have opened) and simulated pore pressure distributions, respectively. The model suggests that the high angle of approach between the induced fracture and natural fracture results in a crossing interaction because the   46 pressure required to redirect the hydrofrac and open in a direction non-parallel to the minimum in-situ stress is greater than that for the hydrofrac to continue along the same path as if the high angled natural fracture wasn?t there (despite the natural fracture being weaker than the incipient fractures). At the same time, the injection pressure is sufficient despite leak-off for the hydraulic fracture to continue propagating instead of arresting. The induced fracture in Figure 2.18 indicates its propagation is bi-wing, propagating in both directions from the wellbore. Further interrogation of the pore pressure distribution (Figure 2.19) likewise indicates a bi-wing pattern with pore pressures increasing on both sides of the wellbore. The tensile and opening response of the model to fluid injection is shown in Figure 2.20. 2.5.3.2 Offset interaction The second scenario where a hydraulic fracture approaches, propagates along and then re-continues along an offset path from the tip of the natural fracture was observed when the approach angle was less than 30 degrees and the differential in-situ stresses were low (Figure 2.21). The results presented in Figures 2.22 and 2.23, show the induced fractures and simulated pore pressure distribution, respectively. The model suggests that the low angle of approach between the hydrofrac and natural fracture under a low in-situ differential stress (?1 ? ?3 < 7 MPa) favored a path re-routed along the lower strength natural fracture (zero cohesion and tensile strength) even though it was not optimally oriented with respect to the minimum principal stress. This compares to the alternative crossing path involving an incipient fracture more optimally oriented relative to the minimum principal stress but with a   47 finite cohesive and tensile strength. The approach angle was low enough in this case that the hydraulic fracture favored the weakest path. The fracture pattern also suggests that the opening of the natural fracture acts to limit the effectiveness of the hydraulically induced fracture as a result of fracture fluid being diverted to the natural fracture. The hydraulic fracture (Figure 2.22) and pore pressure distribution (Figure 2.23) also indicate that the response is asymmetrical as a result of this interaction, with leak-off occurring into the natural fracture resulting in a localized increase in the dilated zone shown in Figure 2.24. 2.5.3.3 Arresting interaction Simulations for the third scenario indicating an arresting response was observed when the hydraulic fracture approached the natural fracture with an angle of 45 degrees under intermediate differential in-situ stress conditions (Figure 2.25). Results presented in Figures 2.26 and 2.27, show the induced hydrofrac and simulated pore pressure distributions, respectively. The fracture pattern in Figure 2.26 suggests that the 45 degree angle of approach is too sharp for the hydrofrac to cross but is not inclined at a suitable angle relative to the minimum principal stress for the fluid pressures to open the natural fracture and allow the hydrofrac to continue by initiating at the crack tip (i.e., offset interaction). The result is that the arresting of the hydraulic fracture acts to significantly limit the effectiveness of the treatment. Figure 2.27 shows that the propagation of the hydraulic fracture is asymmetrical as a result of the arrested interaction; a continuum analysis would inherently predict a symmetric hydrofrac around the wellbore. Further interrogation of the pore pressure distribution (Figure 2.27) indicates higher pore pressures increase around the   48 wellbore in the dilated zone which ultimately could cause branching and growth of the hydrofrac in that zone. The tensile and opening response of the model to fluid injection is shown in Figure 2.28. The branching would be facilitated by the presence of the ubiquitous natural fracture network as shown in Figure 2.14.  In this case, hydrofrac arrest occurred under intermediate differential stresses (?1 ? ?3 ? 9 MPa) and an intermediate angle of approach (45 degrees). Arresting could be a temporary interaction, and by continued pressure increase, re-routing (offset) or crossing could occur. When the injection pressure exceeds the normal stress acting across the natural fracture (i.e., zero effective normal stress), the natural fracture will open, slip and dilate, and fluid pressure will leak off into the natural fracture. With continued injection, the hydrofrac may then re-initiate once the pressure required to open the adjacent incipient fractures is exceeded. This staged interaction results in the hydraulic fracture length on the two sides of the wellbore differing as well as the degree of hydrofrac branching in the dilated zone. 2.5.3.4 Bounding relationships as a function of approach angle and stress state  Based on the scenarios described above, a series of sensitivity analyses were performed to provide a basis for selecting fracturing treatments that is most effective in connecting the most natural fracture system to the wellbore. Figure 2.29 shows the results of a series of numerical experiments applying several different angles of approach and states of stress using the bisection approach, which is a variation of the incremental search method (Chapra, 2008). The results show clear zones of opening and crossing type interactions, separated by a zone of arresting interaction in between the opening and crossing limits. The cases representing the three scenarios described in the previous sections are highlighted (Figure   49 2.29). The results of UDEC numerical experiment for each point presented in Figure 2.29 is provided in the Appendix. It should also be emphasized that the analyses carried out here were based on a conceptualized model of a shale gas reservoir. The same numerical study can be performed for any given reservoir to determine site specific correlations between the planned hydrofrac and natural fracture sets present in the basin under study. In different zones of a single basin, different fracture orientations can be expected and corresponding geo-risk assessed based on the reservoir in-situ stress conditions. The numerical modeling procedure outlined here would provide a valuable tool for decision makers to assess the potential for crossing, offsetting or arresting interactions, and the potential for asymmetrical hydrofrac treatments. 2.6 Discussion and conclusions The purpose of this study has been to investigate the application of distinct element numerical modeling techniques to model the interactions between a hydraulic fracture and a network of natural fractures. The discontinuum-based approach showed value in modeling the complex behavior of hydraulic fracture propagation in naturally fractured reservoirs. The model was capable of simulating the tensile failure and shear dilation of incipient fractures in the intact rock and along reactivated natural fractures. Furthermore, the models were also able to simulate the interconnections that develop adjacent to the hydrofrac in the invaded and dilated zones. The value of using a discontinuum- over a continuum-based approach was demonstrated through its ability to model asymmetric hydraulic fracture behaviors in response to the influence of geological discontinuities. The importance of accounting for natural discontinuities in the invaded zone was investigated and the results show a   50 considerable influence in that zone with respect to enhanced permeability arising from branching and dilation of fractures adjacent to the propagating hydrofrac.  Interactions between the induced hydrofrac and existing fractures were further examined with respect to different angles of approach and states of stress. The numerical results showed that hydraulic fractures tend to cross pre-existing fractures under high differential stresses and high angles of approach. At low differential stresses and angles of approach, pre-existing fractures tend to open and divert fracturing fluid into the pre-existing fracture resulting in a re-routing of the hydrofrac such that it offsets and re-initiates and propagates from the tip of the natural fracture. An intermediate differential stress and angle of approach was seen to result in the arrest (propagation stopped) of the hydraulic fracture.  The implication of this study for hydraulic fractures in the field is that the orientation of the fracture network relative to that of the in-situ stress field plays an important role with respect to the extent and symmetry of hydraulic fracture propagation, and thus the overall success of achieving the reservoir coverage desired. Symmetrical fracture extension on both sides of a wellbore is unlikely in the presence of a larger discrete fracture such as fault because of the different stress and coupled hydromechanical interactions that arise. Exceptions may arise where high differential stresses and a high angle of approach results in the hydraulic fracture crossing the natural fracture. The presence of natural fractures may further act to redirect fluid pressures requiring longer pressurization intervals to achieve the desirable fracture length, with consideration given to pressure/fluid leak-off into adjacent natural fractures. In the worst case scenario, the presence of a natural fracture (under intermediate differential stresses and intermediate angle of approach) may cause the hydraulic fracture propagation to stop. As a result, a longer pressurization and/or higher   51 pumping pressures may be required to re-initiate the hydrofrac to achieve the desirable fracture length. In different parts of a basin, the dominant fracture network angle could spatially vary  (e.g. Appalachian Basin; Cliffs Minerals Inc., 1982). As a result, different interactions can be expected across different zones in the basin. This introduces an element of geo-risk; the orientation of the natural fractures relative to the stress field could adversely affect the propagation of the hydraulic fracture, possibly limiting its extension beyond the first set of natural fracture it encounters. The extension of connection between wellbore and natural fractures can be achieved by further pumping, which may be of value in some reservoir treatments (Blanton 1986). Therefore, in different zones of a basin, different decisions may be required regarding operational factors based on an understanding of the type of potential interactions that may develop between the hydrofrac and natural fractures.  The purpose of this study was to provide a conceptual basis for selecting fracturing treatments that will result in an effective connection between the hydrofrac, natural fracture system and wellbore. Future work will involve applying these techniques to an actual gas shale reservoir, ground-truthing the modeling results using field data (e.g., microseismic data). Future work will also involve a closer examination of other parameters such as natural fracture property values and other limiting assumptions such the role of fluid viscosity, natural fracture aperture, tortuosity, the use of proppant, and extending the models to three dimensions.     52 Table 2.1. Properties assigned to the modeled planes of weakness and natural fractures. Discontinuity property Incipient fractures Natural fractures Units Friction angle 30 25 Degrees Residual  friction angle 25 20 Degree Cohesion 1.0 0.0 MPa Residual cohesion 0.0 0.0 MPa Tensile strength 0.5 0.0 MPa Residual tensile strength 0.0 0.0 MPa Normal stiffness 1x104 1x103 MPa/m Shear stiffness 1x103 1x102 MPa/m   Table 2.2. Assigned stresses and angle of approach to investigate hydraulic fracture behavior with respect to natural fracture. Interaction behavior Smax (MPa) Smin (MPa) Differential stress (MPa) Angle of approach (Degrees) Crossing 30 19 11 60 Offsetting 30 23 7 30 Arresting 30 21 9 45   53  Figure 2.1. Examples of natural fractures present in shales: a) Devonian shale outcrop from northeastern British Columbia; b) televiewer image from a horizontal well through the Maskwa/Otter Park formation; c) open and partially open fractures in retrieved core; and d) closed fractures in core samples from the Muskwa/Otter Park formation. After Reine and Dunphy, 2011.     54   Figure 2.2. Discontinuum model: distinct element discretization using Voronoi tessellation.   Figure 2.3. Rock mass model with background incipient fractures only, generated using Voronoi tessellation (20 x 20m).    55   Figure 2.4. Increased fracture aperture showing simulated hydraulic fracture initiation pattern for a reservoir rock comprised of incipient fractures only.   Figure 2.5. Conceptualization (plan view) of the invaded and dilated zones (after Dusseault and McLennan, 2011).     56  Figure 2.6. Simulated pore pressure distribution around the borehole.   Figure 2.7. Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red.     57  Figure 2.8. Displacement vectors around the pressurized borehole.    Figure 2.9. Rock mass model including background incipient fractures superimposed by a persistent fault.    58  Figure 2.10. Displacement vectors around the pressurized borehole and along the fault.   Figure 2.11. Corresponding shear displacement along the fault in response to the effective stress change.      59     Figure 2.12. Distinct element rock mass model with the distributed natural fractures. Zoomed in window provided above.            60    Figure 2.13. Simulated pore pressure distribution showing influence of natural fractures on development of the invaded and dilated zones. Zoomed in window provided above.                61     Figure 2.14. Modeled apertures showing interactions that develop between the induced hydrofrac and natural fractures, including the development of invaded and dilated zones. Zoomed in window provided above.    62  Figure 2.15. Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red.  Figure 2.16. Diagram of hydraulic fracture intersecting natural fracture.   63   Figure 2.17. Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 60 degrees angle of approach.          64   Figure 2.18. Modeled response of the induced hydrofrac (i.e., incipient fractures that have opened) as it intersects a natural fracture at a 60 degree approach angle under high differential stresses. In this case the hydrofrac crosses the natural fracture.  Zoomed in window provided above.          65    Figure 2.19. Simulated pore pressure distribution for scenario assuming a natural fracture at a 60 degree approach angle under high differential stresses. Zoomed in window provided above.        66   Figure 2.20. Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above.    67  Figure 2.21. Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 30 degree angle of approach.           68     Figure 2.22.  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 30 degrees approach angle under low differential stresses. In this case the hydrofrac path is offset by the natural fracture. Zoomed in window provided above.          69   Figure 2.23. Simulated pore pressure distribution for scenario assuming a natural fracture at a 30 degree approach angle under low differential stresses.  Zoomed in window provided above.               70    Figure 2.24. Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above.   71  Figure 2.25. Rock mass model (200 x 200m) with background of incipient fractures superimposed with a natural fracture at a 45 degree angle of approach.                  72   Figure 2.26.  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 45 degree approach angle under intermediate differential stresses. In this case the hydrofrac is stopped by the natural fracture. Zoomed in window provided above.           73   Figure 2.27. Simulated pore pressure distribution for scenario assuming a natural fracture at a 45 degrees approach angle under intermediate differential stresses.  Zoomed in window provided above.                74    Figure 2.28. Joints that open in response to the modeled hydrofrac  (zero normal stress), shown in green, and fractures that shear and dilate in response to the change in the effective stresses, shown in red. Zoomed in window provided above.       75  Figure 2.29. Correlation between the angle of approach and differential stress based on UDEC based numerical experiments.                 Opening upper bound Crossing lower bound     76 Chapter  3: Investigation of the Influence of Stress Shadows on Multiple Hydraulic Fractures from Adjacent Horizontal Wells  3.1 Introduction Unconventional reservoirs such as gas shale require technology-based solutions for optimum development. The successful exploitation of these reservoirs has relied on technological advances in horizontal drilling, multiple stage completions, innovative fracturing, and fracture mapping to engineer economic completions.  Towards this purpose, hydraulic fracturing serves as the primary means for improving well productivity.  Simulations demonstrate that shale with ultra-low permeability requires an interconnected fracture network with relatively small spacing between fractures to obtain a reasonable recovery factor (Warpinski et al., 2008). Hence, multiple hydraulic fractures from one or more horizontal wellbores provide an effective means to maximize the fracture network surface area. Horizontal wells provide greater contact with a larger volume of reservoir rock than vertical wells due to the limited vertical thickness of the reservoir compared to its lateral extent.  Recent studies have suggested that simultaneous or near simultaneous hydraulic fracturing of adjacent wells results in better well performance than fracturing adjacent wells sequentially (Vermylen and Zoback, 2011; Kennedy et al., 2012). This has evolved into the drilling of multiple horizontal wells in an attempt to maximize the stimulated volume of reservoir rock through hydraulic fracturing. However, changes to the in-situ stress field caused by an earlier hydrofrac on subsequent hydraulic fractures, referred   77 to here as the ?stress shadow? effect, are not accounted for in conventional hydraulic fracturing design calculations. These stress shadow effects are potentially critical to the design of multiple hydrofracture treatments. Because of the stress shadow effects, multiple hydraulic fractures, whether off a multi-stage single well or adjacent horizontal wells should not be designed identical to a single horizontal well treatment. Thus the effects of stress shadows around an initial horizontal well must be considered with respect to their influence on the performance of a hydraulic fracture from a neighboring well (e.g. Wong et al. 2013).  This chapter describes a series of numerical experiments investigating the influence of stress shadow on multiple hydraulic fracture treatments. One objective is to evaluate the utility of 2-D distinct element techniques to model shear slip and aperture dilation along existing rock joints and tensile rupture of intact rock in response to fluid injection and subsequent changes to the effective stress field.  The analysis is carried out for different completion scenarios to investigate their effect on the propagation of hydraulic fractures. The changes in fluid pressure and corresponding effective stress changes around each wellbore during different completion techniques are examined. The effects of multiple hydraulic fractures on stress shadowing are also studied as a function of the reservoir depth, in-situ stress ratio, wellbore spacing, and injection rate.     3.2 Influence of stress shadows arising from multiple hydrofracs Altered-stress fracturing is a concept whereby a hydraulic fracture in one well is affected by another in a nearby adjacent well. One of the earliest studies was done by Warpinski and   78 Branagan (1989), where they presented field tests and finite element calculations examining the modified stress field around a wellbore. Termed ?stress shadows?, the disturbance of the stress field is especially important when considering a multiple stage hydrofracing design. When sequential hydrofrac stages are initiated in horizontal wells that are close to one another, the stress perturbation caused by one may affect subsequent hydraulic fractures.  Different authors have examined stress perturbation in multi-stage fracturing from horizontal wellbores to optimize horizontal completion techniques. For example, Fisher et al. (2004) presented results investigating the influence of multi-stage fractures with a wellbore separated into equal sections in Barnett shale using microseismic data. They showed that stress diversion is present when the reservoir has been supercharged by a previous fracture treatment stage. They also showed that stress in this region is increased due to locally higher fluid pressures, which influence subsequent stages. Morrill and Miskimins (2012) performed a series of numerical simulations of stresses around a fracture tip in a multiple hydraulically fractured horizontal well in order to determine the optimal fracture spacing to avoid stress field interactions and allow for predictable fracture geometries and conductivities in shale gas. Daneshy et al. (2012) reported the results of pressure measurements carried out in four horizontal wells where two of the wellbores served as observation wells throughout the project while the other two were actively being fractured. The motivation for these measurements was the detection of fracture shadowing created through their extension. They showed that field measurements can be incorporated into the development plan and concluded that real-time monitoring gives the operator time to optimize the treatment and modify future designs, accordingly.   79 A recent consideration of stress shadow effects is when closely spaced multiple horizontal wells are used. Vermylen and Zoback (2011) studied stress shadow effects in multiple horizontal wells in the upper Barnett shale for different completion procedures (simulfrac and zipperfrac) to test the effectiveness of different fracture methods. Simulfracs involve pressurization of two adjacent horizontal wellbores at the same time; zipperfracs involve first injecting (fracing) from one wellbore while the neighboring wellbore is not active and then injecting from the neighboring wellbore after injection into the first wellbore was completed. Vermylen and Zoback (2011) compared the activity level of a fracture stage for the different completion techniques using microseismic events. They found significant differences in stimulation outcome for the different fracing procedures owing to stress shadow effects. Nagel and Sanchez-Nagel (2011) performed a numerical evaluation of the effect of multiple hydraulic fractures on stress shadowing as a function of in-situ stress and operational factors. Wu et al. (2012) presented results of their study on stress shadow effects of multi stage fracturing from a horizontal well, showing that fractures can either enhance or suppress each other depending on their initial relative positions.  They concluded that accounting for these factors and their effects provides a means to optimize shale completions. Previous studies suggest only the influence of stress-shadows generated by hydraulic fracturing on subsequent hydrofracs, whether off a multi-stage single well or adjacent horizontal wells. Further study is required to investigate the role of stress perturbation in multi-stage fracturing from multiple horizontal wellbores towards optimization of horizontal completion techniques by comparing the stimulated volume for different completion procedures. Also, investigation of the influence of other factors, including the local in-situ   80 stress and operational factors, on the hydraulic fracturing effectiveness is still required when the reservoir has been supercharged by a previous fracture treatment stage. 3.3 Numerical modeling methodology The Distinct Element Method (DEM) is a Lagrangian numerical technique in which the problem domain is divided through by discontinuities of variable orientation, spacing and persistence (Cundall and Hart, 1992). Figure 3.1 provides an idealization of a DEM discretization of a problem domain and representation of the hydromechanical interactions between neighboring blocks. One fundamental advantage of the DEM is that pre-existing joints in the rock mass can be directly incorporated, providing the freedom for the problem domain to undergo large deformations through shear or opening along the discontinuities. The 2-D commercial code UDEC (Universal Distinct Element Code; Itasca Consulting Group, 1999) is used here to simulate the response of a jointed rock mass subjected to static loading and hydraulic injection. One advantage of UDEC over other available continuum-based hydraulic fracture codes is that the rock mass is treated as discontinuum therefore it helps to treat the geology in a more realistic way. Another advantage is the code is not limited to small strain behavior of joints and allows large strain behavior. UDEC is capable of modeling the progressive failure associated with crack propagation through the breaking of contacts between the pre-defined joint bounded blocks. The blocks are deformable but remain intact. Key for simulating hydrofracturing, UDEC has the capability to model fluid flow through the defined fracture network. A fully coupled hydro-mechanical analysis can be performed in which the mechanical deformation of joint apertures changes conductivity and, conversely, the connectivity changes the joint water   81 pressure, which affects the mechanical computations of joint aperture. The blocks in this assemblage are treated as being impermeable, and fracture flow is calculated using a cubic law relationship for joint aperture:                                                         lPkaq ?? 3                                                               (3.1) where, k is a joint conductivity factor (dependent on the fluid dynamic viscosity), a is the contact hydraulic aperture, P? is the pressure difference between the two adjacent domains, and l is the length assigned to the contact between the domains. Since the UDEC formulation is restricted to the modeling of fracture flow, leak-off along the fractures diffusing into the intact rock matrix is assumed to be zero (only leak-off into other fractures is considered). Furthermore, the cubic law flow assumption disregards tortuosity. When an incipient fracture contact is broken, the fluid flows into the new fracture. 3.4 Model set-up and simulation scenarios  The rock mass modeled in this study is represented by two orthogonal planes of weakness (Figure 3.2). These serve as incipient planes along which hydrofrac propagation is restricted (as previously noted, the blocks are otherwise non-divisible). The size of the blocks is increased towards the outer edges of the model. This is done for two reasons. First, it is important to have smaller blocks in the area of interest (where rock is fractured) in order to eliminate the influence of bock size on the induced fracture length. Second, the gradation to larger blocks towards the outer boundary helps to balance computing memory requirements against minimizing boundary effects on the results. The time required for executing each model in this study was typically on the order of 7 to 10 days.   82 The incipient planes of weakness were modeled assuming a Coulomb-slip constitutive model with both peak and post-peak properties; these are given in Table 3.1. The blocks were modeled as being elastic with a Young?s modulus of 30 GPa and Poisson?s ratio of 0.25. The stress condition in this study assumes a horizontal to vertical in-situ stress ratio, K, of 1.5 unless otherwise stated, representing a situation where locked-in horizontal stresses are present (as discussed in more detail below). The out of plane stress is assumed equal to the minimum stress and constant. The model (Figure 3.2) represents a 2-D vertical plane where the out-of-plane stress is horizontal. A vertical stress of 20 MPa was assumed (i.e. 1000 m depth). The gravitational variation of vertical stress from top to bottom of the model is neglected because the variation is small in comparison with the magnitude of stress acting on the volume of rock to be modeled. The initial background pressure of 10 MPa was applied. The wellbores pressurized in the model are horizontal and parallel to the out-of-plane direction (Figure 3.3). As a result, the alignment of the major principal stress with one of the two orthogonal planes of weakness facilitates the generation of a horizontal hydraulic fracture in the direction of the maximum stress. Boundary conditions are specified for the external boundaries of the model assuming constant stress and constant fluid pressure conditions. Three different hydrofrac scenarios were tested to investigate the influence of stress shadowing on hydraulic fracture stimulation:  ? Conventionally fractured well, where a single well (wellbore A or wellbore B) is pressurized.   83 ? Zipperfrac well, where first one well (wellbore B) and then an adjacent well (wellbore A) is pressurized. ? Simulfrac well, where the two wells (wellbore A and wellbore B) are pressurized simultaneously. The hydraulic fracture is simulated by applying a fluid injection assuming a constant flow rate of 5 m3/min into the wellbore. The fluid injection is discontinued at the shut-in time by setting flow rate to zero. 3.4.1 Stress shadowing simulation The in-situ stress field that exists at depth originates from gravitational, tectonic and remnant diagenesis stress components. These dictate the magnitudes and orientations of the principal stresses, which are often assumed to be aligned with the horizontal and vertical axes relative to the horizontal borehole. The in-situ stresses are known to control the direction of fracture propagation, and therefore represent the key boundary conditions influencing oil and gas reservoir stimulations. Significant perturbations to the stress field result from any process that changes the reservoir pressure and/or initiates and dilates fractures in the rock. To determine the influence of stress perturbations arising from an initial hydraulic fracture stimulation on subsequent hydraulic fractures from a neighboring injection well, three interaction scenarios as described by Vermylen and Zoback (2011) are modeled here. 3.4.1.1 Conventional hydrofrac (single well) The first scenario involved pressurizing a single horizontal wellbore (A in Figure 3.2) to generate a hydraulic fracture. Figure 3.4 shows the changes in fluid pressure at different points in the rock mass between the injection and a neighboring wellbore (Points 1 to 8; see   84 Figure 3.2). Points 8, 7, 6 and 5 are located at 5 m intervals from wellbore A, and points 1, 2, 3 and 4 are located at 5 m intervals from wellbore B. The latter (wellbore B) is inactive for this scenario but is included for future comparison with the other scenarios. Figure 3.4 shows that during injection, the pressure values recorded at the monitoring points increase and are highest nearest to the injection wellbore. After injection is discontinued, a gradual pressure decline is seen in the pressure history plots for all points. The declining pressure occurs due to dissipation of fluid pressure back to the wellbore and leak-off to neighboring fractures.   The initial pressure increase shows that the corresponding stress state around the neighboring well changes from the virgin in-situ stress condition and increases considerably, as shown in Figure 3.4. The perturbation of stress around the pressurized wellbore (points 5, 6, 7 and 8) as well as the neighboring non-active wellbore (points 1, 2, 3 and 4) is seen in the horizontal stress histories (Figure 3.5). Here the horizontal stress is seen to locally increase through a mechanism in which the vertical opening (dilation) of the horizontally generated hydraulic fracture results in the neighboring blocks (matrix) being compressed vertically, expanding horizontally, generating higher stresses horizontally since the blocks are confined. The stress perturbation decays as the pressure declines at these points, but still does not return to the virgin horizontal in-situ stress of 30 MPa. This means that if the second wellbore is pressurized, the stress field influencing the propagation of the second hydrofrac has been altered from that of the initial in-situ stress condition.  The stress shadow effect is manifest in the elevated principal stresses that develop around the pressurized wellbore and radiate outwards into the reservoir for a significant distance towards the second wellbore (Figure 3.6). The figure shows some small boundary effects on the right side of wellbore A, but the measurements are made on the left side of the   85 wellbore A therefore they are not significantly affected by that boundary effect. Also, the three different scenarios are equally affected by that boundary effect therefore the relative comparison made between different scenarios is likely affected non-significantly by the boundary effects. 3.4.1.2 Two alternating injection wells (zipperfrac) The second treatment scenario involves first injecting (i.e., fracturing) from wellbore B, and then shutting-in at wellbore B and injecting into wellbore A. The results are presented in Figures 3.7 and 3.8. Figure 3.7 shows the pressure increase at different points within the rock mass between the two wellbores. Points 1 to 4, located closer to the pressurized wellbore B, detect the pressure rise first with points 5 to 8 only detecting a muted response away from the injection. After shut-in a switch in from wellbore B to A, the monitoring points closer to wellbore A then respond and those closer to the first wellbore start showing pressure decay. The pressure decay around the first wellbore continues until the pressure front from the second wellbore arrives resulting in a lower rate of depressurization. The corresponding stress changes at the monitoring points are shown in Figure 3.8.  The perturbation of the in-situ stress field after the first and second stimulations can be seen in the maximum horizontal stress contours in Figure 3.9a and b, respectively. The stress perturbation decays after the injection pressure ceases, but still does not return to the lower virgin in-situ stress. This means that if the second wellbore is pressurized soon after the hydrofrac injection in the first wellbore, the hydrofrac generated will be controlled by a disturbed stress field and not the original in-situ stress field.    86 3.4.1.3 Two simultaneous injection wells (simulfrac) The simulation of hydraulic fracturing for the third scenario is conducted where both horizontal wellbores (A and B) are pressurized at the same time. The results are presented in Figures 3.10 and 3.11. Figure 3.10 shows the changes of fluid pressure at different points between the two wellbores. Here, simultaneous pressurization results in all six monitoring points showing a fast response to the hydrofrac injection.  After the injection is discontinued in both wellbores, a gradual pressure decline is seen in the pressure history plot for all points.  The corresponding stress responses around the two wellbores (Figure 3.11) are similar to those around wellbore A for the single well injection scenario (i.e. Figure 3.5), but with higher stress values. This indicates that the individual pressure fronts induced at the two neighboring wellbores communicate with each other and affect the stresses that develop at each adjacent wellbore. Further inspection of the modeled stresses through the horizontal stress contours shows the development of a stress shadow. Initially separate elevated stress zones develop, locally centered around each wellbore (Figure 3.12a). The elevated stress zones grow during the pressurization as show in Figures 3.12b and c and eventually, the two stress perturbations merge into a larger disturbed stress field (Figure 3.12d). Thus, the stress distribution around one wellbore in this case is altered by the other and vice versa. 3.4.2 Scenario comparison Figure 3.13 shows the change in induced hydraulic fracture aperture as a function of distance away from the injection well for the three different scenarios, where injection well is at zero.   87 Figures 3.14 and 3.15 compare the maximum hydraulic fracture apertures and lengths for each, respectively. These indicate that the hydraulic fracture aperture and length increase from conventional to simulfrac to zipperfrac.  The reason for these increases is the elevated fluid pressures resulting from a previous treatment which changes the effective stresses acting normal to the propagating hydraulic fracture. The method of interpretation here is based on superposition of the pressure effect of the neighboring well at the well in question (Matthews, 1961). Examining wellbore B for the conventional case, the pressure histories shown in Figure 3.4 suggest that the fluid pressures adjacent to the inactive wellbore (wellbore B) are affected by the neighboring active wellbore (wellbore A). The maximum formation pressure arising from the injection is approximately 19 MPa (relative to an initial background pressure of 10 MPa). The changes for the simulfrac scenario (shown in Figure 3.10) which indicates that the pressure fronts induced at both wellbores during simultaneous injection overlap, amplifying the respective pressures relative to those for the conventional case. The maximum formation pressure arising from the simultaneous injection is approximately 20 MPa. The pressure histories show an even larger increase for the zipperfrac scenario (Figure 3.7), specifically around wellbore A. Here the pressure histories show that the fluid pressures adjacent to wellbore A increase, despite it being inactive, due to the injection into the neighboring wellbore B. After wellbore B is shut-in and wellbore A becomes active, the pressure increase caused by pressurizing wellbore A is superimposed on the slowly diffusing pressure front already developed around wellbore A. This pressure response can be approximated by superimposing the individual pressure responses induced at each well in isolation (Matthews, 1961). The result is that this superposition of fluid pressures results in the highest formation fluid pressure between the   88 three scenarios; the maximum formation pressure arising from the alternating injection is approximately 22 MPa.   In both cases (zipperfrac and simulfrac), the higher fluid pressures arising from the stress interaction between wellbores results in larger fracture apertures, more concentrated fluid flow (via the hydro-mechanical coupling between aperture and flow), and decreasing effective stresses, which in turn produces longer hydraulic fracture lengths assuming equal injection times. These results are partly dependent on the model size, which is limited here by computing times required as well as the pumping rate and duration chosen for the analysis. For a typical hydrofrac stage, fluids may be pumped at a rate of 50 to 60 BPM (barrels per minute) for a duration of approximately 150 minutes (Vermylen and Zoback, 2011). For the analysis carried out here, numerical limitations restricted the fluid pump rate to 5 m3/min (30 BPM) for a duration of 90 minutes. Regardless, the results still show the relative hydraulic fracture length increase from conventional to simulfrac to zipperfrac scenarios.  3.5 Influence of in-situ stress and operational factors The influence of in-situ stress and operational factors on the stress shadow effect should be considered when designing a hydrofrac for maximum extent and recovery factor. In the following sections, the stress shadow effect is further examined as a function of reservoir depth (i.e. target formation depth), in-situ stress ratio, well spacing, and injection rate.      89 3.5.1 Influence of reservoir depth The influence of reservoir depth (i.e., pay zone depth) was examined whereby the vertical stress was increased while maintaining a constant in-situ stress ratio (K=1.4). Thus the sensitivity of the hydraulic fracture propagation to confining stress represented by the overburden depth was tested.  Figure 3.16 shows the hydraulic fracture length for the three different scenarios as a function of depth; each is normalized to the same conventional fracture length for all depths. Comparison of the fracture lengths as a function of both scenario and depth is not compatible as different fluid injection times were applied for the different modeled depths. This was necessary to enable the stress shadow effects from the neighboring wellbore to develop roughly equally for each depth, given that the higher in-situ stress conditions with depth adversely influenced the penetration rate of the injection. Thus, increased injection times were imposed for increasing reservoir depths to generate the same fracture length under the conventional scenario for all depths. The results confirm that the zipperfrac technique continuously produces the longest hydrofrac at each depth. More specific to the influence of reservoir depth, the results show that the relative difference between the conventional, simulfrac and zipperfrac fracture lengths increases as pay zone depth increases. This is because the perturbation to the in-situ stress field caused by the first hydrofrac in the zipperfrac sequence decays more slowly with increasing depth as the higher stress magnitudes dampen the rate of leak-off. As a result, the second wellbore is pressurized while there is still a significant stress shadow remaining from   90 the first wellbore hydrofrac, resulting in higher fluid pressures driving the hydrofrac propagation.   3.5.2 Influence of in-situ stress ratio Starting from the general stress state assumption that the horizontal to vertical in-situ stress ratio, K, is greater than one (required if assuming a horizontal hydrofrac to have the stress shadows, from the two neighboring wellbores, travel towards each other), the sensitivity of the modeled hydraulic fracture to K was tested for the three different stimulation scenarios (conventional, zipper and simul-frac). In these models, the vertical stress was held constant and the horizontal stress increased; often, the magnitude of the horizontal stress involves significant uncertainty and therefore added geo-risk to the success of the hydraulic fracture operation. The out of plane stress is assumed equal to the minimum stress and constant. Figure 3.17 shows the hydraulic fracture length for each scenario as a function of stress ratio, K, normalized to the longest fracture length for all model simulations. The figure shows hydraulic fracture length increases with increasing horizontal stress (as vertical stress is kept constant). Here, the higher horizontal stress facilitates fracture opening in the vertical direction, and as a result, the tensile rupture of the rock in the horizontal direction occurs with less resistance to dilation. An analogy can be drawn with Jaeger and Cook?s (1979) analytical solution for stress distribution around an elliptical opening of width W and height H subject to a vertical (p) and horizontal (Kp) stress, as shown in Figure 3.18.    91 The tangential stress at the tip of the elliptical opening (point A in Figure 3.18) can be calculated by:                                                  ???????? ???AAWKp ??21                                                       (3.2) where, the radius of curvature at the tip of the elliptical opening,A? , is equal toWH 2 . In the above equation, ????????AW?2 is a geometric term and depends on the geometry of the ellipse. With internal water pressure acting against p, the above equation demonstrates that the stress resisting fracture opening at point A decreases when K increases. Thus, the induced fracture length for a horizontal hydrofrac increases as the K ratio increases.   3.5.3 Influence of wellbore spacing The importance of wellbore spacing in ultimate economic recovery of reservoir has been studied by Holditch et al. (1978). They performed a series of numerical modeling to study the optimum well spacing for different reservoirs and concluded that determination of the optimum well spacing is not a ?common sense? type problem and different in-situ and operational factors must be considered before the optimum development plan can be formulated.  Simulation results in this study using different wellbore spacing suggest that, as would be expected, the distance between adjacent horizontal wells (i.e. stress shadow communication) affects the extension of the hydraulic fracture. Figure 3.19 shows the change in fracture aperture as a function of fracture length from the injection wellbore for different   92 wellbore spacing. These results are for the simulfrac scenario, which were chosen for this demonstration to omit the asymmetric influence of time on stress perturbation decay; i.e., with the two wellbores being pressurized at the same time, the respective leak-off of the fluid pressure front from each is symmetric. Figure 3.19 shows that the closer the neighboring wells, the stronger the stress shadow effect in terms of promoting longer induced hydrofracs. This arises because the pressure fronts reach each neighboring well faster and have less time to decay. A similar response is observed for the zipperfrac scenario, with increased fracture lengths developing for closer wellbores because the pressure around the first wellbore has less time to decay before injection is started from the second wellbore. The shorter the distance between the two wellbores, the more pressure from the earlier treatment is available to enhance any subsequent treatments. Figure 3.20 shows the fracture lengths for the different wellbore spacings, normalized to the longest hydrofrac modeled (i.e., zipperfrac with 60 m wellbore spacing) . The results show that although the closeness of neighboring wellbores has a strong positive inflence on both the zipper- and simul-fracs, the simulfrac scenario is more sensitive to wellbore spacing. This is because the interaction between the neighboring wellbore pressure fronts develops faster (relative to the zipperfrac) when they initiate at the same time. This is not as pronounced for the zipperfrac scenario because of the timing difference between wellbore pressurization schedules. For the zipperfrac, the fluid pressure from the first wellbore treatment diffuses towards the second wellbore but also leaks-off reducing the elevated pressure around that wellbore compared to the simulfrac case.    93 3.5.4 Influence of injection rate The influence of fluid injection rate on hydraulic fracturing from a single well has been investigated by de Pater and Beugelsdijk (2005), with focus on the characteristic time scales of the fracture process. They found from numerical modeling results that different fluid injection rates change the stress concentration around the wellbore, which in turn affect the hydraulic fracture propagation.  The fluid injection rate was also varied here between 5 m3/min and 10 m3/min, to model its influence for the three different scenarios. The results show that by increasing the injection rate, the stimulated fracture length increases for all three scenarios (Figure 3.21), and is most pronounced for the zipperfrac. This supports the findings presented earlier with respect to the sensitivity of the zipperfrac scenario to the pressure front. 3.6 Discussion and conclusions It is generally accepted that multiple hydraulic fracturing generated from neighboring horizontal wells provides the best completion option for stimulating unconventional and low permeability oil and gas reservoirs (e.g, Yost and Overbey, 1989). This approach has provided significant production results and is largely responsible for the successful development of several oil and gas fields in the United States and Canada, and is rapidly seeing acceptance worldwide. While the industry has made significant advancement in developing mechanical systems for carrying out these completions, progress in the design of multiple hydrofracs with respect to wellbore spacing, reservoir characterization, in-situ stress state, injection rate and overall production optimization requires a more thorough mechanistic understanding of their influence. Understanding the effect of stress shadowing   94 has significant benefits with respect to impact and risk mitigation on the cost and profitability of these operations.  This study has been carried out to contribute to the existing body of knowledge, specifically with respect to the modeling of stress shadow effects resulting from conventional, zipper- and simul-frac from neighboring wellbores. The results were obtained using distinct element numerical modeling techniques based on the hydro-mechanical coupled response of a fracture flow system. The results elaborate on the concept of stress perturbation through hydraulic fracturing, including the effect of multiple hydraulic fractures on stress shadowing as a function of in-situ stress and other operational factors, as a means to optimize shale completions by understanding the influence of these factors on hydraulic fracture propagation. The results simulating different fracturing techniques/schedules suggest that zipperfrac is the most effective technique where multiple horizontal wells are stimulated, producing the longest fractures for the same injection rate and volume compared to conventional and simulfrac techniques. In addition to adopting the most effective technique for fracturing multiple horizontal wells, improved understanding of the in-situ and operational conditions will allow designers and operators to help control the hydraulic fracture in the area being treated. In this study, the effects of in-situ stress as well as some operational factors are investigated to identify those factors that have the most influence on the effectiveness of the treatment. It was found that the influence of reservoir depth has the biggest influence on how much the zipperfrac outperforms the other two techniques. The horizontal to vertical stress ratio was also seen to have a considerable influence. Regarding the operational factors, wellbore spacing has a significant effect on the extension of hydraulic fractures promoting longer induced fractures   95 with closer wellbore spacings. The model responses were less sensitive to injection rate, but its influence was still significant. The simulations show great potential in providing a deeper understanding of the influence of stress shadow effects on the propagation of hydraulically induced fractures. Future work will include a three-dimensional simulation and use of microseismic data to help calibrate and ground-truth the model.                    96 Table 3.1. Strength and deformation properties assigned to the modeled planes of weakness. Incipient fracture property Value Units Friction angle 30 Degrees Residual  friction angle 25 Degrees Cohesion 1.0 MPa Residual cohesion 0.0 MPa Tensile strength 0.5 MPa Residual tensile strength 0.0 MPa Normal stiffness 1x104 MPa/m Shear stiffness 1x103 MPa/m             97   Figure 3.1. Discontinuum model: distinct element discretization.   Figure 3.2. Rock mass model assuming two orthogonal sets of incipient fracture planes (i.e., planes of weakness), parallel and perpendicular to bedding. Shown are the locations of the two horizontal wellbores and 8 monitoring points referred to in subsequent figures.    98  Figure 3.3. Three dimensional configuration of the two horizontal wellbores.   99  Figure 3.4. Conventional hydrofracing scenario (single well): Fluid pressure histories for different monitoring points between the two wellbores, with injection occurring from Wellbore A.    100   Figure 3.5. Conventional hydrofracing scenario (single well): Stress histories for different monitoring points between the two wellbores, with injection occurring from Wellbore A.    Figure 3.6. Stress distribution, resulting from the pressurization of wellbore A for a period of 90 minutes.   101  Figure 3.7. Zipperfrac scenario between two alternating injection boreholes: Fluid pressure histories for different points between the two wellbores. At 90 minutes, the injection at wellbore B is shut-in and injection started at wellbore A.   102   Figure 3.8. Zipperfrac scenario between two alternating injection boreholes: Stress histories for different points between the two wellbores. Point 1 is closest to the initial injection (wellbore B), with point 8 being closest to the subsequent injection (wellbore A).               103   Figure 3.9. Stress distributions arising from the zipperfrac scenario, after injection into: (a) wellbore B for a period of 90 minutes, and (b) then into wellbore A, for a period of 90 minutes.  (a) (b)   104  Figure 3.10. Simulfrac scenario between two simultaneously injecting hydrofrac wellbores: Fluid pressure histories for different points between the two wellbores. Point 1 is closest to injection well B and point 8 closest to injection well A.    105   Figure 3.11. Simulfrac scenario between two simultaneously injecting hydrofrac wellbores: Stress histories for different points between the two wellbores. Point 1 is closest to injection well B and point 8 closest to injection well A.                106     Figure 3.12 Stress distributions around simultaneously pressurized wellbores during pressurization (simulfrac scenario): (a) after 20 minutes of injection, (b) after 40 minutes of injection, (c) after 60 minutes of injection, (d) after 90 minutes of injection.  (a) (b) (c) (d)   107 0.10.20.30.40.50.6-30 -20 -10 0Induced Hydraulic Fracture Aperture (mm)Fracture Half Length (m)Zipperfrac Simulfrac ConventionalFigure 3.13. Hydraulic fracture comparison between the three modeled scenarios.  0.500 0.505 0.510 0.515 0.520 0.525 0.530 0.535Induced Fracture Aperture (mm)Zipperfrac Simulfrac ConventionalFigure 3.14. Maximum induced hydraulic fracture aperture comparison between the three modeled scenarios.   Injection well   108 Figure 3.15. Hydraulic fracture length comparison between the three scenarios. 0.94 0.96 0.98 1.00 1.02 1.04 1.06 1.08 1.10 1.12300025002000Normalized Fracture Length Reservoir Depth (m)Conventional Simulfrac ZipperfracFigure 3.16. Normalized hydraulic fracture length comparison between the three scenarios for different reservoir depths.   109 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.001.51.41.31.2Normalized Fracture LengthStress RatioConventional Simulfrac ZipperfracFigure 3.17. Normalized hydraulic fracture length comparisons between different completion scenarios for different horizontal to vertical in-situ stress ratios.   Figure 3.18. Geometry of an elliptical opening in a biaxial stress field.   110 0.10.20.30.40.50.6-30 -20 -10 0Induced Hydraulic Fracture Aperture (mm)Fracture Half Length (m)100908070Figure 3.19. Hydraulic fracture comparison for different wellbore spacing.  0.78 0.80 0.82 0.84 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.0010090807060Normalized Fracture LengthWellbore Spacing (m)Zipperfrac SimulfracFigure 3.20. Normalized hydraulic fracture length for simulfrac and zipperfrac scenarios with different wellbore spacing.  Wellbore  Spacing (m)  Injection well   111 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00ConventionalSimulfracZipperfracNormalized Fracture LengthCompletion ScenariosHigh Injection Rate Low Injection RateFigure 3.21. Normalized hydraulic fracture length for simulfrac and zipperfrac scenarios with different injection rate.                112 Chapter  4: A numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing 4.1 Introduction Hydraulic fracture operations involve the injection of fluids from a wellbore into a formation to maximize the extraction of oil and gas resources from the reservoir. These stimulations allow the development of previously uneconomical reserves, for example shales gas. However, it has been well established that these fluid injections also induce seismicity (e.g., Hollister and Weimer (1968), Ohtake (1974), Fletcher and Sykes (1977), Pearson, (1981), Talwani and Acree, (1985), Simpson et al. (1988), Zoback and Harjes (1997)).  The occurrence of induced seismicity, whether through hydraulic fracturing, enhanced geothermal systems, or carbon dioxide sequestration, represents an important consideration in the risk profile of an injection well. Zoback and Gorelick (2012) suggested that because earth crust is naturally critically stressed, by injecting fluid in deep wells the pore pressure increases and causes decrease in effective stress in the vicinity of the pre-existing faults. When the pore pressure increases, the shear resistance decreases which causes fault to slip. This allows the release of elastic energy, which is already stored in the surrounding rocks (Healy et al., 1968).  This paper describes a numerical experiment investigating the influence of pore pressure diffusion on the modeling of seismicity after a hydraulic fracture treatment. One objective is to evaluate the utility of 2-D distinct element techniques to model both pore pressure build-up and diffusion along existing rock weakness planes in response to fluid   113 injection and subsequent changes to the effective stress field. A second objective is to estimate the corresponding maximum seismic magnitude that is likely to occur in response to the hydraulic fracture treatment. 4.2 Effect of fluid injection on induced seismicity The occurrence of injection-induced seismicity is usually confined in both space and time, with pressure build-up and diffusion controlling the spatial and temporal pattern of the seismicity after a hydrofrac treatment. Spatially, the problem can be conceptualized as a scenario in which a wellbore in the vicinity of a critically-stressed fault is pressurized and the injected fluids diffuse towards the fault, increasing pore pressures and reducing the effective normal stress until slip is triggered along a portion of the fault. The transient nature of the fluid front as it radiates outward from the injection well allows for the triggering of events after injection and shut-in, referred to here as post-injection seismicity.   The importance of fluid pressure in generating induced seismicity was demonstrated in 1962 in Denver, Colorado, when it was observed that a series of seismic events were being generated shortly after a chemical weapons plant had begun injecting contaminated wastewater down a deep well. Detailed studies by Hollister and Weimer (1968) and Healy et al. (1968) confirmed the causal observations of Evans (1966) that the seismic events were induced by the waste fluid injection. Similar relations between local seismicity and fluid injection were likewise observed at the Rangely Oil Field, Colorado by Pakiser et al. (1969), Healy et al. (1972) and Gibbs et al. (1973), and in the Los Angeles basin by Teng et al. (1973).    114 More recently, post-injection seismicity associated with gas shale and geothermal projects have entered the public spotlight with heightened sensitivity directed towards hydraulic fracturing practices and induced seismicity. Between 2000 and 2005, numerous seismic events were recorded at the European geothermal project in Soultz-sous-For?ts, France, with events as large as M=2.5 being recorded during injection and M=2.6 several days after shut-in (Charlety et al., 2007). For stimulations carried out in 2000, 2003 and 2004 the largest seismic events occurred just after or after shut-in. Similar experiences were encountered at Basel, Switzerland where induced seismicity reached a Richter magnitude of 3.4 despite precautionary reductions of the injection rate, leading to suspension of its hot dry rock enhanced geothermal systems project (Haring et al., 2008).  Experiences with hydraulic fracturing treatments in gas shales have generated similar seismicity, generally viewed as low risk due to the remote nature of the locations. In the Bowland shale, Lancashire, UK, two notable events were recorded with Richter magnitude of 1.5 and 2.3. These occurred almost 10 hours after shut-in (de Pater and Baisch, 2011). de Pater and Baisch (2011) suggested that reducing the treatment volume is one means to mitigate induced seismicity. In a detailed study in the Horn River Basin in northeastern British Columbia, Canada, approximately 38 events ranging between Richter magnitudes 2.2 and 3.8 were recorded between 2009 and 201l, with only one being felt on surface in this remote region (BCOG, 2012). In their 2012 report of the findings from this investigation, the B.C. Oil and Gas Commission found that several of the Horn River events greater than M 2.0 were located along faults intersecting the wellbores. They also note that there were other instances where faults intersected wellbores without anomalous events being detected (BCOG, 2012).   115 4.3 Numerical modeling methodology The Distinct Element Method (DEM) applies a Lagrangian formulation to compute the motion and interaction between a series of discrete deformable blocks, representing the problem domain, via compliant contacts and Newton?s equation of motion (Cundall and Hart, 1993). This enables the problem domain to be divided through by one or more discontinuity sets of variable orientation, spacing and persistence. One fundamental advantage of the DEM is that pre-existing (natural) joints in the rock mass can be modeled explicitly and allow for joints to undergo large deformations in shear (slip) or opening (dilation). The 2-D commercial code UDEC (Universal Distinct Element Code; Itasca Consulting Group, 1999) is used here to model the response of a jointed rock mass subjected to static loading and hydraulic injection. UDEC is capable of modeling the behavior of weak jointed rock masses in which both the deformation and yielding of weak rock and slip along pre-existing discontinuities are important controlling factors. Progressive failure associated with crack propagation and fault slip can be simulated by the breaking of pre-existing contacts between the pre-defined joint bounded blocks, which although deformable, remain intact.  Key for simulating hydrofracturing, UDEC has the capability to model fluid flow through the defined fracture network. A fully coupled hydro-mechanical analysis can be performed, in which fracture conductivity is dependent on mechanical deformation of joint apertures and, conversely, joint water pressures can affect the mechanical computations of joint aperture. The blocks in this assemblage are treated as being impermeable, and fracture flow is calculated using a cubic law relationship for joint aperture:   116                                                                  lPkaq ?? 3                                                         (4.1)  where, k is a joint conductivity factor (dependent on the fluid dynamic viscosity), a is the contact hydraulic aperture, ?P is the pressure difference between the two adjacent domains, and l is the length assigned to the contact between the domains. Since the UDEC formulation is restricted to the modeling of fracture flow, it should be noted that leak-off along the fractures diffusing into the rock matrix is assumed to be zero (only leak-off into the incipient fractures is considered). Furthermore, the cubic law flow assumption limits tortuosity. When a joint contact is broken, the fluid flows into the joint.  4.4 Simulation setup and parameters The rock mass modeled in this study is represented by two persistent orthogonal planes of weakness (Figure 4.1). These serve as incipient planes along which hydrofrac propagation is restricted. The simulation of induced seismicity is executed through the inclusion of a fault, which extends across the model. An injection well is located such that a significant portion of the fluid injected diffuses towards the fault and eventually penetrates the fault following shut-in. The fluid pressure decays slowly after the injection and the disturbed pressure front diffuses through the surrounding rock mass. Although the fluid pressures decrease with time and distance, there is still sufficient pressure to trigger fault slip. The fault model is based on interactions between neighboring fault segments allowing the model to simulate slippage on a single contact together with the subsequent interactions and responses of its neighboring contacts.    117 The input material parameters include both those for the rock matrix and incipient planes of weakness and fault. The rock matrix was modeled as being elastic, assuming typical values for shale (Young?s modulus=30 GPa, Poisson?s ratio=0.25). The incipient planes of weakness and fault were modeled assuming a Coulomb-slip constitutive model with both peak and post-peak properties. These are given in Table 1. The fault plane represented by the model has 1000 m depth. The maximum and minimum stresses were assumed to be 30 and 20 MPa, respectively. Boundary conditions are specified for the external boundaries of the model assuming constant stress and constant fluid pressure conditions. The time required for executing the model in this study was typically on the order of 8 days. 4.5 Post-injection seismicity simulations Numerical simulations were performed using the model developed to investigate the influence of hydraulic fracturing on fluid pressure changes around the neighboring fault and any subsequent shear slippage along the fault. Both fluid pressure build-up during the injection until the time of shut-in as well as fluid pressure diffusion after the shut-in were considered. 4.6 Fluid pressure build-up during injection and diffusion after shut-in Experience gained from mapping hundreds of hydraulic fracturing treatments with downhole geophones has shown that the occurrence of seismic event induced during a treatment is greatly influenced by the injection volume and rate used. Here, the hydraulic fracturing simulation was conducted by pressurizing the wellbore in the vicinity of a critically stressed fault (Figure 4.1). Figure 4.2 shows the joint fluid pressure distribution in the rock mass at   118 the time that fluid injection is stopped (i.e., shut-in). It can be seen that the fluid pressure has not reached the fault at the end of the hydraulic fracture treatment and shut-in. Here the fluid pressure treatment was applied for a period of 50 minutes. After shut-in, the injected fluids continue to diffuse and radiate towards the fault, eventually penetrating it (Figure 4.3). The response of the fault to the elevated fluid pressures then depends on the spatial and temporal characteristics of the diffusion pulse, with a series of slip events being produced as opposed to a single event. The strongest slip event was typically observed late in the sequence, sometimes long after injection had stopped (up to 150 minutes for the model simulations performed here). 4.7 Shear slip of critically stressed fault The shear slip distribution along the fault, after 50 minutes of injection and 150 additional minutes of shut-in, is presented in Figures 4.4. The figure shows a distribution of slip magnitudes along the length of the fault, with a maximum fault slip of approximately 10 mm. Displacements between 5 and 10 mm were observed along 680 meters of the total 800 meter fault length, with an average fault slip of 8 mm. This average slip magnitude was then used to calculate the maximum earthquake magnitude, as presented in the following section. 4.8 Seismic moment and moment magnitude calculations The seismic moment is a direct measurement of the energy released by a seismic event and is related to the strength characteristics of the fault. It can be calculated as follows (Kramer, 1996):      DAMo ??                                                                         (4.2)   119 where: oM is the seismic moment (dyne?cm),  ? is the rupture strength of the rock along the fault (dyne/cm2),  A is the rupture area (cm2), and  D  is the average amount of slip (cm). Here, the rupture strength of the fault is equal to 1E10 dyne/cm2, the length of the fault that slipped is equal to 6.8E4 cm (680 m), and the average amount of slip is equal to 0.8 cm. Assuming a unit depth of fault slip due to the 2-D nature of the model, for a fault with 1000 m depth, the rupture area is equal to 6.8E9 cm2 and the resulting seismic moment is equal to 5.4E19 dyne.cm Seismic moment can be converted into a moment magnitude using the following relationship (Hanks and Kanomori, 1979):                               7.10log3/2 ?? ow MM                                                              (4.3) where: wM  is moment magnitude (dyne.cm), and    oM  is seismic moment (dyne.cm). The moment magnitude calculated from equation 3 is 2.45. Figure 4.5 is used to convert the calculated moment magnitude to the Richter local magnitude. As seen in Figure 4.5, for earthquake magnitudes smaller than 6, the moment and local magnitudes are near equal; the local magnitude for the calculated moment magnitude is equal to 2.45.   120 Figure 4.6 uses a well-established seismological relationship that correlates earthquake magnitude to the size of the fault that slipped and seismic moment (Stein and Wysession, 2002). This suggests that only faults that are at least tens of kilometers long are capable of producing large seismic events with magnitudes exceeding 6 (Zoback and Gorelick, 2012). Zoback and Gorelick (2012) note here that the fault size in Figure 4.6 is a lower bound value that refers to the size of the fault segment that slips in a given earthquake. As a geological feature, the total fault length is generally much longer than the part (segment) that slips during an individual event. As shown in Figure 4.6, an active fault slip segment of 680 meters, as produced in the UDEC model, is capable of producing an earthquake with a magnitude between 2.3 and 3.8 for a fault slip displacement between 1mm and 1 cm, depending on the magnitude of stress released by the event, or stress drop. For the case modeled here, the 8 mm of average slip produced corresponds to an event with a magnitude equal to 2.45, which is in the range shown in Figure 4.6. 4.9 Discussion and conclusions A numerical modeling study has been carried out to investigate the application of the distinct element technique to the simulation of fault slip and induced seismicity resulting from a hydraulic fracture treatment via a wellbore in the vicinity of a natural fault. The model was able to predict the occurrence of post-injection seismicity in response to diffusion of the injected fluids, in a system governed by fracture permeability, long after the hydraulic fracturing treatment is finished and shut-in is initiated.  In most areas, regional-scale faults should be easily identified during geological site characterization studies. Smaller-scale faults or those that are shallow dipping (i.e., that do   121 not daylight on surface and therefore would not be detected through surface mapping) may be more difficult to locate a priori. Experience gained from monitoring hundreds of hydraulic fracturing treatments with downhole geophones has shown that the occurrence of seismic events induced by the treatments is greatly influenced by the injection volume used during the operation. Improved understanding of this condition will allow designers and operators to control the amount of injection during a hydraulic fracturing treatment to minimize fluid pressure diffusion and subsequent slip of a neighboring fault. The simulations presented here show great potential in providing a deeper understanding of the effect of natural fracture systems and pressure diffusion of fracing fluids into a neighboring fault causing shear slip and induced seismicity after a hydraulic fracture treatment. Future work will include modeling of hydrofrac fluid diffusion into a fault in the vicinity of the hydrofrac treatment site in a three-dimensional model.            122 Table 4.1. Properties assigned to the modeled planes of weakness and fault. Discontinuity property Incipient fractures Fault Units Friction angle 30 20 Degrees Residual  friction angle 25 20 Degrees Cohesion 1.0 0.0 MPa Residual cohesion 0.0 0.0 MPa Tensile strength 0.5 0.0 MPa Residual tensile strength 0.0 0.0 MPa Normal stiffness 1e4 10 MPa/m Shear stiffness 1e3 1 MPa/m                     123  Figure 4.1. Rock mass model with two sets of weakness planes and a fault.   Figure 4.2. Pore pressure distribution at the time of shut-in.   124   Figure 4.3. Pore pressure distribution 2 ? hours after shut-in.    Figure 4.4. Fault shear displacement 2 ? hours after shut-in.    125  Figure 4.5. Saturation of various magnitude scales: Mw (moment magnitude), ML (Richter local magnitude), Ms (surface wave magnitude), mb (short-period body wave magnitude), mB (long-period body wave magnitude), and MJMA (Japanese Metrological Agency magnitude). After Idriss (1985) and Kramer (1996).                         126   Figure 4.6. Relationship among various scaling parameters for earthquakes. The larger the earthquake, the larger the fault and amount of slip, depending on the stress drop in a particular earthquake. Observational data indicate that earthquake stress drops range  between 0.1 and 10 MPa (modified after Zoback and Gorelick, 2012).             127 Chapter  5: Thesis Discussion and Conclusions This thesis presents the results from a series of detailed numerical investigations examining the behavior of hydraulic fracturing in a discontinuous rock mass using 2-D distinct element techniques. First, findings from numerical simulations of hydraulic fracture initiation and propagation within a naturally fractured rock mass were presented to demonstrate the influence that pre-existing fractures may potentially have on hydrofrac propagation. Results were presented for different hydrofrac approach angles to the natural fractures under different in-situ conditions and showed that the presence of natural fractures may cause a hydrofrac to arrest or offset. Next, results were presented examining the influence of stress shadow effects on the extent of hydraulic fracture for reservoir treatments involving multiple hydraulic fractures from adjacent horizontal wells. Comparisons were made for models simulating conventional, simulfrac and zipperfrac treatment scenarios. Sensitivity analyses were performed for design variables such as wellbore spacing, in-situ stress state, reservoir depth and injection rate for overall production optimization. Lastly, the effect of hydraulic fracturing and pressure diffusion on triggering slip along a neighboring fault was numerically investigated. The results demonstrated the influence of fluid diffusion generated after shut-in on the activation of a critically stressed fault.  Together, the results from this research help to further: 1) the evaluation of 2-D distinct element techniques that account for both shear and dilation along existing joints and tensile rupture of intact rock in response to fluid injection, and subsequent changes to the effective stress field for modeling hydraulic fracturing; 2) understanding of the influence of two important factors, natural fractures and stress shadow effects, on hydraulic fracture   128 propagation;  3) the assessment of geo-risk arising from the spatial and temporal interaction between a hydraulic fracture and a critically stressed fault, resulting in induced seismicity after shut-in.  5.1 Chapter summaries 5.1.1 Influence of natural fractures and in-situ stresses on hydrofrac propagation A comprehensive numerical study was performed to investigate the influence of natural fractures on hydrofrac propagation, aided by a Voronoi tessellation procedure to model a background of incipient fractures superimposed with a network of fluid conducting joints. Emphasis was placed on the examination of the behavior of a hydraulic fracture approaching a natural fracture under different angles of approach relative to different in-situ stress conditions. The influence of the proximity of the natural fracture to the wellbore during the initial stages of hydraulic fracture propagation was also investigated.   The results from these numerical simulations demonstrated that the interaction between the natural fractures and hydraulically induced fracture may result in varied responses. These include those that are neutral in terms of geo-risk, with the hydrofrac crossing or offsetting at its intersection with the natural fracture, but still serving to connect and enhance permeability to the wellbore. Alternatively, other responses were seen to be more adverse with respect to geo-risk, with the hydrofrac arresting at its intersection with the natural fracture.  The implication of this study is that the relationship between the orientation and magnitudes of the in-situ principal stresses relative to the orientation of the natural fracture sets present in the rock play an important role on the extent, symmetry and effectiveness of the hydraulic fracture propagation with respect it?s reservoir coverage.   129 Predicting the interaction behavior between a hydraulic fracture and the natural fracture system plays an important role in planning and executing a successful hydraulic fracture stimulation. 5.1.2 Influence of stress shadow on fracture propagation in multi-horizontal wells A detailed and comprehensive numerical study was performed to investigate the effects of stress shadowing on the propagation of multiple hydraulic fractures from adjacent horizontal wells. These stress shadow effects are potentially critical to the design of multiple hydraulic fracturing treatments. Understanding the effect of stress shadowing has significant benefits with respect to quantifying and mitigating the impact and risk related to cost and effectiveness of these operations.  The analysis is carried out for three different completion scenarios (conventional, simulfrac, and zipperfrac) to investigate their effect on the propagation of hydraulic fractures. The results elaborate on the concept of stress perturbation through hydraulic fracturing, including the effect of multiple hydraulic fractures on stress shadowing as a function of in-situ stress and other operational factors, as a means to optimize shale completions by understanding the influence of these factors on hydraulic fracture propagation.      5.1.3 Influence of hydraulic fracturing on fault slip and post-injection seismicity Results were reported from 2-D distinct element techniques that modeled both pore pressure build-up and diffusion through leak-off into a system of interconnected fractures. The model response to fluid injection and the subsequent changes to the effective stress field were of particular interest. The model was able to predict the occurrence of post-injection seismicity   130 in response to diffusion of the injected fluids, in a system governed by fracture permeability, long after the hydraulic fracturing treatment was finished and shut-in was initiated. Using the numerical outputs, the corresponding maximum seismic magnitude that is likely to occur in response to the hydraulic fracture treatment was estimated. The simulations show great potential in providing a deeper understanding of the effect of natural fracture systems and pressure diffusion of fracing fluids into a neighboring fault causing shear slip and induced seismicity after a hydraulic fracture treatment. 5.2 Key conclusions and scientific contributions The key scientific contributions arising from this study are as follows: 1. Hydraulic fracturing theory is largely founded on linear elastic fracture mechanics, with many standard hydraulic fracturing simulators adopting assumptions such as the generation of bi-wing, perfectly planar, symmetrical fractures. Whether empirical or numerical, these methods inherently treat the rock mass as an elastic continuum; the resulting hydraulic fractures modeled/predicted are symmetric, extending in equal distances in opposite directions from the wellbore. The results presented here adopt an alternative approach based on discontinuum mechanics that treats the geology in a more realistic manner with respect to accounting for both the presence of natural fractures and also incipient fractures that represent planes of weakness and potential fracture pathways. The use of distinct element model applied here clearly demonstrated the value of accounting for the presence of natural discontinuities with respect to designing a hydraulic fracture treatment, identifying that the natural fracture system may result in significant asymmetry in the final hydrofrac. The   131 impact of an arrested hydrofrac or unequal radial distribution of fluid pressure around the wellbore can lead to significant underperformance compared to that which would be predicted by treating the rock mass as a continuum.  2. The application of advanced numerical modeling in this study demonstrated the value of such tools in helping to investigate the possible range of reservoir rock behavior during hydraulic fracturing where data/experience may be limited or biased, or where complex interactions arise, for example in the case of multi-frac treatments in the same volume of rock. Results provided here examining the influence of stress shadows on hydraulic fracture propagation, where limited experience exists, and the subsequent sensitivity analyses, represent value-added contributions to the understanding of disturbed stress fields with respect to their impact on the cost and effectiveness of these operations.   3. The selected discontinuum modeling technique is used to further demonstrate the utility of the method for modeling fault slip triggered by hydrofracing activities. The model was able to predict the occurrence of post-injection seismicity in response to diffusion of the injected fluids, in a system governed by fracture permeability, long after the hydraulic fracturing treatment was finished and fluid pressures shut-in. Using the numerical model outputs, the corresponding maximum seismic magnitude was estimated. Improved understanding of this condition will allow designers and operators to test different diffusion responses and the potential for fault slip after shut-in for different injection pressures and schedules. Improved understanding of this condition will allow designers and operators to control the amount of injection   132 during a hydraulic fracturing treatment to minimize fluid pressure diffusion and subsequent slip of a neighboring fault.   5.3 Further research 1. This study highlighted the practical value of 2-D discontinuum analyses, using the commercial code UDEC, for performing numerical modeling of hydraulic fracturing by accounting for the presence of both natural and insipient fractures in the reservoir rock. Inherent in the 2-D analysis were simplifications regarding the out-of-plane conditions. It should be noted that the 2-D models executed in this study required considerable running times (on the order of 7 to 10 days). This restricted the use of 3-D modeling using the same distinct-element formulation as this would have required significantly longer run times. However, with continuous advances in computation speed and pre-and post-processing software capabilities, future research should explore the application of 3-D discontinuum codes like 3DEC. 2. Another state-of-practice may involve use of 2-D and 3-D hybrid FEM-DEM brittle fracture analyses using ELFEN to incorporate the role of intact rock progressive failure associated with hydraulic fracturing. However, the use of these codes is limited at present by the considerably long computing run times required by these codes. This limitation will also be overcome by continuous advances in computation speed and pre-and post-processing software capabilities.  3. The modeling performed in this study demonstrated the predicted behavior of hydraulic fracturing based on model dimensions that pushed the limits of acceptable run times. However, the dimensions of the selected models were smaller than those   133 encountered in actual hydraulic fracture treatments. Both computing memory requirements limiting the number of blocks and elements that could be used together with required run times (several days to weeks), forced restrictions on the size of the models solved. For the same reason the injection rate and injection times selected were smaller than those typically used during a hydrofrac treatment. These limitations will also be avoidable in future as software and computing capabilities are advanced. 4. The model scenarios presented were largely based on conceptual realizations. Further work should be carried out using actual site specific data to properly constrain and calibrate the models. The use of microseismic data, for example, represents one avenue for further investigation of hydrofrac behavior and model calibration, by simulating synthetic seismicity to compare to specific hydrofrac performance and rock mass response in an actual gas shale reservoir.              134 Bibliography  Adachi, J. I., Siebrits, E., Peirce, A., and Desroches, J., 2006. Computer simulation of hydraulic fractures. International Journal of Rock Mechanics and Mining Sciences. 44: 739-757. Aggour, T. M., and Economides, M. J., 1998. Optimization of the performance of high-permeability fractured wells. Paper presented at Proceedings of the 1998 International Symposium on Formation Damage Control, February 18-19. Al-Busaidi, A. 2005. Distinct element modeling of hydraulically fractured Lac du Bonnet granite. Jounal of Geophysical Research, 110: B06302. Balen, R. M., Meng, H. Z., and Economides, M. J., 1988. Applications of the net present value (NPV) in the optimization of hydraulic fractures. Paper presented at Proceedings: 1988 Eastern Regional Conference, November 1. Barree, R. D., 1983. A particular numerical simulator for three-dimensional fracture propagation in heterogeneous media. Society of Petroleum Engineers. SPE 12273, 1983. BCOG, 2012. Investigation of observed seismicity in the Horn River Basin. British Columbia Oil and Gas Commission Report 2012. Blanton, T. L., 1982. An experimental study of interaction between hydraulically induced and pre-existing fractures. Society of Petroleum Engineers. SPE/DOE 10847.   135 Blanton, T. L., 1986. Propagation of hydraulically and dynamically induced fractures in naturally fractured reservoirs. Society of Petroleum Engineers. SPE 15261. Bobet, A., Fakhimi, A., Johnson, S., Morris, J., Tonon, F., and Yeung, R., 2009. Numerical Models in Discontinuous Media: Review of Advances for Rock Mechanics Applications. Geoenvironmantal Engineering. 135(11): 1547-1561. Carter, B. J., Desroches, J., Ingraffea, A. R., and Wawrzynek, P. A., 2000. Simulating fully 3D hydraulic fracturing. Modeling in Geomechanics. Ed. Zaman, Booker, and Gioda, Wiley Publishers, 730.  Charlety, J., Cuenot, N., Dorbath, L., Dorbath, C., Haessler, H., and Frogneux, M., 2007. Large earthquakes during hydraulic stimulations at the geothermal site of Soultz-sous-Forets. International Journal of Rock Mechanics and Mining Sciences. 44: 1091-1105. Choi, S. O., 2003. Numerical study on the estimation of the shut-in pressure in hydraulic fracturing test. Geosystem Engineering. 55-62. Chopra, S. C., 2008. Applied Numerical Methods with MATLAB. Mc Graw Hill, 2nd  edition, 588 pp. Cliff Minerals Inc., 1982. Analysis of the Devonian shale in the Appalachian basin: Final report to US Department of Energy under contract DE-AC21-80MC14693, V.1 and 2, Morgantown WV, 56p. Cramer, D. D., 2008. Stimulating unconventional reservoirs: Lessons learned, successful practices, areas for improvement. Paper presented at Unconventional Reservoirs Conference, February 10-12, 2008.   136 Cundall, P. A. and Hart, R. D., 1992. Numerical modeling of discontinua. Engineering Computations. 9:101-113. Cundall, P. A., and Hart, R. D., 1993. Numerical modeling of discontinua. Comprehensive Rock Engineering. 2: 231-243. A. Hudson, ed. Oxford: Pergamon Press Ltd. Damjanac, B., Gil, I., Pierce, M., and Sanchez M., 2010. A new approach to hydraulic fracturing modeling in naturally fractured reservoirs. 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium, June 27-30, 2010, Salt Lake City, Utah. Daneshy, A., Au-Yeung, J., Thompson, T., and Tymko, D., 2012. Fracture shadowing: A direct method for determining of the reach and propagation pattern of hydraulic fractures in horizontal wells. Society of Petroleum Engineers. SPE 151980. de Pater, C. J. and Beugelsdijk, L. J. L., 2005. Experiments and numerical simulation of hydraulic fracturing in naturally fractured rock. Paper ARMA/USRMS 05-780, Alaska Rocks 2005, 40th US RockMechanics Symposium on Rock Mechanics, Anchorage, Alaska, USA, 25-29 June. de Pater, C. J., and Baisch, S., 2011. Geomechanical study of Bowland Shale seismicity, Synthesis Report 2011. Dershowitz, W. S., Cottrell, M. G., Lim, D. H., and, Doe, T. W., 2010. A discrete fracture network approach for evaluation of hydraulic fracture stimulation of naturally fractured reservoirs. 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium. June 27 - 30, 2010, Salt Lake City, Utah.   137 Detournay, E., and Cheng, A. H. D., 1991. Plane strain analysis of a stationary hydraulic fracture in a poroelastic medium. International Journal of Solids and Structures. 27(13):1645-1662. Dong, C. Y., and de Pater, C. J., 2002. Numerical modeling of crack orientation and link-up. Advances in Engineering Software. 33, 577-587. Dusseault, M., and McLennan, J., 2011. Massive multistage hydraulic fracturing. 45th US Rock Mechanics/Geomechanics Symposium, San Francisco, USA. Economides, M. J., Hill, A. D. and Ehlig-Economides, C. A., 1994. Petrolum production systems, Englewood Cliffs, New Jersey, USA, Prentice-Hall. Economides, M. J., and Nolte, K. G., 2000. Reservoir Stimulation, 3rd edition. John Wiley & Sons, Ltd., New York. Economides, M. J., Oligney, R. and Valko, P., 2002. Unified Fracture Design-Bridging the Gap between Theory and Practice, Orsa Press, Alvin, TX. Economides, M. J., and Martin, T., 2007. Modern Fracturing, ET Publishing, Houston, TX. Economides, M. J., and Demarchos, A., 2008. Benefits of a p-3D over a 2D model for unified fracture design. Paper presented at SPE International Symposium and Exhibition on Formation Damage Control 2008, Febrary 13-15, 2008. Evans, D. M., 1966. The Denver area earthquakes and the Rocky Mountain arsenal disposal well, Mountain Geologist. 3: 23-36.   138 Fisher, M. K., Heinze, J. R., Harris, C. D., Davidson, B. M., Wright, C. A., and Dunn, K. P., 2004. Optimizing horizontal completion techniques in the Barnett Shale using microseismic fracture mapping. Society of Petroleum Engineers, SPE 90051. Fletcher, J. B. and Sykes, L. R., 1977. Earthquakes related to hydraulic mining and natural seismic activity in western New York state. Journal of Geophysical Research. 82(26): 13767-3780. Geertsma, J. and de Klerk, F., 1969. A rapid method of predicting width and extent of hydraulic induced fractures,? paper SPE 2458, Journal of Petroleum Technology. 21, 1571-1581. Gibbs, J. F., Healy, J. H., Raleigh, C. B., and Coakley, J., 1973. Seismicity in the Rangely, Colorado Area: 1962-1970. Bulletin of Seismological Society of America. 63: 1557-1570. Gil, I., Damjanac, B., and Nagel, N. 2010.  Geomechanical evaluation of solids injection. 44th U.S. Rock Mechanics Symposium and 5th U.S.-Canada Rock Mechanics Symposium, June 27 - 30, 2010, Salt Lake City, Utah. Goodman, R. E., 1976. Methods of Geological Engineering in Discontinuous Rocks, West publishing, St. Paul. Haimson, B. C., and Cornet, F. H., 2003. ISRM suggested methods for rock stress estimation; part 3, hydraulic fracturing (HF) and/or hydraulic testing of pre-existing fractures (HTPF); rock stress estimation. International Journal of Rock Mechanics and Mining Sciences. 40, (7-8): 1011-1020.   139 Hamidi F., and Mortazavi A., 2012. Three dimensional modeling of hydraulic fracturing process in oil reservoirs. 46th U.S. Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association (ARMA), Chicago, United States. Hanks, T.C., and Kanamori, H., 1979. A moment magnitude scale. Journal of Geophysical Research. 84: 2348-2350. Haring, M. O., Schanz, U., Ladner, F., and Dyer, B. C., 2008. Characterization of the Basel 1 enhanced geothermal system. Geothermics. 31: 469-495. Healy, J. H., Rubey, W. W., Griggs, D. T., and Raleigh, C. B., 1968. The Denver earthquakes, Science. 161: 1301-1310. Healy, J. H., Lee, W. H. K., Pakiser, L. C., Raleigh, C. B., and Wood, M. D., 1972. Prospects for earthquake prediction and control, Tectonophysics. 14: 319-332. Hoek, E., Grabinsky, M. W., and Diederchs, M. S., 1991. Numerical modeling for underground excavation design. Transactions of the Institution of Mining and Metallurgy (Section A: Mining Industry), 100, A22?A30. Holditch, S. A., Jennings, J. W., and Neuse, S. H., 1978. The optimization of well spacing and fracture length in low permeability gas reservoirs. Society of Petroleum Engineers. SPE 7496. Hollister, J. C. and Weimer, R. J., 1968. Geophysical and geological studies of the relationship between the Denver earthquakes and the Rocky Mountain arsenal well. Quaterly, Colorado School of Mines. 63(1): 1-251.   140 Hubbert, M. K., and Willis, D. G. 1957. Mechanics of hydraulic fracturing. Journal of Petroleum Technology. 9(6): 153-66. Idriss, I. M., 1985. Evaluating seismic risk in engineering practice. Conference Proceeding, 11th International Conference on Soil Mechanics and Foundation Engineering. 1:255-320. Itasca Consulting Group Inc., 1999. UDEC/3DEC (Universal Distinct Element Code in 2/3 Dimensions), Version 5.0. Minneapolis, MN: ICG. Itasca Consulting Group Inc., 1999. PFC2D/3D (Particle Flow Code in 2/3 Dimensions), Version 2.0. Minneapolis, MN: ICG. Ito, T. and Hayashi, K., 1993. Analysis of crack reopening behavior for hydrofac stress measurement. International journal of rock mechanics and mining sciences & geomechanics. 30: 1235-1240. Jaeger, J. C., and Cook, N. G. W., 1979. Fundamentals of Rock Mechanics, 3rd ed., Chapman & Hall, London. Kennedy, R. L., Gupta, R., Kotov, S., Burton, A., Knecht, W. N., and Ahmed, U., 2012. Optimized shale resource development: Proper placement of wells and hydraulic fracture stages. Society of Petroleum Engineers. SPE 162534. Kramer S. L., 1996. Geotechnical Earthquake Engineering. Prentice-Hall, Englewood Cliffs, N. J., 653 pp.   141 Marongiu-Porcu, M., Economides, M. J., and Holditch, S. A., 2008. Economic and physical optimization of hydraulic fracturing. Paper presented at SPE International Symposium and Exhibition on Formation Damage Control, Febrary 13-15. Matthews, C. S., 1961. Analysis of pressure build-up and flow test data. Journal of Petroleum Technology. 13 (9): 862-870. McLennan, J., Tran, D., Zhao, N., Thakur, S., Deo, M., Gil, I., Damjanac, B., 2010. Modeling of fluid invasion and hydraulic fracture propagation in naturally fractured rock: A three-dimensional approach. Society of Petroleum Engineers. SPE 127888. Meng, H. Z., and Brown, K. E., 1987. Coupling of production forecasting, fracture geometry requirements and treatment scheduling in the optimum hydraulic fracture design. Paper presented at Proceedings 1987 SPE/DOE Joint Symposium on Low Permeability Reservoirs. Meyer, B. R., 1989. Three-dimensional hydraulic fracturing simulation on personal computers: Theory and comparison studies. Society of Petroleum Engineers. SPE 19329. Meyer, B. R., 2011. Meyer Fracturing Simulators, User?s Guide, Ninth Edition. Mohaghegh, S., Balanb, B., Platon, V., and Ameri, S., 1999. Hydraulic fracture design and optimization of gas storage wells. Journal of Petroleum Science and Engineering. 23, (3-4): 161-171.   142 Morales, R. H., 1989. Microcomputer analysis of hydraulic fracture behavior with a pseudothree dimensional simulator. Society of Petroleum Engineers Production Engineering, SPEPE, February, 69-74. Morrill, J. C., and Miskimins, J. L., 2012. Optimizing hydraulic fracture spacing in unconventional shales. Society of Petroleum Engineers. SPE 152595. Munjiza, A., 2004. The Combined Finite-Discrete Element Method. Chichester: J. Wiley & Sons. 348 pp.  Nagel, N. and Sanchez-Nagel, M., 2011. Stress shadowing and microseismic events: A numerical evaluation. Society of Petroleum Engineers, SPE 147363. Nelson, R. A. 1985. Geological analysis of naturally fractured reservoirs. Houston, TX, Gulf Publishing, 320 pp. Ohtake, M., 1974. Seismic activity induced by water injection at Matsushiro, Japan. Journal of Physics of the Earth. 22: 163-176. Olson, J. E., Bahorich, B., and Holder, J., 2012. Examining hydraulic fracture-natural fracture interaction in hydrostone block experiments. Society of Petroleum Engineers. SPE 15261. Pakiser, L. C., Eaton, J. P., Healy, J. H., and Releigh, C. B., 1969. Earthquake prediction and control, Science. 166: 1467-1474. Pearson, C., 1981. The Relationship between microseismicty and high pore pressures during hydraulic stimulation experiments in low permeability granitic rocks. Journal of Geophysical Research. 86: 7855-7864.   143 Perkins, T., and Kern, L., 1961. Widths of hydraulic fractures. Journal of Petroleum Technology. 222:937-949. Potyondy, D. O., and Cundall, P.A., 2004. A bounded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences. 41:1329-1364. Rahman, M. M., Rahman, M. K., and Rahman S. S., 2001. An integrated model for multiobjective design optimization of hydraulic fracturing. Journal of Petroleum Science and Engineering. 31(1): 41-62. Rahman, M. M., Rahman, M. K., and Rahman S. S., 2003. Multicriteria hydraulic fracturing optimization for reservoir stimulation. Petroleum Science and Technology. 21, (11-12): 1721-58. Reine, C. A., and Rory, B.D., 2011. Weighing in on the seismic scale: The use of seismic fault measurements for constructing discrete fracture networks in Horn River Basin. Paper presented at CSPG CSEG CWLS Convention. With permission. Reynolds, M., Shaw, J., and Pollock, B., 2004. Hydraulic fracture design optimization for deep foothills gas wells. Journal of Canadian Petroleum Technology. 43, (12): 5-11. Rockfield Software, 2009. ELFEN Version 4.3.3, Swansea, UK. Rogers, S., Elmo, D., Dunphy, R., and Bearingr, D., 2010. Understanding hydraulic fracture geometry and interactions in the Horn River Basin through DFN and numerical modeling. . Canadian Society for Unconventional Gas and Society of Petroleum Engineers. SPE 137488.   144 Settari, A. and Cleary, M. P., 1986. Development and testing of a pseudo-three-dimensional model of hydraulic fracture geometry. SPE, Trans. AIME, (281): 449-466. Shaffer, R. J., Thorpe, R. K., Ingraffea, A. R., and Heuze, F. E., 1984. Numerical and physical studies of fluid driven fracture propagation in jointed rock. SPE Unconventional Gas Recovery Symposium, Pittsburgh, Pennsylvania. Shah, K. R., Carter, B. J. and Ingraffea, A. R., 1997. Hydraulic fracturing simulation in a parallel computing environment. International Journal of Rock Mechanics and Mining Sciences. 34:(3?4), p 282. Simpson, D. W., Leith, W. S., and Scholz, C. H., 1988. Two types of reservoir induced seismicity. Bulletin of the Seismological Society of America. 78: 2025-2040. Stein, S. and Wysession, M., 2002. An Introduction to Seismology, Earthquakes and Earth Structure. Blackwell, Oxford. Talwani, P. and Acree, S., 1985. Pore pressure diffusion and the mechanism of reservoir-induced seismicity,  PAGEOPH. 122: 947-965. Teng, T. L., Real, C. R., and Henyey, T. L., 1973. Microearthquakes and water flooding in Los Angeles, Bulletin of Seismological Society of America. 63: 859-875. Valko, P., and Economides, M. J., 1995. Hydraulic Fracture Mechanics, John Wiley & Sons, Ltd., New York. Vandamme, L., and Curran, J. H., 1989. A three-dimensional hydraulic fracturing simulator. International Journal of Numerical Methods in Engineering. 28, 900-927.   145 Ventura, J. L., 1985. Slator ranch fracture optimization study. Journal of Petroleum Technology.  37, (8): 1251-1262. Vermylen, J. P., and Zoback, M. D., 2011. Hydraulic fracturing, microseismic magnitudes, and stress evolution in the Barnett Shale, Texas, USA. Society of Petroleum Engineers, SPE 140507. Warpinski, N. R., Teufel, L. W., 1987. Influence of geologic discontinuities on hydraulic fracture propagation. Journal of Petroleum Technology. 39(2): 209-220. Warpinski, N. R., and Branagan, P. T., 1989. Altered-stress fracturing. Society of Petroleum Engineers, SPE 17533. Warpinski, N. R., Moschovidis, Z. A., Parker, C. D., and About-Sayed, I. S., 1994. Comparison study of hydraulic fracturing models-Test case: GRI staged field experiment no. 3. SPE Prod. Facil., Feb., 7-16. Warpinski, N. R., Mayerhofer, M. J., Vincent, M. C., Cipolla, C. L., Lolon, E. P., 2008. Stimulating unconventional reservoirs: maximizing network growth while optimizing fracture conductivity. Society of Petroleum Engineers, SPE 114173. Warren, W. E. and Smith, C. W., 1985. In situ stress estimations from hydraulic fracturing and direct observation of crack orientation. Journal of Geophysical Research, 90:6829-6839. Wong, S., Geilikman, M., Xu, G., 2013. Interaction of multiple hydraulic fractures in horizontal wells.   146 Wu, R., Kressee, O., Weng, X., Cohen, C., and Gu, H., 2012. Modeling of interaction of hydraulic fractures in complex fracture. Society of Petroleum Engineers, SPE 152052. Yost, A. B., and Overbey, W. K., 1989. Production and stimulation analysis of multiple hydraulic fracturing of a 2000-ft horizontal well. Society of Petroleum Engineers, SPE 19090. Zoeback, M. D., and Harjes, H. P., 1997. Injection induced earthquakes and the crustal stress at 9 km Depth at the KTB deep drilling site, Germany. Journal of Geophysical Research. 102(B8): 18477-18491. Zoback, M. D., 2007. Reservoir Geomechanics. 1st ed. United Kingdom: Cambridge University Press, Cambridge. Zoback, M. D., Gorelick, S.M., 2012. Earthquake triggering and large-scale geologic Storage of carbon dioxide. Proceeding of the National Academy of Science of the United States of America (PNAS). 109(26): 10164-10168.With permission.         147 Appendix The following presents the UDEC models on interaction of hydraulic fracture with natural fracture for the sensitivity analyses performed to find a bounding relationship as a function of approach angle and stress state shown in Chapter 2, Section 2.5.3.4.                      148  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 23 degrees approach angle under differential stress of 13 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 14 degrees approach angle under differential stress of 10 MPa.   149  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 26 degrees approach angle under differential stress of 12 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 32 degrees approach angle under differential stress of 6 MPa.   150  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 29 degrees approach angle under differential stress of 10 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 40 degrees approach angle under differential stress of 5 MPa.   151  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 24 degrees approach angle under differential stress of 10 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 48 degrees approach angle under differential stress of 5 MPa.    152  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 30 degrees approach angle under differential stress of 9 MPa.  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 60 degrees approach angle under differential stress of 3 MPa.     153  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 40 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 70 degrees approach angle under differential stress of 3 MPa.    154  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 50 degrees approach angle under differential stress of 4 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 80 degrees approach angle under differential stress of 3 MPa.    155  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 57 degrees approach angle under differential stress of 4 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 85 degrees approach angle under differential stress of 3 MPa.    156  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 60 degrees approach angle under differential stress of 4 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 43 degrees approach angle under differential stress of 6 MPa.    157  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 30 degrees approach angle under differential stress of 11 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 34 degrees approach angle under differential stress of 10 MPa.    158  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 50 degrees approach angle under differential stress of 8 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 33 degrees approach angle under differential stress of 13 MPa.    159  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 38 degrees approach angle under differential stress of 12 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 27 degrees approach angle under differential stress of 13 MPa.    160  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 40 degrees approach angle under differential stress of 11 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 50 degrees approach angle under differential stress of 8 MPa.    161  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 30 degrees approach angle under differential stress of 11 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 46 degrees approach angle under differential stress of 6 MPa.   162  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 50 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 31 degrees approach angle under differential stress of 10 MPa.    163  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 58 degrees approach angle under differential stress of 7 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 52 degrees approach angle under differential stress of 5 MPa.   164  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 58 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 56 degrees approach angle under differential stress of 5 MPa.   165  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 60 degrees approach angle under differential stress of 5 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 85 degrees approach angle under differential stress of 5 MPa.    166   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 65 degrees approach angle under differential stress of 4 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 80 degrees approach angle under differential stress of 4 MPa.   167   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 70 degrees approach angle under differential stress of 5 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 85 degrees approach angle under differential stress of 4 MPa.   168   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 67 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 40 degrees approach angle under differential stress of 7 MPa.   169  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 70 degrees approach angle under differential stress of 4 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 69 degrees approach angle under differential stress of 6 MPa.    170  Modeled response of the induced hydrofrac as it intersects a natural fracture at an 80 degrees approach angle under differential stress of 5 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 40 degrees approach angle under differential stress of 13 MPa.    171  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 50 degrees approach angle under differential stress of 9 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 51 degrees approach angle under differential stress of 12 MPa.    172  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 60 degrees approach angle under differential stress of 8 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 85 degrees approach angle under differential stress of 6 MPa.    173  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 41 degrees approach angle under differential stress of 12 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 85 degrees approach angle under differential stress of 8 MPa.    174  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 75 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at a 52 degrees approach angle under differential stress of 9 MPa.   175  Modeled response of the induced hydrofrac as it intersects a natural fracture at a 71 degrees approach angle under differential stress of 6 MPa.   Modeled response of the induced hydrofrac as it intersects a natural fracture at an 80 degrees approach angle under differential stress of 7 MPa.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0074048/manifest

Comment

Related Items