INVESTIGATION OF INDUCED SEISMICITY MECHANISMS AND MAGNITUDE DISTRIBUTIONS UNDER DIFFERENT STRESS REGIMES, GEOMECHANICAL FACTORS, AND FLUID INJECTION PARAMETERS by Afshin Amini M.Sc., University of Tehran, 2011 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate and Postdoctoral Studies (Geological Engineering) The University of British Columbia (Vancouver) November 2020 © Afshin Amini, 2020 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: Investigation of Induced Seismicity Mechanisms and Magnitude Distributions under Different Stress Regimes, Geomechanical Factors, and Fluid Injection Parameters submitted by Afshin Amini in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Geological Engineering Examining Committee: Dr. Erik Eberhardt - Professor at Department of Earth, Ocean and Atmospheric Sciences – University of British Columbia Supervisor Dr. Marc Bustin - Professor at Department of Earth, Ocean and Atmospheric Sciences– University of British Columbia Supervisory Committee Member Dr. Roger Beckie - Professor at Department of Earth, Ocean and Atmospheric Sciences– University of British Columbia University Examiner Dr. Davide Elmo - Associate professor at Norman B. Keevil Institute of Mining Engineering– University of British Columbia University Examiner Additional Supervisory Committee Members: Dr. Michael Bostock- Professor at Department of Earth, Ocean and Atmospheric Sciences– University of British Columbia Supervisory Committee Member Dr. Lori Kennedy- Associate professor at Department of Earth, Ocean and Atmospheric Sciences – University of British Columbia Sciences Supervisory Committee Member iii Abstract Hydraulic fracturing is the primary means of developing and completing production wells targeting unconventional gas reservoirs. However, these operations have been subject to public concern and regulatory oversight over their rapid growth in use, environmental footprint and potential hazards. This thesis is motivated by one of these potential hazards, induced seismicity, and the need to better understand the level of hazard present with respect to the event magnitudes possible, recognizing that these are not equal across different shale gas plays due to regional differences in the geological conditions or even within the same play due to local differences. The research presented focuses on three objectives: 1) to investigate the effect of different tectonic stress regimes on the magnitude distribution of induced seismicity events; 2) to analyze the influence of fluid injection rate and volume on induced seismicity; and 3) to compare the influence of different geological and operational parameters on induced seismicity in the Montney play in northeastern British Columbia. The methodology integrates statistical and machine learning analysis techniques with advanced 3-D numerical modelling. The empirical analyses are applied to compiled databases of recorded seismicity, hydraulic fracturing well data, and available geology and in situ stress data. Focus is placed on the Montney, but comparisons are also made to other shale gas basins. These analyses are complemented by 3-D numerical modelling used to provide mechanistic understanding to the trends observed in the empirical analyses. The main findings of this thesis are as follows. 1) Thrust faulting stress regimes have lower b-values than strike-slip stress regimes and therefore are more susceptible to larger induced seismicity events. 2) Specific to Montney, injection volume influences susceptibility to induced iv seismicity for wells that target naturally fractured formations, such as Middle Montney formation, whereas injection rate was seen to be an influencing factor for wells targeting the Upper Montney formation. 3) Machine learning provides a valuable means to assess the importance of different geological and operational parameters on induced seismicity for individual shale gas plays, allowing for susceptibility maps and mitigation options to be determined in advance. v Lay Summary Oil and gas industry activities to produce energy can cause earthquakes. In order to protect people and infrastructure threatened by these earthquakes it is necessary to improve our understanding of the parameters that causes these earthquakes to occur and influence how destructive they can get. In order to do so, two novel approaches can be utilized: computer-based simulations and statistical analysis of the available data. Many examples of these studies exist. However, these studies either focus on a single parameter which is highly unlikely or fail to provide a mechanistic understanding for the statistical observations. This thesis, through the development of a new numerical model, as well as statistical analysis of the recorded earthquake and oil and gas activity data, addresses these two knowledge gaps. vi Preface A portion of Section 1.3.1 has been previously published in the Encyclopedia of Engineering Geology: • Eberhardt E., Amini A. (2018). Hydraulic fracturing. In: Bobrowsky P., Marker B. (eds) Encyclopedia of Engineering Geology. Encyclopedia of Earth Sciences Series. Springer, pp. 489-495. https://doi.org/10.1007/978-3-319-12127-7_159-1 A version of Chapter 2 has been published in the Journal of Petroleum Science and Engineering: • Amini, A., Eberhardt, E. (2019). Influence of tectonic stress regime on the magnitude distribution of induced seismicity events related to hydraulic fracturing. Journal of Petroleum Science and Engineering, 182 (106284). https://doi.org/10.1016/j.petrol.2019.106284 I am the first and corresponding author of this research and completed this research under the supervision of Prof. Erik Eberhardt, who provided research guidance. A version of Chapter 3 has been orally presented and published at the 2018 SPE Canada Unconventional Resources Conference: • Amini, A., Eberhardt, E. (2018). Empirical and numerical investigation of the effects of hydraulic fracturing injection rate on the magnitude distribution of induced seismicity events. Society of Petroleum Engineers. doi:10.2118/189820-MS vii Chapter 4 has not yet been published. The work represented in it was completed by myself under the supervision of Prof. Erik Eberhardt. viii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ....................................................................................................................... viii List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................. xiii List of Abbreviations ................................................................................................................. xix Acknowledgements ......................................................................................................................xx Chapter 1: Introduction ................................................................................................................1 1.1 Research Motivation ....................................................................................................... 1 1.2 Thesis Objectives ............................................................................................................ 7 1.2.1 Influence of Tectonic Regime ..................................................................................... 7 1.2.2 Influence of Operational Factors ................................................................................ 8 1.2.3 Susceptibility to Induced Seismicity ........................................................................... 8 1.3 Background Information ................................................................................................. 9 1.3.1 Hydraulic Fracturing ................................................................................................... 9 1.3.1.2 History............................................................................................................... 12 1.3.1.3 Applications ...................................................................................................... 14 1.3.2 Empirical Experiences .............................................................................................. 15 1.3.2.1 Tectonic Stress Regime..................................................................................... 15 1.3.2.2 Operational Parameters ..................................................................................... 18 ix 1.3.2.3 Geological Factors ............................................................................................ 21 1.3.3 Numerical Modelling of Induced Seismicity ............................................................ 24 1.4 Thesis Structure and Outline ......................................................................................... 29 Chapter 2: Influence of Tectonic Stress Regime on the Magnitude Distribution of Induced Seismicity Events Related to Hydraulic Fracturing .................................................................31 2.1 Introduction ................................................................................................................... 31 2.2 Background ................................................................................................................... 33 2.3 Empirical Database Study ............................................................................................. 35 2.3.1 Methodology and Databases Compiled .................................................................... 35 2.3.2 Empirical Analysis and Results ................................................................................ 42 2.3.3 Sensitivity of Results ................................................................................................ 45 2.4 Numerical Modelling .................................................................................................... 47 2.4.1 Modelling Method and Setup.................................................................................... 47 2.4.2 Modelling Results ..................................................................................................... 51 2.5 Discussion ..................................................................................................................... 54 2.6 Conclusions ................................................................................................................... 62 Chapter 3: Empirical and Numerical Investigation on the Influence of Fluid Injection Volume and Rate on Induced Seismicity in the Montney Formation, Northeastern British Columbia .......................................................................................................................................64 3.1 Introduction ................................................................................................................... 64 3.2 Geological Setting of Montney ..................................................................................... 67 3.3 Empirical Study ............................................................................................................ 68 3.3.1 Data Compilation and Preparation ............................................................................ 68 x 3.3.2 Statistical Analysis .................................................................................................... 72 3.4 Numerical Modelling Study .......................................................................................... 80 3.4.1 Model Setup .............................................................................................................. 80 3.4.2 Model Results ........................................................................................................... 85 3.5 Discussion ..................................................................................................................... 90 3.6 Conclusions ................................................................................................................... 95 Chapter 4: Machine learning analysis of factors influencing induced seismicity susceptibility in the Montney Formation of northeastern British Columbia .........................98 4.1 Introduction ................................................................................................................... 98 4.2 Background ................................................................................................................. 100 4.3 Data Compilation and Preparation .............................................................................. 105 4.4 Machine Learning Algorithm Development ............................................................... 109 4.5 Results ......................................................................................................................... 113 4.5.1 Feature Importance ................................................................................................. 113 4.5.2 Model Validation .................................................................................................... 116 4.6 Discussion ................................................................................................................... 120 4.7 Conclusions ................................................................................................................. 123 Chapter 5: Conclusions .............................................................................................................125 5.1 Summary ..................................................................................................................... 125 5.1.1 Influence of Tectonic Stress regime ....................................................................... 126 5.1.2 Influence of Operational Parameters ...................................................................... 128 5.1.3 Susceptibility to Induced Seismicity ....................................................................... 130 5.2 Key Scientific Contributions....................................................................................... 131 xi 5.3 Future Works .............................................................................................................. 133 References ...................................................................................................................................136 Appendix A: Data and Resources .............................................................................................156 Appendix B: Feature Definitions ..............................................................................................157 xii List of Tables Table 2.1: Breakdown of identified induced seismicity events for different major shale gas plays in North America. ......................................................................................................................... 41 Table 2.2: Sensitivity of the spatial and temporal filter used on the b-value of induced events in the Montney and Woodford shale gas plays. ................................................................................ 46 Table 2.3: Discontinuity strength and deformation properties used in the 3DEC numerical simulations. ................................................................................................................................... 49 Table 2.4: Horizontal to vertical in situ stress ratios modelled for each tectonic stress regime. .. 50 Table 2.5: Annual total injection volumes reported for the Woodford and Montney. ................. 58 Table 3.1: Associated wells with Induced Seismicity Separated by Hydraulic Fracturing Target Formation ...................................................................................................................................... 75 Table 3.2: Discontinuity strength and deformation properties used in numerical simulations .... 82 Table 4.1: Input features used for machine learning analysis. .................................................... 108 Table 4.2: Relative ranking of feature importance based on the different machine learning results. ......................................................................................................................................... 116 Table 4.3: Variance of some of the features used in this study .................................................. 122 xiii List of Figures Figure 1.1: Factors that have shaped the world energy market (“World Energy Scenarios,” 2016).......................................................................................................................................................... 2 Figure 1.2: Recoverable unconventional natural gas plays around the globe (WEC, 2013). ......... 3 Figure 1.3: Shale gas play in North America (DOE/EIA, 2016). ................................................... 4 Figure 1.4: Schematic of a hydraulic fracturing operation and formation breakdown. ................ 10 Figure 1.5: Relationship between hydraulic fracture opening and propagation directions relative to the in situ principal stresses, assuming these are aligned with vertical and horizontal (after Eberhardt & Amini, 2018). ........................................................................................................... 11 Figure 1.6: First commercial fracturing treatment on March 17, 1949 (after Smith and Montgomery, 2015). ..................................................................................................................... 13 Figure 1.7: Location and outline of the Western Canadian Sedimentary Basin (WCSB). ........... 16 Figure 1.8: Stress map of North America derived from the World Stress Map (Heidbach et al. 2010). ............................................................................................................................................ 16 Figure 1.9: a) Plot of events in southern California, from the SCSN catalogue. In all frames the green, red and blue lines mark the b values of mainly normal, strike-slip and thrust events, respectively; the grey line marks the average b value and the vertical bars indicate the standard error. b) Plot of events in southern California. The circles at the top of the frame show the range of rake λ used for computing the b values. Inset, circle explaining the rake values and the corresponding colours of the classes of events. c) Frequency–magnitude distributions for pure normal (green) and pure thrust (blue) events of the SCSN and Harvard catalogues. d) Same plot xiv as in a) but using data from the Harvard earthquake catalogue. e) Same plot as in b) but using data from the Harvard earthquake catalogue. From Schorlemmer et al. (2005). .......................... 18 Figure 1.10: Maximum seismic moment and magnitude as functions of total volume of injected fluid from the start of injection until the time of the largest induced earthquake. The equation along the solid line relates the upper bound seismic moment to the product of the modulus of rigidity and the total volume of injected fluid. From McGarr (2014). ......................................... 20 Figure 1.11: Force-displacement and constitutive behaviour of interacting blocks used in the distinct-element method codes UDEC and 3DEC. ....................................................................... 28 Figure 1.12: Fully-coupled hydro mechanical solution procedure used in 3DEC to model potential fracture slip in response to injected fluid pressures. ...................................................... 29 Figure 2.1: Stress map of North America using data from the World Stress Map (Heidbach et al., 2018). ............................................................................................................................................ 36 Figure 2.2: Compilation of reported focal mechanisms for North American earthquakes since 2009 used here to indicate the general tectonic stress regimes for major North American shale basins outlined in green. ............................................................................................................... 38 Figure 2.3: Induced seismicity events and their magnitudes related to hydraulic fracturing and waste water injection, identified by cross-correlating the compiled databases and applying a 10-km spatial filter and 3-month temporal filter relative to well locations and activities, respectively. These are superimposed over a map of major North American shale basins outlined in green. .. 40 Figure 2.4: Magnitude distribution and calculated b-values for induced seismicity events in the Montney versus the Woodford...................................................................................................... 43 Figure 2.5: Histogram of induced events in Montney and Woodford and their distribution curves........................................................................................................................................................ 45 xv Figure 2.6: General modelling procedure used for 3DEC simulations. ........................................ 48 Figure 2.7: 3DEC distinct-element model developed for the numerical simulations, representing a rock mass block with edge lengths of 400 m, cut by a persistent fault (blue plane). The injection point is indicated in yellow and is located 45 meters above and normal to the fault (see left plot). The location of the monitoring points on the fault are shown in the right plot. ........... 51 Figure 2.8: Results of numerical simulations for injection induced fault slip under four different stress regimes. ............................................................................................................................... 52 Figure 2.9: Modelled shear stress changes (i.e., stress drop) for representative monitoring points along the fault (see Figure 2.7) for each stress regime modelled. ................................................ 54 Figure 2.10: Histograms comparing the depth of wells in the Montney and Woodford. ............. 57 Figure 2.11: Average monthly injection volumes per well (solid lines) and overall averages (dashed lines) for the Montney and Woodford basins. ................................................................. 58 Figure 2.12: Sketch of Carboniferous-Permian Dawson Creek Graben Complex (Barclay et al., 1990). ............................................................................................................................................ 60 Figure 2.13: Focal mechanisms of earthquakes in the Montney basin relative to the axis of the Fort St. John Graben (FSJG). ........................................................................................................ 61 Figure 2.14: Location and magnitude of induced seismicity events in the Montney relative to the Fort St. John Graben (FSJG). ........................................................................................................ 62 Figure 3.1: Spatial distribution of earthquakes (red circles) and hydraulic fracturing wells (black circles) included in the respective databases compiled for the Montney play in northeastern British Columbia. .......................................................................................................................... 70 Figure 3.2: Timing and magnitude of induced seismicity events (red circles) superimposed over the daily fluid injection volumes for hydraulic fracturing wells in the Montney (black bars). .... 71 xvi Figure 3.3: Number of induced seismicity events versus cumulative injection volume for hydraulic fracturing wells associated with induced seismicity in the Montney. Compared are the data (red points) and the linear regression fit (blue line). The goodness-of-fit is reported with respect to the R2 value (i.e. coefficient of determination). ........................................................... 71 Figure 3.4: Histograms showing the distribution of injection rates (left), injection volumes (middle) and injection pressure (right) for those wells associated with induced seismicity events (shaded red) compared to the parent data set encompassing all hydraulic fracturing wells in the Montney (shaded grey). ................................................................................................................ 73 Figure 3.5: Histograms showing the distribution of injection rates (left) and injection volumes (middle) and maximum injection pressure (right) in the Montney for those wells associated with induced seismicity events (shaded red) compared to those that are non-seismogenic (shaded grey). The dashed lines indicate the median for each distribution. .............................................. 74 Figure 3.6: Histogram of operational parameters, including distribution of injection volumes and rates, for the Upper and Lower Montney. Plotted are the data for the seismogenic wells (shaded red) and all wells (shaded grey). The P values for comparing these two distributions are inset [KS, Kolmogorov-Smirnov; MW, Mann Whitney]. The dashed lines show the respective average values for the seismogenic wells (red) and for all wells (black). .................................................. 76 Figure 3.7: Box plot of injection volumes and rates for all hydraulic fracturing wells in the Montney as well as the subsets of seismogenic wells in both the Upper and Middle Montney. For each box, the bottom and top edges indicate the 25th and 75th percentiles, respectively, and the central mark indicates the median. The whiskers extend to the most extreme data points not considered outliers. ....................................................................................................................... 78 xvii Figure 3.8: 3-D distinct element model developed for the numerical simulations, representing a rock mass with edge lengths of 1600 m, cut by a persistent thrust fault and a network of orthogonal incipient and natural fractures (lower diagram). The upper diagram shows a cut away highlighting the through going dipping fault (shaded red) and the location of the fluid injection point small (yellow sphere), which is located 90 meters above and normal to the fault surface. 84 Figure 3.9: 3-D distinct-element results showing the simulated hydraulic fracture and invaded zone for the two modeled reservoir geometries: (a) without natural fractures and incipient fractures only (Upper Montney type case), and (b) with both natural and incipient fractures (Middle Montney type case). The fluid injection rate for both models was 6 m3/min, which is approximately the average injection rate used in the Montney for hydraulic fracturing operations. (c-d) Conceptual illustration of hydraulic fracturing, leak-off and secondary fracturing resulting in an invaded zone of elevated pore pressures (shaded grey) for scenarios cases (c) without and (d) with the presence of an interconnected natural fracture network. ........................................... 86 Figure 3.10: (a) Pore pressure distribution on fault plane and (b) corresponding shear slip response to a fluid injection rate of 6 m3/min for 30 minutes for the model scenario without natural fractures (i.e., incipient fractures only). (c) Pore pressure distribution on fault plane and (d) corresponding shear slip response to the same fluid injection rate and duration for the model scenario with natural fractures. ..................................................................................................... 88 Figure 3.11: Magnitude distribution of events calculated from numerical simulations of fault slip for the different fluid injection rates modeled. ............................................................................. 90 Figure 3.12: Injection volume analysis of wells associated with induced seismicity in the Montney for events with M<3, M>3 and considering all events. The dashed line marks the average injection volume of all wells in the database. .................................................................. 94 xviii Figure 3.13: Injection rate analysis of wells associated with induced seismicity in the Montney for events with M<3, M>3 and considering all events. The dashed line marks the average injection volume of all wells in the database. ............................................................................... 95 Figure 4.1: Splitting of the compiled dataset into training and testing sub-sets, with the testing data to be used for validation. ..................................................................................................... 111 Figure 4.2: Results of the 50-fold cross-validations to determine the optimum tree depth for each algorithm, showing from top to bottom: Decision Tree, Random Forest, and Gradient Boost. . 112 Figure 4.3: -Feature importance calculated using three different supervised machine learning algorithms: Decision Tree, Random Forest, and Gradient Boost. .............................................. 115 Figure 4.4: Calculation and use of a confusion matrix to compare test data against correlations derived from machine learning using training data. ................................................................... 117 Figure 4.5: Confusion matrixes comparing predicted versus actual results for the trained decision tree, random forest and gradient boost models against the validation test data. IS: induced seismicity events; No-IS: no induced seismicity events. ............................................................ 118 Figure 4.6: Results for a SHAP feature importance analysis using the Random Forest trained model........................................................................................................................................... 120 xix List of Abbreviations DEM Distinct Element Method EGS Enhanced Geothermal Systems FCN Fully Convolutional Network GWP Gross World Productivity KS Kolmogorov-Smirnov MW Mann-Whitney NGL Natural Gas Liquid NEBC Northeast British Columbia NCSN Northern California Seismic Network SHAP SHapley Additive exPlanations SCSN Southern California Seismic Network WCSB Western Canadian Sedimentary Basin WEC World Energy Council xx Acknowledgements I am very fortunate to work closely with my supervisor and advisor Erik Eberhardt without whom this research project would not have been possible. His guidance helped me to see and not to lose track of the big picture in these years while giving me the opportunity to explore different ideas. His wealth of knowledge and expertise, valuable advices and suggestions helped me bring this study into success. I would like to thank my lovely wife Maral for her constant support and encouragement throughout these years. You are my best friend and unshakable wall which I can lean on. Also, thanks to my parents and brother for their support and friendship. I am very grateful to have the opportunity to collaborate with Steve Rogers, Stuart Venables, Michelle Gaucher. Thank you for your great support and guidance. Thanks to my thesis and examining committee: Marc Bustin, Michael Bostock and Lori Kennedy. Your comments and guidance are greatly appreciated. I would like to thank my office mates: Masoud Rahjoo, Jordan Aaron, Andrew Mitchell, Negar Ghahramani, Ali Mehrabifard, Mostafa Gorjian, Sophia Zubrycky, Sahar Ghahremani, Erika Schmidt, John Whittall, Justin Roy, Giona Preisig, Siobhan Whadcoat, and Christina Brueckman for their great ideas, friendship and support. xxi Dedication To my Mother Your love is always present 1 Chapter 1: Introduction 1.1 Research Motivation According to World Energy Council (WEC) report in 2016 (“World Energy Scenarios,” 2016), the global population is one of the most significant factors that has shaped world energy supply and performance alongside new technologies and productivity, environmental priorities and international governance and geo-political relationships (Figure 1.1). The world’s population has doubled from 3.7 billion people in 1970 to 7.4 billion in 2015 and it is expected to increase by another 25% over the coming decades. At the same time, the period of 1970 to 2015 was one of remarkable world economic growth—Gross World Productivity (GWP) grew 4.4 fold, or 3.3% per annum. High population growth and rapid growth in the labour force was complemented with a high rate of productivity growth, which meant an increase in worldwide demand for energy. The WEC’s 2013 report (WEC, 2013) showed that in 2011, fossil fuel resources met 82% of the world’s energy needs, and although this share was projected to decrease to 76% by 2020, it is clear that it is still the primary contributor to the world’s energy supply. The continued increase in worldwide demand for energy has pressured innovation in the natural gas industry, specifically that related to making low permeability reservoirs (e.g., shales, siltstones, etc.) that were previously seen as being too difficult and costly to produce, feasible; these reservoirs are termed unconventional reservoirs. There are approximately 700 known shale formations worldwide in more than 150 basins (Figure 1.2). At present only a small number of these have been properly assessed regarding their production potential, most of which are in North America 2 (Figure 1.3). The potential volumes of shale gas are enormous and their development has already begun to significantly reshape the natural gas markets worldwide. Figure 1.1: Factors that have shaped the world energy market (“World Energy Scenarios,” 2016). 3 Figure 1.2: Recoverable unconventional natural gas plays around the globe (WEC, 2013). 4 Figure 1.3: Shale gas play in North America (DOE/EIA, 2016). The key technological advances that have made unconventional shale gas reservoirs (also referred to as “tight” gas) accessible and economically feasible include those related to horizontal drilling and hydraulic fracturing. Hydraulic fracturing is the primary means of increasing and maintaining well productivity. The process involves pumping of fluids under high pressure into isolated stages of a horizontal wellbore to first break the reservoir rock and create new fractures, and then by pumping proppant into these and other existing fractures, to prop them open in order to increase the permeability of the rock mass and therefore enhance resource 5 recovery. With the technological advances in multistage hydraulic fracturing and horizontal drilling over the past 10 years, there has been a sharp increase in natural gas production from shale gas plays, which are expected to continue to grow (DOE/EIA, 2016). Although hydraulic fracturing has application in many industries and is used worldwide, it has attracted public concerns over its pace and environmental footprint. These environmental impacts include induced seismicity, potential contamination of groundwater resources, and methane emissions during and after hydraulic fracturing operations (Uddameri et al. 2015; Zoback and Arent 2014). It is noted that induced seismicity is also a key hazard faced by the development of enhanced geothermal systems (EGS), similarly related to fluid injection, and has equally generated public and regulator scrutiny(BC Oil and Gas Commission, 2014). Induced seismicity refers to earthquakes and tremors that are caused by anthropogenic activities, in most cases due to altering the effective stresses acting on a fault, either through a pore pressure increase and/or by reducing the normal stresses. Experiences with induced seismicity in response to fluid injection, and its study, can be traced back to early research in the 1960’s. Healy et al. (1968) presented statistical evidence for a correlation between fluid injection and earthquakes in the Denver region. Since then, numerous studies have studied induced seismicity related to hydraulic fracturing and fluid injection for different purposes: unconventional gas development (Atkinson et al., 2020; Schultz et al., 2020; Davies et al., 2013; Hummel & Shapiro, 2013; Shapiro & Dinske, 2009; Vermylen & Zoback, 2011; Warpinski, Du & Zimmer, 2012; Zang et al., 2013), waste water injection (W. L. Ellsworth, 2013; Healy et al., 1968; Horton, 2012b), and enhanced geothermal development (Preisig et al., 2015; Weidler et al., 2002). Many 6 researchers have also studied different aspects of hydraulic fracturing in the context of shale gas production separate from the question of induced seismicity, such as: effects of in-situ stresses on hydraulic fracture propagation (Vermylen & Zoback, 2011), influence of stress shadows from adjacent hydraulic fractures (Zangeneh et al., 2015), monitoring the extent of hydraulic fracture propagation (Jeffrey et al., 2009), the effects of natural fractures on hydraulic fractures propagation and their interaction (Jeffrey et al. 2009; Ouenes et al. 2013; Zangeneh et al. 2015), microseismic monitoring (Van Der Baan et al. 2013; Maxwell et al. 2002), magnitude distribution of microseismic events (Warpinski et al., 2012), how to distinguish between natural earthquakes and induced seismicity (Cesca et al. 2012) and earthquake physics (Galis et al., 2017). Ultimately, the local geology and differences in operational parameters for a given hydraulic fracturing treatment results in considerable variability in induced seismicity experiences and susceptibility. Key controlling factors such as the presence of a critically stressed fault, or network of faults, and their geometric and strength characteristics relative to the hydraulic fracturing fluid injection volume, rate and pressures, represent significant uncertainty that limit deterministic evaluations of fault slip and maximum seismic magnitudes. With this in mind, assessing magnitude distributions is suggested here in this thesis as being more appropriate compared to considering only maximum magnitude, allowing analyses to be constrained by reasonable lower and upper bounds of controlling parameters based on site specific data. From a statistical perspective, this would include the maximum potential event. Thus, key research questions can be posed regarding induced seismicity susceptibility and sensitivity, including how do different i) tectonic stress regimes (e.g., thrust vs strike-slip), ii) geological factors (e.g., 7 presence of natural fractures, distance to basement, etc.), and iii) completion parameters (e.g., injection volumes and rates) affect the magnitude distribution of induced seismicity responses. Improved understanding of these sensitivities relative to local site conditions would allow designers and operators to better assess the likelihood and therefore risk of large induced seismicity events, and to mitigate this through adjusting the injection volume and/or rate during a hydraulic fracturing treatment to minimize fluid pressure diffusion and subsequent potential of slip along near-field faults. 1.2 Thesis Objectives The overarching objective of the research presented in this thesis is to improve our understanding of the geological and operational parameters affecting the magnitude and magnitude distribution of induced seismic events. To do so, three research objectives are defined as described in the following sub-sections. 1.2.1 Influence of Tectonic Regime The first objective is to enhance our understanding of fault slip and induced seismicity under different tectonic stress regimes. Three in situ stress regimes will be investigated: normal, strike-slip and thrust. This objective will involve two parallel studies. The first will focus on empirical data pertaining to recorded seismicity related to hydraulic fracturing operations for which source mechanisms are reported and/or can be derived. To complement the empirical analysis, a second study focusing on numerical modelling using the 3-D distinct element code 3DEC is performed to simulate fault slip under different stress regimes. Techniques will be developed to model induced seismicity mechanisms and sensitivity, including means to calculate from the model 8 event magnitudes. A comparative analysis will be carried out simulating induced seismicity under different stress field boundary conditions applied relative to the geometrical and operational input parameters assumed. 1.2.2 Influence of Operational Factors Different operators conducting fluid injection, use different operational parameters (e.g., fluid injection volumes and rates) tailored for the local geological and tectonic settings. As such, they are of special interest as they offer a means to mitigate induced seismicity hazards through their control or manipulation. The second objective of this thesis will investigate the sensitivity of fault slip and induced seismicity to operational factors. This will be carried out by integrating empirical data related to hydraulic fracturing activity in northeastern British Columbia (i.e., distribution of typical injection volumes and rates), with numerical modelling using the 3-D distinct-element modelling techniques developed for the previous objective. 1.2.3 Susceptibility to Induced Seismicity A key knowledge gap in understanding induced seismicity and assessing hazard potential is determining the susceptibility of wells to these events. The third and final objective of this thesis will investigate the influence of different geological and operational factors on the spatial and temporal distribution of induced seismicity events. This is a complex task as there are many parameters that can potentially affect induced seismicity, and these can vary for different shale gas plays. Besides the complexity of the task, massive amounts of geological, operational, and seismic data are being collected during hydraulic fracturing activities thus requiring robust analysis methods capable of handling large data sets. The rapid development of data science, and 9 specifically machine learning techniques, offers a potentially useful tool for this purpose. Thus, this objective will investigate the application of different machine learning algorithms in determining the relative importance of geological and operational parameters on triggering induced seismicity, with focus placed on the Montney region in northeastern British Columbia. 1.3 Background Information Despite the considerable attention hydraulic fracturing and induced seismicity have received in the research literature, questions specifically related to magnitude distributions of induced seismicity events have not yet been addressed completely and represent a knowledge gap. For the research presented in this thesis to address this (and build on previous studies), several foundational topics are relevant to be built on as presented in the following sub-sections. 1.3.1 Hydraulic Fracturing 1.3.1.1 Definition Hydraulic fracturing is a well stimulation technique used to generate and open fractures in the rock to increase its permeability. The process involves injecting 'fracking’ fluids (containing primarily water, proppants and other additives) under high pressures into a wellbore to initiate and propagate a fracture or system of fractures (King, 2012). During injection, if fluid is pumped into a well faster than the fluid can escape into the formation, inevitably pressure rises, and at some point the rock strength (plus minimum horizontal stress) will be exceeded and the formation will break (Figure 1.4). The hydraulic fracture generated will open against the minimum principal stress, propagating in the direction of the maximum principal stress. This will result in the generation of vertical fractures where the minimum principal stress aligns with one 10 of the horizontal stresses (Shmin). Where the minimum principal stress aligns with the vertical stress, the hydraulic fracture will propagate in the horizontal plane. This is demonstrated in Figure 1.5. Figure 1.4: Schematic of a hydraulic fracturing operation and formation breakdown (Eberhardt & Amini, 2018). 11 Figure 1.5: Relationship between hydraulic fracture opening and propagation directions relative to the in situ principal stresses, assuming these are aligned with vertical and horizontal (after Eberhardt & Amini, 2018). Today, most hydraulic fracturing treatments involve multi-million dollar operations undertaken by service companies. All of the materials to be used for the hydraulic fracturing treatment (water, proppant, and chemicals) are fed into a blender that mixes the materials into a slurry that is distributed by a manifold to the high-pressure pumps. These high-pressure pumps intensify the fluid pressure to that required to frac the well, where they are then forced down the well through the wellhead and well “frac” string to the perforations opening into an isolated interval of reservoir rock. Key operational factors include the use of proppants, fluid viscosity, pump rates 12 and injection volumes (Eberhardt & Amini, 2018). During pumping, the hydraulic fracture is held open by the fluid pressure. However, once pumping stops and the injection pressure dissipates, the minimum principal stress will act to close the hydraulic fracture created. To prevent this, a propping agent, typically sand, is then pumped with the hydraulic fracturing fluid to maintain an open, conductive fracture. Montgomery and Smith (2015) note that today oil and gas reservoir treatments average approximately 45 metric tonnes of propping agent and 200 m3 of fluid, with large treatments exceeding 2000 metric tonnes of propping agent and 4000 m3 of fluid. Fluid viscosity and pump rate work in unison to control the net pressure to attain the desired hydraulic fracture height and length, as well as to ensure sufficient opening to allow proppant to enter the fracture and carrying velocity to transport the proppant deep into the hydraulic fracture. 1.3.1.2 History Hydraulic fracturing was first introduced in the United States in the late 1940s by Stanolind Oil and later commercialized in the early 1950s by Halliburton to increase production from oil and gas wells (Figure 1.6). The early technique involved injecting a blend of crude oil and gasoline at treatment pressures that would fracture the reservoir rock, together with sand to prop the fractures open ((Smith & Montgomety, 2015)). The first wells treated saw an average increase in production of 75%. This led to a rapid growth in the application of hydraulic fracturing, and by the mid-1950s, more than 100,000 individual treatments had been performed (Hubbert & Willis, 1957). 13 Figure 1.6: First commercial fracturing treatment on March 17, 1949 (after Smith and Montgomery, 2015). For much of its history, hydraulic fracturing treatments were applied in vertical wellbores. However, in the 1990s, hydraulic fracturing saw a step change in its use when it was combined with advancements in horizontal drilling. This greatly expanded the viability of low permeability geological plays, particularly shale gas (Eberhardt & Amini, 2018). As the industry has become more proficient at drilling longer horizontal wells, the lateral length increased steadily, with laterals reaching out several kilometers and numerous hydraulic fracturing stages being placed along the length of the well in contact with the targeted formation. This has resulted in hydraulic fracturing treatments evolving into substantial multistage fracturing programs. In northeastern British Columbia, the average injection volume reported by Johnson and Jonson (2012) varies from 250 m3 per fracture stage and 1,900 m3 per well for the Montney Trend where energized treatments are favoured, to 1,000 m3 per fracture stage and 7,800 m3 per well in the Montney 14 North Trend where hybrid-energized slickwater treatments are favoured. For the Horn River region, where slickwater treatments are favoured, this increases further to 2,000 m3 per fracture stage and 34,900 m3 per well (Allen et al., 2019). 1.3.1.3 Applications Although hydraulic fracturing is commonly associated with the development of unconventional shale gas reservoirs, it is also used in other industries too where there is a need to increase the permeability of the rock, such as water well production enhancement (Williamson & Woolley, 1980), conventional oil and gas production, carbon sequestration, coalbed methane development, coal mine methane drainage, and more recently, enhanced geothermal (hot dry rock). In fact, induced seismicity events of up to magnitude 4.6 have been recorded at The Geysers enhanced geothermal project in Northern California, and induced seismicity has contributed to the cancellation of the Basel Deep Heat Mining enhanced geothermal project in Switzerland. Separate from increasing rock permeability, the fractures generated during hydraulic fracturing also weakens and breaks up stronger, more massive rock. Accordingly, hydraulic fracturing is also used in the mining industry to precondition an orebody to ensure suitable fragmentation for caving and mining (van As & Jeffrey, 2000), and to soften the rock mass to reduce its ability to store strain energy for rock burst mitigation (Feng & Shi, 2012). Lastly, because the fluid injection pressures must overcome the near- and far-field stresses present in the rock to generate and propagate a hydraulic fracture, hydraulic fracturing is commonly used as an in situ stress measurement technique for geotechnical rock engineering design (e.g.; mines, tunnels, dams, foundations, etc.) (ASTM International, 2004). 15 1.3.2 Empirical Experiences 1.3.2.1 Tectonic Stress Regime The main shale basin in Canada is the Western Canadian Sedimentary Basin (WCSB), which extends into northeastern British Columbia (Figure 1.7). The WCSB is a large and tectonically diverse region and encompasses some of the most, and least, seismically active areas in Canada. The basin is bounded by the foreland fold and thrust belt to its west and extends from the Northwest Territories in the north down into the U.S. in the south. The main influences on the tectonics of the basin are the relative motion of the Pacific and North American plates where they meet along the west coast of North America, and the subduction of the Juan de Fuca plate beneath the North American plate (Ristau et al. 2007), also known as the Cascadian subduction zone. The Cascadia subduction zone is a high seismicity region of active plate convergence. The compressive stress regime caused by this subduction zone is considered to be the dominant stress regime in northeastern British Columbia (Figure 1.8) (Heidbach et al., 2010). 16 Figure 1.7: Location and outline of the Western Canadian Sedimentary Basin (WCSB). Figure 1.8: Stress map of North America derived from the World Stress Map (Heidbach et al. 2010). 17 Figure 1.8 also demonstrates that the overall stress regime varies across North America. For example, on the west coast of the U.S., south of the Cascadian subduction zone, a strike-slip regime is dominant (e.g., the San Andreas fault in California). This transitions to normal fault and strike-slip regimes in the midcontinent (e.g., Basin and Range), and finally evolves to a thrust fault regime in the east. This would suggest that shale gas plays in different basins in North America are subject to different tectonic stress regimes. It is also known that natural earthquakes obey a power law size distribution, the slope of which is known as the ‘b value’ (Gutenberg & Richter, 1945). This is commonly used to describe the relative occurrence of small and large events; a higher b value indicates a larger proportion of small earthquakes, and vice versa, a smaller b value indicates a larger proportion of large earthquakes. It has also been shown that b-value depends systematically on earthquake focal mechanism, with smaller b-values being associated with thrust fault earthquakes, intermediate values with strike-slip earthquakes, and larger values with normal fault earthquakes (Schorlemmer et al., 2005) (Figure 1.9). Considering that different shale basins are located in different stress regimes and earthquake studies show that stress regimes (represented by different faulting styles) correspond to different b-values, it is hypothesized that the tectonic stress regime will influence the likelihood and magnitude distribution of induced seismicity. This hypothesis is investigated in Chapter 2. 18 Figure 1.9: a) Plot of events in southern California, from the SCSN catalogue. In all frames the green, red and blue lines mark the b values of mainly normal, strike-slip and thrust events, respectively; the grey line marks the average b value and the vertical bars indicate the standard error. b) Plot of events in southern California. The circles at the top of the frame show the range of rake λ used for computing the b values. Inset, circle explaining the rake values and the corresponding colours of the classes of events. c) Frequency–magnitude distributions for pure normal (green) and pure thrust (blue) events of the SCSN and Harvard catalogues. d) Same plot as in a) but using data from the Harvard earthquake catalogue. e) Same plot as in b) but using data from the Harvard earthquake catalogue. From Schorlemmer et al. (2005). 1.3.2.2 Operational Parameters There are numerous operators conducting hydraulic fracturing and wastewater injection activities in North America, each using different operational parameters (e.g., fluid injection volumes and rates). These are tailored for the local geological and tectonic settings. This raises the question of how do different operational parameters influence the magnitude distribution of induced seismicity events? Typical injection operations for northeastern BC involve average injection 19 rates of 8 m3/min and injection volumes of 800 m3 (BC Oil & Gas Commission, 2018). Operational factors are not limited to injection volume and rates but also involve other parameters such as injection pressure and fluid viscosity. It has been argued that the moment release attributable to induced seismicity is related to the net volume of the injected fluid, with empirical trends established that link an upper limit for the moment magnitude to injection volume (Hallo et al., 2014; McGarr, 2014) (Figure 1.10). Weingarten et al. (2015) investigated this combining information from public databases on injection wells with available earthquake catalogues, and concluded that injection rate is the most important operational parameter affecting induced seismicity. It has been further suggested that stress transfer to nearby faults following an induced seismicity event may lead to subsequent events and thus the volume of injected fluid may not serve to limit event magnitude (Sumy et al., 2014). 20 Figure 1.10: Maximum seismic moment and magnitude as functions of total volume of injected fluid from the start of injection until the time of the largest induced earthquake. The equation along the solid line relates the upper bound seismic moment to the product of the modulus of rigidity and the total volume of injected fluid. From McGarr (2014). When considering induced seismicity in the WCSB, different studies have come to the conclusion that injection volume is a key factor in triggering induced seismicity. Schultz et al. (2018) investigated the relationship between injection parameters and induced seismicity in the Duvernay shale play in Alberta and concluded that events are associated with completions that used larger injection volumes and that seismic productivity scales linearly with injection volume. Their analysis also showed that injection pressure and rate have an insignificant association with seismic response and volume and geological factors account for the majority of the variability in the induced seismicity rate in this region. Other studies such as Kao et al. (2018) showed that 21 injection in areas with moderate tectonic strain may temporarily increase the local seismic hazard but tectonic strain rate is an important control on injection induced earthquakes. It was also shown that the maximum magnitude appears to increase with the total injected volume, but large injection volumes do not always result in large magnitude events and the inference of large injected volume is a necessary, but not a sufficient/sole condition for the occurrence of large magnitude induced seismicity events. Although the above studies draw correlations between induced seismicity and operational factors such as injection volume and rate, it is observed that there are different opinions as to which operational parameter is more influential between different shale gas plays. Also, knowledge gaps and uncertainties remain regarding their influence on triggering induced seismicity and their magnitude distributions. Yet, understanding the influence of these parameters is of primary importance as they offer a means to possibly mitigate induced seismicity hazards by controlling or adjusting the injection pressure, rate and/or volume. The influence of operational parameters on induced seismicity likelihood and magnitude distribution specific to the Montney region of northeastern British Columbia is investigated in Chapter 3. 1.3.2.3 Geological Factors In the early years of studying induced seismicity, many empirical and numerical studies focused on the influence of operational factors such as injection volume and rate. This focus was largely due to a lack of geological data being available for these systems. However, in recent years as more data has become available, subsequent studies have shown the importance of geology and geological factors when investigating induced seismicity. Göbel (2015) compared different fluid 22 injection operations in California and Oklahoma and examined their temporal and spatial seismicity variations. His results suggest that operational parameters for fluid injection are likely of secondary importance and that the primary controls on induced seismicity are the site-specific geology and geological setting. Van der Baan and Calixto (2017) compared current and historic seismicity rates in six U.S. states and three Canadian provinces to past and present hydrocarbon production. Their study showed that in addition to injection volumes, local and regional-scale geology and tectonics influenced increased induced seismicity hazard susceptibility. Other detailed studies include Shah and Keller (2017), who combined geophysical and drill hole data to map subsurface geologic features in the crystalline basement in Oklahoma and showed that most induced seismicity events are located where the crystalline basement is likely composed of fractured intrusive or metamorphic rock. In contrast, areas with extrusive rocks or thick (>4 km) sedimentary cover exhibit little seismicity, even in high injection rate areas, similar to deep sedimentary basins in Michigan and western North Dakota. They concluded that the differences in seismicity may be due to variations in permeability structure: within intrusive rocks, fluids become more narrowly focused in fractures and faults, causing an increase in local pore fluid pressure, whereas more distributed pore space in sedimentary and extrusive rocks may relax fluid injection pressures. Hincks et al. (2018) developed an advanced Bayesian Network to model joint conditional dependencies between spatial, operational, and seismicity parameters in Oklahoma. They found that injection depth relative to the crystalline basement most strongly correlated with seismic moment release and that the combined effects of depth and volume are critical, as injection rates become more influential near the basement interface. This study showed the importance of injection depth when studying induced seismicity. 23 Skoumal et al. (2015) investigated induced seismicity due to hydraulic fracturing operations in Ohio and also showed that based on comparisons between relocated hypocenters, well paths and estimated basement depths, events most likely occurred near the basement contact. Similarly, Currie et al. (2018) investigated the relationships between disposal well locations and operational histories, spatio-temporal patterns of seismicity, and mapped subsurface structures in Washington County, Ohio. Their results indicate that seismicity occurred along faults below the injection interval in the Precambrian crystalline basement. From seismic reflection lines they showed the presence of basement hosted fault systems that cut through the injection intervals targeted by an injection well, and in doing so, provided permeability pathways for fluid pressure increases that initiated slip. In western Canada, Schultz et al. (2016) studied induced seismicity in Alberta and found that hypocenters of induced seismicity clusters coincide with the margins of Devonian carbonate reefs and interpreted this spatial correspondence as the result of geographically biased activation potential, possibly as a consequence of reef nucleation preference to paleobathymetric highs associated with Precambrian basement tectonics. Their work provided evidence that in some areas Paleozoic and Precambrian strata are likely in hydraulic communication with each other, thus pointing to the important role of regional- and local-scale geological factors influencing induced seismicity. Eaton and Schultz (2018) also point to natural/geological processes, including transformation of organic material (kerogen) into hydrocarbon and cracking to produce gas, that can cause fluid overpressures. They present two examples from the WCSB where induced seismicity caused by hydraulic fracturing was strongly clustered within areas characterized by high pore-pressure gradients. This suggests that high in situ overpressure of 24 shale formations is also a controlling factor influencing induced seismicity. Eyre et al. (2019) proposed a new model for fault activation in shales which takes the clay content of the rock hosting the fault into account. In this model, unstable regions of a fault progressively load stable regions by means of aseismic slip making them more susceptible to induced seismicity. Lastly, Pawley et al. (2018) used machine learning model to evaluate the likelihood of a well generating an induced seismicity event based on inputs related to the underlying geology, tectonics, geomechanical properties, and hydrological proxies. Their results suggested that proximity to basement, in situ stress, proximity to fossil reef margins, lithium concentration, and rate of natural seismicity are among the strongest model predictors. Although this study did not provide mechanistic explanations for the parameters determined by the machine learning analysis as being the most influential, it demonstrated the power of machine learning to analyze large complex data sets combining geological, operational and induced seismicity data sets. These techniques are investigated in more detail in Chapter 4 using data for the Montney region of northeastern British Columbia. 1.3.3 Numerical Modelling of Induced Seismicity Data analysis techniques are capable of determining correlations hidden in large, complex data sets. However, they are limited in providing mechanistic insights that explain these correlations. Numerical models on the other hand do provide mechanistic insights by applying physics-based mathematical relationships and solving to obtain the modelled behavior of a complex system in response to a simulated perturbation. The application of numerical modelling to the analysis of rock engineering problems represents a challenging task due to the heterogeneity present in the 25 rock mass, and therefore complexity in the rock mass behaviour. In this context, the available numerical approaches are typically classified as either continuum- or discontinuum-based methods. Continuum-based methods, such as the finite difference and finite element methods, are more commonly used due to their need for less user input regarding the geology (i.e. any discontinuities present are treated implicitly rather than explicitly) and thus faster computing times. An example is that reported by Rutqvist et al. (2013) who used the finite-difference code TOUGH-FLAC (FLAC coupled with the TOUGH2 multiple flow and heat transport simulator) to simulate induced seismicity. They studied the effects of fault permeability on the magnitude of seismic events and showed that for a critically stressed fault which is also permeable, the rupture zone can be more extensive than an impermeable case. Also, their results indicated that the peak and residual coefficients of friction are important parameters that control the potential seismic magnitude. In Rutqvist et al. (2015), they extended their work into 3-D using the same continuum codes. Their simulations showed moment magnitudes ranging from -2.5 to 0.5 with the exception of a M2.3 for a very brittle fault. Their investigation also showed that microseismic event magnitudes increased with simulated depth. Continuum-based methods derive their name from their underlying treatment of the problem domain (and geology) as an interconnected continuum. Thus, the simulation of a fault, which forms a discontinuity, either needs to be treated implicitly or as a discrete interface with separate constitutive relationships that model the interactions that occur along the interface. The latter lends itself to discontinuum-based methods, such as the distinct-element method, which treats the 26 presence of geological discontinuities (e.g., faults, bedding planes, joints, etc.) as a network of interfaces that can slip, open or close in response to the stresses and pore pressures acting across the problem domain. Due to the nature of induced seismicity as deriving from slip along critically stressed discontinuities, discontinuum codes (e.g., UDEC/3DEC, Irazu, PFC, XSite) possess specific advantages over continuum codes (Lisjak & Grasselli, 2014). Amongst the discontinuum methods, the distinct-element code UDEC (Itasca Consulting Group Inc., 1999) has been recently investigated to simulate induced seismicity associated with hydraulic fracturing. Zangeneh et al. (2013) used moment magnitude and seismic moment relationships to calculate the magnitude of induced seismicity from shear displacement on a critically stressed fault. It was shown that even after shut-in, the subsequent diffusion of the pore pressure front can eventually reach and reduce the effective stresses acting on a nearby fault, triggering slip. In their work, they were able to reproduce a M2.45 induced seismicity event. Pirayehgar and Dusseault (2015) extended this work using the same approach to show that the magnitude of induced seismicity events can change with the rock mass fabric pattern and in situ stress ratio. A limitation of these earlier numerical studies by Zangeneh et al., (2013) and Pirayehgar and Dusseault (2015), is that they were restricted to 2-D analyses and 2-D simplifications of the effects of geology relative to the in-situ stress state. For the purpose of this research, numerical analyses are performed in Chapters 2 and 3 using the 3-D distinct-element code 3DECTM (Itasca Consulting Group Inc., 1999). The discontinuum-based formulation of 3DEC uses a Lagrangian numerical technique in which the problem domain is divided through by discontinuities of 27 variable orientation, spacing and persistence. The resulting body of deformable blocks and their contacts may then undergo large deformations in shear and opening in response to a change to the modelled stresses and/or joint water pressures. The blocks, although deformable, remain intact. A detailed discussion of the Distinct Element Method (DEM) is presented by Cundall and Hart (1992). In modelling the behavior of jointed rock masses, 3DEC calculates both the yielding of the rock blocks in shear and tension, together with slip along the network of discontinuities bounding the blocks (Figure 1.11). Progressive failure associated with crack propagation can be modelled by the breaking of contacts pre-defined along the fracture interfaces. To incorporate degrees of freedom in how a fracture can propagate in UDEC/3DEC, outside of the main discontinuity network, voronoi tessellation in 2-D and tetrahedral blocks in 3-D can be used to define smaller potential fracture pathways between the network of natural discontinuities and these smaller segments can be assigned intact rock properties susceptible to fracturing if the induced stresses or fluid injection pressures are high enough (e.g., Zangeneh et al., 2015b). 28 Figure 1.11: Force-displacement and constitutive behaviour of interacting blocks used in the distinct-element method codes UDEC and 3DEC. Key for simulating hydraulic fracturing, 3DEC has the capability to model fluid flow through the fracture network (note that the blocks in this assemblage are treated as being impermeable). A fully coupled hydro-mechanical solution procedure is applied where fracture conductivity is dependent on the mechanical deformation of the fracture apertures via a cubic law relationship, and conversely, the mechanical deformation of the apertures is dependent on the joint water pressures and dilation arising in the event of shear slip or joint opening (Figure 1.12). Thus, the influence of fluid injection on aperture changes and effective stresses can be modelled together with the simulation of hydraulic fracturing creating new fracture permeability and connectivity. Effective stresses accounting for changes in pore pressures are calculated relative to fracture shear strengths to model fracture slip based on Coulomb slip friction laws. 29 Figure 1.12: Fully-coupled hydro mechanical solution procedure used in 3DEC to model potential fracture slip in response to injected fluid pressures. 1.4 Thesis Structure and Outline This dissertation is prepared as a manuscript-based thesis and consists of an introduction chapter, three research chapters, and a conclusion chapter. Chapter 1 introduces the motivation of the research, explains the thesis objectives, and reports the initial literature review carried out. Additional literature review specific to each of the three research chapters are included in each of these chapters. The three research chapters that follow correspond to the three research objectives stated in Section 2 of Chapter 1. 30 Chapter 2 investigates the effect of different tectonic stress regimes on the magnitude and magnitude distribution of induced seismicity in response to hydraulic fracturing operations. Three databases are compiled and analyzed, with focus on North American data, combining earthquake catalogues, hydraulic fracturing well data, and in situ stress data based on dominant focal mechanisms. In parallel, a series of 3-D distinct element numerical simulations are performed to provide mechanistic insights to the interpretation of the results of the empirical analyses. Chapter 3 examines the effect of hydraulic fracturing operational parameters on triggering induced seismicity, with focus placed on data for the Montney play in northeastern British Columbia. An empirical analysis is performed to analyze the influence of injection pressure, rate and volume with induced seismicity, including differences between the upper and lower parts of the Montney formation. Again, the interpretations of the empirical analyses are supported by 3-D distinct-element numerical simulations that model fault slip under different injection rates and volumes. Chapter 4 presents the results of a machine learning analysis conducted for the classification problem of wells either associated or not associated with seismicity. Data is compiled integrating recorded seismicity in the Montney region and different geological parameters and hydraulic fracturing operational parameters. Three different machine learning algorithms are applied to determine and rank the importance of the different geological and operational parameters on triggering induced seismicity in the Montney region. 31 Chapter 2: Influence of Tectonic Stress Regime on the Magnitude Distribution of Induced Seismicity Events Related to Hydraulic Fracturing 2.1 Introduction Hydraulic fracturing is a key technology that is used to increase production from unconventional oil and gas reservoirs. When performing hydraulic fracturing, high pressure fluid is pumped into different sections of a wellbore to first break the reservoir rock and create new fractures, and then to prop these and other existing fractures open in order to increase the permeability of the rock mass and therefore enhance resource recovery. With the technological advances in multistage hydraulic fracturing and horizontal drilling over the past 10 years, there has been a sharp increase in natural gas production from shale gas plays, which are expected to continue to grow (DOE/EIA, 2016). Hydraulic fracturing, however, is not without its impacts on the surrounding natural environment. Among the public and regulator concerns are the hazards and risks related to induced seismicity, either generated directly in response to the hydraulic fracturing treatment or after, during the deep well disposal of wastewater (i.e., flowback) produced from the treatment. Induced seismicity however, is not exclusive to shale gas production; other activities such as salt water disposal (Ellsworth, 2013; Horton, 2012; Keranen et al., 2014), enhanced geothermal systems (Kraft & Deichmann, 2014; Majer et al., 2007; Zang et al., 2014), mining (Hasegawa et al., 1989; J. Holmgren, 2015; Sainoki & Mitri, 2014), and CO2 sequestration (Akber Hassan & Jiang, 2012; Cappa & Rutqvist, 2011; Mazzoldi et al., 2012) are also affected by induced 32 seismicity. In fact, the impoundment of dam reservoirs following their construction has resulted in some of the largest induced seismicity events reported (Guglielmi et al., 2015; Gupta & Chadha, 1995). These each contribute to the experience base and study of induced seismicity. One of the topics related to induced-seismicity hazard assessments is the maximum event magnitude possible. Induced seismicity events related to shale gas and hydraulic fracturing generally have low magnitudes (< 3). However, there have been events reported that reach magnitudes and intensities that can be felt on surface and have even caused damage to buildings (Tagliabue, 2013). In places such as British Columbia, hydraulic fracturing operations are restricted to remote or low population areas and therefore the risk involves civil infrastructure or the infrastructure of the companies conducting the hydraulic fracturing. Therefore, it is important to study the parameters that influence the magnitude of these events and their distribution, and the maxima possible. There are several uncertainties and knowledge gaps regarding our understanding of induced seismicity and its magnitude distribution. Of specific relevance to the investigation presented in this chapter is the question of mechanistic uncertainty – Do we properly account for in situ stress regime in assessing fault slip susceptibility and response? It is postulated that a critical fault in a thrust fault stress regime should respond differently than that in a strike-slip or normal fault stress regime. In this work, I explore and investigate this question by first cross-correlating data in an existing database of stress observations (e.g., borehole breakouts, focal mechanisms) with empirical data 33 I compiled pertaining to recorded seismicity related to hydraulic fracturing operations. This is done relative to the major shale gas plays in North America for which data is publically available. Next, this empirical analysis is supplemented with numerical modelling using the commercial 3-D distinct element code 3DEC (Itasca Consulting Group Inc., 1999) to simulate fault slip under different tectonic stress regimes to further investigate the mechanistic causes for differences observed in the empirical data. 2.2 Background Prior experiences with earthquake magnitudes suggest they obey a power law size distribution, which can be expressed on a frequency-magnitude plot as: 𝐿𝑜𝑔 𝑁 = 𝑎 − 𝑏𝑀 (2.1) where N is the number of earthquakes with magnitudes equal or greater than M, and 𝒂 and 𝒃 are constants related to the fitted line (Gutenberg & Richter, 1945). The parameter 𝒂 describes the total number of earthquakes (y-axis intercept) and 𝒃, known as the ‘b value’, is commonly used to describe the relative frequency of large and small events via the slope of the line. A high b-value indicates a larger proportion of smaller magnitude events, and conversely, a smaller b-value indicates a larger number of higher magnitude events. Studies of b-value have proven to be an effective means of characterizing earthquake magnitude likelihood across multiple scales, from laboratory (Amitrano, 2003; Scholz, 1968) to tectonic regime (Gulia & Wiemer, 2010; Westerhaus et al., 2002) studies. For example, acoustic emissions from microfracturing during uniaxial and triaxial compression rock testing 34 experiments were found to obey the same form of size distribution, and further, that b values decrease linearly with increasing differential stress (σ1-σ3) (Scholz, 1968). However, studies at the tectonic regional scale also show that the total stress magnitude is important, with the b value for earthquakes in the continental crust decreasing at an approximately linear rate with depth (Mori & Abercrombie, 1997; Spada et al., 2013), even though the differential stress with depth decreases. Higher stresses mean more elastic strain energy stored. Van der Elst ( 2016) suggested that nucleation of induced earthquakes are controlled by injection parameters whereas the magnitude of these events are controlled by the tectonics of the region. It has also been shown that b-values for natural earthquakes depend on the earthquake focal mechanism (El-Isa & Eaton, 2014; Gulia & Wiemer, 2010; Schorlemmer et al., 2005). Schorlemmer et al. (2005) calculated the b‐value as a function of rake angle of events from regional and worldwide earthquake catalogues including the Northern California Seismic Network (NCSN), Southern California Seismic Network (SCSN), and the Harvard catalogue, and found that in all catalogues normal faulting corresponds with higher b values (i.e., distribution weighted towards smaller events), strike-slip faulting with intermediate values, and thrust faulting with smaller b-values. For example, the range of b-values given by the Harvard catalogue is 1.12 to 1.17 for normal fault stress regimes, 0.98 to 1.13 for strike-slip, and 0.81 to 0.93 for thrust fault stress regimes (Schorlemmer et al., 2005). The range of b-values can also be seen to vary between the different catalogues, for example in the Harvard catalogue, the range for thrust regimes is 0.81 to 0.93, whereas for the SCSN the range is 0.72-0.79, and for the NCSN it is 0.67-0.73 (Schorlemmer et al., 2005). 35 Complementing these, numerical modelling investigations have been widely used to study induced seismicity mechanisms. For example, continuum modelling techniques (finite-element, finite-difference, etc.) have been used to study the effects of fault permeability on the magnitude of seismic events (Rutqvist et al., 2013), showing that reactivation of a fault can result in a permeability change of several orders of magnitude affecting the final rupture length and seismic magnitude. A parallel study using similar techniques showed that fault rupture lengths are limited to a few hundred meters due to permeability increases and leak off relative to the fluid volume injected (Rutqvist et al., 2015). Employing discontinuum modelling techniques (distinct-element, discrete-element), where the fault and associated fracture network can be explicitly modelled, Zangeneh et al. ( 2013) were able to use moment magnitude and seismic moment relationships to calculate the event magnitude for a simulated induced seismicity event from the modelled shear slip displacement distribution along a critically stressed fault. It was shown that a felt event can be generated in response to hydraulic injection near the critically-stressed fault, even several hours after injection is stopped due to the diffusion of the pore pressure front eventually reaching the fault and reducing the effective stresses. Pirayehgar and Dusseault (2015) built on this approach to show that the magnitude of induced seismic events changes with the rock mass fabric characteristics and in situ stress ratio. 2.3 Empirical Database Study 2.3.1 Methodology and Databases Compiled From the data available in the World Stress Map database (Heidbach et al., 2018), it can be seen that although local variations exist, certain regional generalizations can be made regarding the tectonic stress regime across North America (Figure 2.1). For example, the U.S. west coast south 36 of the Cascadian subduction zone is generally dominated by a strike-slip regime typified by the San Andreas fault. This transitions to normal faulting and strike-slip regimes in the midcontinent (e.g., Basin and Range), and finally to a thrust fault regime in the east. To the north, western Canada along the Rockies is dominated by a thrust faulting regime. This suggests that shale gas basins in different parts of North America are subject to different tectonic stress regimes and therefore different faulting mechanisms, which raises the question of the influence of the stress regime on induced seismicity magnitudes and magnitude distributions. Figure 2.1: Stress map of North America using data from the World Stress Map (Heidbach et al., 2018). 37 To investigate this, three databases were compiled, cross-correlated and analyzed: 1) Stress Regime Database: A database was developed in order to determine the general stress regime representative for different North American shale gas basins. Data from focal mechanism solutions for 2996 earthquakes was used and combined with other stress state studies reported in the literature (Heidbach et al., 2018; Herrmann et al., 2011; H Kao et al., 2012; Lund Snee & Zoback, 2016; McNamara et al., 2015; USGS, 2018). These are superimposed on top of a map of different North American shale basins in Figure 2.2. The data suggests that the stress regime in northeastern British Columbia where the Montney shale gas play is being developed is predominantly a thrust faulting regime. For comparison, the data shows that Oklahoma where the Woodford shale gas play is being developed is predominantly strike-slip. The Montney and Woodford are selected for comparison in this chapter as representing two of the more active shale gas plays in terms of induced seismicity concerns. It should be noted that we’re considering the stress regime on a regional (shale gas basin) scale and that local stress field variations are not considered here. 2) Well Activity Database: Digital data was compiled using publically available sources reporting hydraulic fracturing and waste water disposal well activities in British Columbia (BC Oil & Gas Commission, 2018) and across the USA (Fracfocous, 2017). The resulting database consists of records for 23,828 wells since 2013 in northeastern British Columbia and 120,247 wells since 2009 in the USA. 3) Located Earthquake Database: A database was developed from publically available earthquake catalogues of located earthquakes in British Columbia and across the USA. It should be noted that these do not specify whether the earthquakes are natural tectonic 38 events or triggered by anthropogenic activities. An initial step that was carried out when exporting the data from the earthquake databases was to filter out non-oil and gas related anthropogenic events such as those from mining and blasting activities. These included attributes in the metadata that made them distinguishable. This was not the case for natural earthquakes and induced seismicity related to oil and gas activities, which were filtered differently as discussed below. In total, records for 891 earthquakes located in northeastern British Columbia after 2013 were compiled (Earthquakes Canada, 2018) together with 206,304 earthquakes recorded by the USGS across the USA since 2009 (USGS, 2018). Figure 2.2: Compilation of reported focal mechanisms for North American earthquakes since 2009 used here to indicate the general tectonic stress regimes for major North American shale basins outlined in green. 39 The compilation of these databases was required given the absence of any comprehensive, publically-available records of induced seismicity events already correlated to hydraulic fracturing activities. As previously noted, earthquake catalogues only report the occurrence of a detected event without making any interpretation as to its cause (natural or anthropogenic). In this study, I attempt to do so by cross-correlating the databases compiled and applying a series of spatial and temporal filters to identify the subset of earthquake events that are likely induced seismicity events related to hydraulic fracturing and waste water disposal activities for natural gas development. For the spatial filter, a 10 km radius was applied to search for all earthquake locations that were within 10 km of an active hydraulic fracturing or wastewater well. In this case, 10 km represents the assumed uncertainty in the event location accuracy. To this, a 3-month temporal filter was applied (Atkinson et al., 2016), which accounts for uncertainty in the well activity records as to when the hydraulic fracturing or waste water injection treatment was actually completed (the records only report the start date). Thus, if an earthquake event was recorded as occurring within 10 km of an active well and within 3 months from the start date of the hydraulic fracturing or waste water disposal activity, it was considered here to be an induced seismicity event. Figure 2.3 shows the identified induced seismicity events we ascribe to hydraulic fracturing and waste water injection activities. 40 Figure 2.3: Induced seismicity events and their magnitudes related to hydraulic fracturing and waste water injection, identified by cross-correlating the compiled databases and applying a 10-km spatial filter and 3-month temporal filter relative to well locations and activities, respectively. These are superimposed over a map of major North American shale basins outlined in green. After applying this filter, 359 events were identified in the Montney, and from the 11,221 events located in the USA, 385 were identified in the Woodford. These results are shown in Table 2.1. 41 Table 2.1: Breakdown of identified induced seismicity events for different major shale gas plays in North America. Shale Play Magnitude (ML) Eagle ford Fayetteville Barnett Woodford Pierre Niobrarra Monterey Marcellus Montney (Canada) 1-2 0 348 0 115 0 6 51 5 148 2-3 4 673 26 251 21 6 25 8 193 3-4 5 32 12 19 17 49 3 0 15 4-5 1 2 0 0 3 0 0 0 3 5-6 0 0 0 0 1 0 0 0 0 Total 10 1055 38 385 42 61 79 13 359 42 2.3.2 Empirical Analysis and Results Based on the number of identified induced events shown in Table 2.1, the Montney and Woodford basins were selected for further analysis. It is noted that the number of identified events in other important basins such as the Barnett and Marcellus were too few to allow the b-value to be calculated without significant uncertainty. The Montney Trend is a 26,000-square-kilometre siltstone formation that stretches from the B.C./Alberta border near Dawson Creek to the B.C. Rocky Mountain foothills 200 kilometers northwest. Its depth ranges from 1,400 to 3,800 meters below surface and its thickness varies from 30 to 300 meters. The majority of the Montney is made up of extremely low permeability highly laminated organic claystone, siltstone and fine sandstones with porosities ranging from 2 to 9 %. The depth of the formation also increases from northeast to southwest, generally along with increasing reservoir pressures from 15 to 50 MPa and decreasing natural gas liquid (NGL) and oil content. Thus, reservoir characteristics vary widely across the formation (BC Oil and Gas Commission, 2012). The Devonian-Mississippian Woodford Shale is an important petroleum source rock and reservoir in the Anadarko Basin Province. This basin extends into Kansas, Oklahoma and Texas. The most widespread and characteristic Woodford Basin lithology is black shale. Other common lithologies include chert, siltstone, sandstone, dolostone and light-colored shale, with hybrid mixtures between them. Porosities are reported to range from 3 to 9% (Soeder, 1988). The depth of the Woodford Basin ranges from 1,800 m on the southern Kansas shelf to 3,400 m in the deep 43 basin of southern Oklahoma. The net thickness of the Woodford Shale ranges from 35 to 70 m (Arthur & Langhus, 2008). The induced seismicity events identified for the Montney and Woodford were next analyzed by plotting their frequency-magnitude distributions and applying a linear fit to calculate their respective b-values. These are shown in Figure 2.4. From this, the b-value for the Montney, which was characterized in Figure 2.2 as a thrust faulting stress regime, was calculated to be 0.99. In contrast, the b-value calculated for the Woodford, which was characterized in Figure 2.2 as being a strike-slip regime, is significantly higher at 1.46. These results suggest that thrust regimes are more likely to produce higher magnitude induced seismicity events than strike-slip regimes. Figure 2.4: Magnitude distribution and calculated b-values for induced seismicity events in the Montney versus the Woodford. The maximum curvature method (Wiemer & Wyss, 2000; Woessner & Wiemer, 2005) was used to calculate the magnitude of completeness of the induced seismicity events in the Montney and b=1.4600.511.522.530 1 2 3 4 5Log NMagnitudeWoodfordb=0.9900.511.522.530 1 2 3 4 5Log NMagnitudeMontney44 Woodford basins. The completeness magnitude was 2.0 for the Montney and 2.3 for the Woodford. It should also be noted that the calculated b-values do not include microseismicity events as these are below the detection thresholds represented in the databases. Thus, the calculated b-values relate more to larger events and are therefore characteristic of the general basin scale response than that at the detailed local scale specific to each fault. These findings are further explored by plotting histograms of the induced seismicity event distributions for the Montney and Woodford plays (Figure 2.5). Each histogram bin is 0.2 in magnitude. The kernel distribution fits of these histograms indicate that the majority of events for both plays involve magnitudes less than 2.4. Both trends will include some degree of censoring of low magnitude events, as noted above, related to the lower detection thresholds used by the different seismic monitoring systems being used. At the upper end, it can be seen that the Montney has experienced a small but significant number of larger magnitude events (> M 4) exceeding the largest event for the Woodford (M 3.6). Although the trends from these plots suggest that the average event magnitude for the Woodford is slightly larger than that for the Montney (M 2.2 versus 2.0, respectively), the Montney has clearly experienced larger events than the Woodford play. 45 Figure 2.5: Histogram of induced events in Montney and Woodford and their distribution curves. 2.3.3 Sensitivity of Results The distribution of b-values for these two plays was further investigated to see the effect of the spatial and temporal filter lengths on the b-value of identified events (Figure 2.6). The sensitivity of the spatial and temporal filters in identifying induced events was tested for each play by varying the spatial cutoff thresholds between 5 and 10 kilometers radius and the temporal cutoff between one week and six months. The results are presented in Table 2.2. 46 Table 2.2: Sensitivity of the spatial and temporal filter used on the b-value of induced events in the Montney and Woodford shale gas plays. b-value (Montney) b-value (Woodford) 5 km radius 10 km radius 5 km radius 10 km radius 1 Week 0.97 0.92 1.02 1.52 2 Weeks 0.99 0.94 1.13 1.68 1 Month 0.92 0.97 1.35 1.51 2 Months 0.88 0.95 1.37 1.42 3 Months 0.91 0.96 1.34 1.46 6 Months 0.95 0.96 1.45 1.59 Mean 0.94 0.95 1.28 1.53 SD* 0.041 0.018 0.165 0.093 *Standard Deviation The results of the sensitivity analysis consistently show a lower b-value for the Montney in all cases relative to the Woodford (i.e., higher probability of larger induced events in the Montney than the Woodford). As seen from Table 2.2, the b-values for the Woodford play are lowest for a one and two week temporal filter combined with a 5 kilometer radius spatial filter. This can be related to the fact that with smaller filtering values, fewer events will be identified as induced events, in this case 124 and 170 events, which appears to censor lower magnitude events. Overall, the results of this empirical study suggest that a thrust faulting stress regime is more prone to larger induced events than a strike-slip stress regime. This finding was further 47 investigated by numerical simulations of induced seismicity as presented and discussed in the next section. 2.4 Numerical Modelling 2.4.1 Modelling Method and Setup The influence of stress regime on induced seismicity was further investigated by carrying out a series of numerical simulations using the commercial 3-D distinct element code 3DEC (Itasca Consulting Group Inc., 1999). The distinct-element method models a problem domain (e.g., a rock mass) as an assemblage of deformable blocks defined by a network of discontinuities of variable orientation, spacing, and persistence. The discontinuities represent interfaces that may undergo large deformations in shear and/or opening in response to changes in the modelled effective stresses. 3DEC is hydro-mechanically coupled and capable of simulating fluid flow along the discontinuity network (i.e., fracture flow) via a cubic law aperture relationship. This relates mechanical deformation occurring in the form of normal closure or opening of a discontinuity to changes in its hydraulic conductivity and the subsequent distribution of fluid pressures; conversely, changing fluid pressures along a discontinuity will result in a corresponding change in mechanical aperture, as well as in the effective stresses acting along the discontinuity thereby creating the potential for slip. The interfaces between blocks can be bonded to represent an incipient fracture that may initiate, propagate, and open in response to an injected fluid pressure. This makes the distinct element method and the 3DEC code highly appealing for modelling hydraulic fracturing and induced seismicity (Zangeneh et al., 2013). The workflow of the 3DEC modelling performed for this study is shown in Figure 2.6. 48 Figure 2.6: General modelling procedure used for 3DEC simulations. The 3DEC analysis was carried out to simulate fault slip under different stress regimes while keeping all other parameters constant. The model consists of a block with 400 m edge lengths, divided by an inclined fault and network of parallel and orthogonal discontinuities (Figure 2.7). The intact rock blocks were modeled as elastic with a density of 2,700 kg/m3 and bulk and shear moduli of 20 and 12 GPa, respectively. The discontinuities were modelled as incipient fractures with a shear strength representative of intact sedimentary rock, but that would drop to residual values if the strength is exceeded and a fracture initiated. The fault was modelled as a discrete, weak interface. Both the incipient fractures and fault were modelled assuming a Coulomb slip relationship and their properties are provided in Table 2.3. 49 The vertical stresses were assigned assuming a depth of 1,300 m at the top of the model and 1,700 m at the bottom. It should be noted that these depths correspond to the shallower limits of hydraulic fracturing operations in Montney. However, for the purpose of the study in this chapter, this does not influence the results of numerical simulations. The horizontal stresses were varied for each model according to the tectonic stress regime being simulated. These are reported in Table 2.4 with respect to the maximum horizontal to vertical stress ratio, K. Table 2.3: Discontinuity strength and deformation properties used in the 3DEC numerical simulations. Discontinuity Properties Incipient Joints Fault Friction angle (deg) 40 30 Cohesion (MPa) 10 0 Tensile strength (MPa) 0.5 0 Residual friction angle (deg) 25 15 Residual Cohesion (MPa) 0.1 0 Fracture Aperture – azero (m) 1.0E-4 1.0E-4 Fracture Aperture – ares (m) 4.0E-5 4.0E-5 50 Table 2.4: Horizontal to vertical in-situ stress ratios modelled for each tectonic stress regime. Stress Regime K_Ratio Thrust 1.6 Reverse 1.6 Strike-Slip 1.1 Normal 0.625 The fluid injection point was located at 1,500 m depth, 45 m above the fault. The injection rate was 1 m3/sec for a duration of 10 minutes. Monitoring points were included across the fault surface to track stress drop in response to slip, which together with the area of slip, were used to calculate the event magnitude associated with the induced fault slip. 51 Figure 2.7: 3DEC distinct-element model developed for the numerical simulations, representing a rock mass block with edge lengths of 400 m, cut by a persistent fault (blue plane). The injection point is indicated in yellow and is located 45 meters above and normal to the fault (see left plot). The location of the monitoring points on the fault are shown in the right plot. 2.4.2 Modelling Results The results of the simulation are provided in Figure 2.8 and show the shear displacement magnitudes that occur on the fault in response to the fluid injection-induced slip. In terms of slip displacement magnitudes, the largest values were seen to occur for the normal fault stress scenario, but these were highly localized and occurred over the smallest fault slip area (Figure 2.8d). In contrast, it can be seen that the thrust and strike slip faulting scenarios have the largest slip area. Although these are comparable in area, the shear displacement magnitudes for the thrust fault case are much larger than those for the strike-slip fault case. 52 Figure 2.8: Results of numerical simulations for injection induced fault slip under four different stress regimes. Shear Displacement Magnitude (m)Shear Displacement Magnitude (m)Shear Displacement Magnitude (m)Shear Displacement Magnitude (m)a) Thrust faultd) Normal faultc) Reverse faultb) Strike-slip fault53 Using the measured slip areas and average shear displacements the seismic moment of these events was calculated (Hanks & Kanamori, 1979). This was further used to calculate the moment magnitude of these events. The resulting magnitudes were 2.2 for the thrust fault, 1.6 for the strike slip, 1.4 for the reverse, and 1.1 for normal fault scenarios. These results agree with the trends of the magnitude distributions calculated from the empirical data. It is important to emphasize that these models assume the same fluid injection depth for each case and therefore were initialized with the same vertical stresses acting on each fault. Where they differ is with respect to the initial horizontal stress imposed, with each fault (i.e., stress regime) scenario representing a different K-ratio (Table 2.4). Thus, when comparing the shear displacement magnitudes in Figure 2.8, it is recognized that the maximum horizontal stress acting on the thrust fault is larger than that acting on the normal fault. This results in a higher elastic strain energy being stored (and then released) for the thrust fault stress regime case and therefore larger stress drop and slip area. The shear stress drop for each scenario is shown in Figure 2.9. This shows that the stress drop observed for the thrust and reverse fault scenarios is approximately twice as much as those for the normal and strike slip fault scenarios. 54 Figure 2.9: Modelled shear stress changes (i.e., stress drop) for representative monitoring points along the fault (see Figure 2.7) for each stress regime modelled. 2.5 Discussion The numerical simulations performed for this study represent a comparative analysis specific to the effects that different stress regimes and fault mechanisms may have on the magnitude of 55 induced seismicity events. The models, outside of the presence of the fault, assume homogeneous and isotropic rock mass conditions. This simplification of the geology serves the purpose of allowing the effects of varying the stress regime to be more directly compared. Similarly, the effect of other parameters such as injection volume and rate are not considered here. Another factor that needs to be addressed is the scale at which both the stress regime and seismicity are investigated. The focus in this study was the influence of the tectonic stress regime and therefore the global behavior of the seismicity at the regional scale. It is also recognized that stress effects can vary at the local scale (Cloos, 1928; Riedel, 1929). The same can be said for seismicity, and so it is recognized that the results presented here focus on large induced events above the magnitude of completeness of the regional seismic network. The calculated b-values will vary if the microseismic events representing the local scale response are included. In addition to the tectonic stress regime, several other factors can also influence the magnitudes of induced seismicity events, such as the local geology, depth of fluid injection, and differences in operational parameters (e.g., fluid injection volumes). Therefore, these too might also partly explain the differences in the event distributions observed in the data between the Montney and Woodford reported in the previous sections. With respect to geology, both the Montney and Woodford reservoir rocks are similar in that they both consist of very fine sedimentary rocks with very low permeabilities. 56 One key difference, though, is related to operational differences and the sources of induced seismicity. For the Montney, the induced seismicity is associated with hydraulic fracturing operations and the interactions between the injected fluids and faults that cross the invaded zone and adjacent formations (BC Oil and Gas Commission, 2014). Little induced seismicity in the Montney has been associated with subsequent wastewater disposal activities, which typically target shallower formations. In contrast, the majority of induced seismicity in the Woodford is related to waste water injection in the highly permeable Arbuckle group at depths well below the formations being targeted by hydraulic fracturing (Langenbruch & Zoback, 2016). This involves interactions between faults in the crystalline basement and pressurized fluids migrating downwards from the disposal formations (Hincks et al., 2018; Walsh & Zoback, 2015). This raises the question of differences in operating depths between the Montney and Woodford, but perhaps similarities with respect to the larger induced seismicity events being associated with deeper wells. Recognizing that the vertical location accuracy of induced seismicity events is generally poor, a histogram of the depths of the wells associated with induced seismicity for the Montney and Woodford was generated and is shown in Figure 2.10. This shows that the associated wells for the Montney have a distribution between 3,400 to 4,800 meters. This corresponds with the depths of the Montney formation and the Belloy/Debolt formations below it. Since the formations targeted for hydraulic fracturing are not significantly thick (BC Oil and Gas Commission, 2014), it can be assumed that the depths of induced seismicity are approximately the same as the injection depths. In contrast, the depths of the wells in the Woodford associated with induced seismicity show a broader distribution of 1,200 to 4,800 meters, as would be expected where events are being generated for both the shallower hydraulic 57 fracturing operations and deeper waste water injection close to crystalline basement. Hincks et al. (2018) determined the depth of the basement interface across Oklahoma, which showed that the depth of the basement increases in the Anadarko and Arkoma basins from 3,000 to 5,500 meters. Thus, the depths of induced events related to wastewater injection in the Woodford appear to be similar to those for hydraulic fracturing in the Montney. Figure 2.10: Histograms comparing the depth of wells in the Montney and Woodford. Next, a comparison can be made between the injection volumes of the wells being considered. Data for the Montney (BC Oil & Gas Commission, 2018) and Woodford (OCC, 2018), reported in Table 2.5, shows that the scale of operations are much higher in the Woodford with twenty times more active wells (i.e., wells reporting annual injection volume > 0 m3). As a result, Table 2.5 also shows that the annual total volumes being injected are two orders of magnitude higher in the Woodford compared to the Montney. However, with respect to the average monthly injection 58 volumes per active well (Figure 2.11), the volumes are in a similar range, with the Woodford being slightly higher likely due to the heavier weighting towards waste water injection wells. Table 2.5: Annual total injection volumes reported for the Woodford and Montney. Woodford Montney Year Injection Volume (m3) Active wells Injection Volume (m3) Active wells 2013 4.35E+08 8531 4.49E+06 353 2014 4.96E+08 8952 8.12E+06 648 2015 4.88E+08 8966 7.66E+06 530 2016 4.35E+08 8757 4.56E+06 309 Figure 2.11: Average monthly injection volumes per well (solid lines) and overall averages (dashed lines) for the Montney and Woodford basins. 59 Based on these data and comparisons, the differences in the magnitude distribution of induced seismicity events observed in the empirical analysis is likely most strongly influenced by the differences in the tectonic stress regimes between the Montney and Woodford. More specific to the Montney, it is also understood that the stress regime is not necessarily uniform throughout the basin and that the presence of certain regional-scale geological structures can contribute to stress heterogeneity. An example of one such regional-scale structure is the Fort St. John Graben (FSJG). The FSJG is a west-trending fault system located at the center of the Peace River Embayment of the Western Canada Sedimentary Basin (Figure 2.12). The FSJG forms the central arm of a larger fault system known as the Dawson Creek Graben Complex which was filled with sediments in several stages during the Carboniferous to Permian periods (Barclay et al., 1990). Although the Dawson Creek Graben Complex was formed during a long period of extension, its bounding fault systems often display alternating zones of both compressional and extensional structural styles. The unique structural characteristics of the graben complex have been attributed to early strike-slip faulting (Carboniferous to Permian) as well as to late reactivation and inversion processes associated with the development of the nearby Rocky Mountain fold and thrust belt (Berger, 1994). 60 Figure 2.12: Sketch of Carboniferous-Permian Dawson Creek Graben Complex (Barclay et al., 1990). The presence of the FSJG in the Montney separates the stress regime in this basin into a thrust fault stress regime in the northern part and a strike slip stress regime in the south (Figure 2.13). Included with this interpretive figure are the focal mechanisms determined for several natural earthquakes in the region. This shows that most of the events in the northwestern and central Montney are thrust faulting earthquakes and those in the southeastern parts of the basin represent a change to strike slip events. 61 Figure 2.13: Focal mechanisms of earthquakes in the Montney basin relative to the axis of the Fort St. John Graben (FSJG). The location and magnitude of the induced seismicity events related to hydraulic fracturing are shown in Figure 2.14. These events are obtained from the compiled database previously presented. As can be seen in this figure, the magnitudes of the induced seismicity events in the thrust stress regime portion of the Montney are higher than those in the strike slip portion. This follows the same trend of induced seismicity magnitude distributions seen in the larger dataset comparing the overall Montney (a higher upper bound due to the presence of a thrust fault stress regime) with the Woodford (a marginally smaller upper bound due to the predominantly strike-slip stress regime). 62 Figure 2.14: Location and magnitude of induced seismicity events in the Montney relative to the Fort St. John Graben (FSJG). 2.6 Conclusions An empirical study supported by 3-D numerical modeling was performed to investigate the influence of different tectonic stress regimes on the magnitude and magnitude distribution of induced seismicity events related to hydraulic fracturing activities. This included the compilation of existing stress, seismicity and well data for key North America shale gas plays. The stress data showed that the stress regime for the Montney play in northeastern British Columbia is predominantly a thrust fault regime (i.e., high horizontal stress relative to the vertical stress with 63 depth). In contrast, the stress regime for the Woodford play in Oklahoma was seen to be that of a strike-slip fault regime (i.e., a horizontal stress similar to the vertical stress with depth). For the seismicity and well data, a spatial and temporal filter was applied to identify those wells associated with induced seismicity events for each play. The resulting frequency-magnitude plots showed a lower b-value for the Montney (i.e., skewed towards larger magnitude events) than the b-value for the Woodford. These findings were supported by the results obtained for the 3-D numerical simulations. The results showed that under similar geological and operational conditions, the tectonic in situ stress regime has a significant influence on the amount of strain energy that can be stored and therefore released through the area and magnitude of slip on the fault. Both the area and magnitude of slip influence the magnitude of the induced seismicity event. Results from the numerical simulations showed that both the magnitude and area of slip are greater in the case of a thrust fault stress regime, compared to the other fault/stress regimes if all other factors are kept equal, owing to the high horizontal stresses for this case and how the stress field resolves to act on the fault. This combination acts to promote the stored strain energy being released as a single large event. In contrast, numerical simulations for the other faulting types showed these to be more susceptible to either smaller slip magnitudes or smaller slip areas with the stored strain being released as multiple smaller slip events. Accordingly, the modelling results show a higher upper bound can be expected in the distribution of induced seismicity events when operating in a thrust fault stress regime compared to a strike-slip or normal fault stress regime. 64 Chapter 3: Empirical and Numerical Investigation on the Influence of Fluid Injection Volume and Rate on Induced Seismicity in the Montney Formation, Northeastern British Columbia 3.1 Introduction In recent years, certain regions across Canada and the United States have experienced a significant increase in seismicity relative to historical baselines (Ellsworth, 2013; Farahbod et al., 2015a,b; Keranen et al., 2014). This increase has been linked to hydraulic fracturing and deep wastewater disposal wells associated with the development of unconventional gas resources (BC Oil and Gas Commision, 2012; BC Oil and Gas Commission, 2014; Horton, 2012a; Tagliabue, 2013). The hydraulic fracturing process involves pumping fluids under high pressure into isolated sections of a wellbore to generate fractures in order to increase the permeability and stimulated volume of the reservoir. At the same time, flowback fluids and produced water arising from these operations (i.e., wastewater) requires disposal, which is usually carried out by injection into deep formations. In both cases, if fluid flow connectivity is established proximal to a critically stressed fault, the fluid injection serves to create localized increases in pore pressures and reductions in the effective normal stresses acting on the fault, resulting in fault slip and induced seismicity. For the most part, these induced seismic events have low magnitudes (M3) that in cases involving higher population densities and sensitive ground conditions have resulted in damage to infrastructure and property (Tagliabue, 2013). As a result, induced seismicity has heightened public and regulator concerns in affected regions, leading to a number of scientific 65 reviews on hydraulic fracturing pointing to the need to study and better understand the factors that influence induced seismicity and event magnitudes (Allen et al., 2019). Induced seismicity is a complex interaction of geomechanics and anthropogenic activities. It is highly unlikely that a single causative factor is solely responsible for an induced seismicity event. Instead, multiple factors, both geological and operational (i.e., completion related), can play an influencing role. Geology and regional tectonics are of primary importance and are the subject of several recent research studies (Amini & Eberhardt, 2019; Göbel, 2015; Kao et al., 2018; Pawley et al., 2018; Shah & Keller, 2017). Proximity to basement is another geological factor shown to have a controlling influence (Currie et al., 2018; Pawley et al., 2018; Skoumal et al., 2015). These factors certainly influence the susceptibility of a targeted formation but have not yet been controlled or manipulated (outside of avoidance). Operational factors on the other hand, such as fluid injection volumes and rates, may be controlled, though may also reduce well production, possibly offers a means to minimize induced seismicity hazard for susceptible formations. In this study, focus is placed specifically on operational factors. In the United States, most reported injection induced earthquakes are associated with wastewater disposal into deep formations beneath the production zones (Ellsworth, 2013; Hincks et al., 2018; Keranen et al., 2014), and some studies have shown that injection volume is associates with triggering induced seismicity. For example, Horton (2012) showed in a study of wastewater disposal in Arkansas a high correlation between injection volume and seismicity, together with time lags of up to 15 weeks between injection and seismicity. Kim (2013) investigating induced seismicity associated with a deep fluid injection well in Ohio, found correlations between periods 66 of seismic inactivity and minima in injection volumes. McGarr (2014) in his global analysis of case histories of earthquake sequences triggered by fluid injection found that maximum event magnitudes appear to be limited according to the total volume of fluid injected. In contrast, other studies have pointed to the importance of injection rate. Weingarten et al. (2015) in their study of empirical data from the eastern and central U.S. concluded that injection rate is the most decisive operational factor affecting the triggering of an induced seismicity event for disposal wells. Furthermore, they found that the cumulative injected volume, monthly wellhead pressure, depth, and proximity to the crystalline basement did not strongly correlate with earthquake association. Curry et al. (2018) compared observed seismicity rates and monthly injection records from disposal wells in Ohio in which they found a close correspondence between earthquake frequency and injection rate. Kwiatek et al. (2015) investigating induced seismicity at the northwestern Geysers Geothermal field in California showed that changes in long-term spatial, temporal, and seismic source characteristics could be clearly attributed to variations in injection rates. In the Western Canadian Sedimentary Basin, the majority of induced seismicity events appear to correlate with fluid injection for hydraulic fracturing and not wastewater disposal (e.g., Mahani et al., 2017; Atkinson et al., 2016; Farahbod et al., 2015 ). Schultz et al. (2018) investigated the relationship between injection parameters and induced earthquakes in the Duvernay shale play in Alberta and concluded that induced earthquakes are associated with hydraulic fracturing operations that used larger injection volumes and that seismic productivity scales linearly with injection volume. Their analysis also showed that injection pressure and rate do not have a significant association with seismic response. Kao et al. (2018) in their study of induced 67 seismicity in western Canada showed that the maximum magnitude appears to increase with the total injected volume, but large injection volumes do not always result in large magnitude events. Although the above studies draw correlations between induced seismicity and operational factors such as injection volume and rate, experiences vary between different shale gas plays regarding which parameter(s) is more influential. These variations are mostly due to the differences in the respective regional and local geology as well as data quality and quantity, but at the same time they present knowledge gaps and uncertainties regarding the influence of these parameters on triggering induced seismicity and their magnitude distributions. In this work, we investigate the association of fluid injection volume and rate with induced seismicity in the Montney shale gas play in northeastern British Columbia (NEBC). To do so, we have compiled an empirical database of induced seismicity events and event magnitude distributions together with well activity and injection volumes and rates. This chapter presents the results from the analyses of these data together with supporting results from 3-D numerical modeling used to aid interpretation of the empirical relationships derived. 3.2 Geological Setting of Montney The Lower Triassic Montney Formation is aerially extensive, covering approximately 130,000 km2 from central Alberta to NEBC (National Energy Board, 2013). It is also thick, typically ranging from 100 to 300 m, but thinning to zero at its eastern and northeastern margins while increasing to over 300 m on its western margin before it begins outcropping in the Rocky Mountains. The Montney Formation unconformably overlies Carboniferous or Permian strata and consists of variable amounts of interbedded shale, siltstone, and sandstone. These strata 68 developed during the first of three major transgression-regression cycles (Edwards et al., 1994; Gibson & Barclay, 1989). Based on the lithostratigraphy, Davies et al. (2018) subdivides the Montney Formation into three members: the Upper, Middle, and Lower Montney. These members are separated by a basin-wide unconformity that developed due to tectonic uplift of the basin margin (Dixon, 2009). 3.3 Empirical Study 3.3.1 Data Compilation and Preparation The database used as the starting point for this study is a comprehensive earthquake catalogue compiled for NEBC and western Alberta (Visser et al., 2017). This catalogue was produced specifically to study induced seismicity in this region and consists of 4916 events for the period of January 2014 to December 2016 with a magnitude of completeness of 1.8 (ML). It should be noted that this catalogue does not specify whether each earthquake is a natural tectonic event or was triggered by anthropogenic activities. A second database was compiled using publicly available sources reporting hydraulic fracturing well activities in NEBC (BC Oil & Gas Commission, 2018). This included parameters such as the well location, targeted injection depth(s), and dates for different well activities, as well as other operational parameters such as the injection rate and volume. Figure 1.1 provides a location map of these wells together with the locations of the seismic events in the earthquake catalogue for the Montney region. As previously noted, the earthquake catalogue only reports the occurrence of a detected event without making any interpretation as to its cause (i.e., natural or anthropogenic). To separate these, we cross-correlated the earthquake catalogue with the well database and applied a series of spatial and temporal filters to identify the subset of earthquake events that are likely induced 69 seismicity events related to hydraulic fracturing and waste water disposal activities. The first step was to clip the data to only include earthquakes located within the boundaries of the Montney area and to filter out events spatially associated with non-oil and gas related anthropogenic activities such as those from mining and construction (i.e., blasting). This step reduced the total number of events being considered from 4916 to 2867. Next, a spatial filter was applied to search for all event locations that were within a 5 km radius of an active hydraulic fracturing well. The 5 km radius used represents the uncertainty in the event location accuracy reported for the earthquake catalogue. To this, a one-week temporal filter was applied (based on feedback received from BC Oil & Gas Commission). Thus, if an earthquake event was recorded as occurring within 5 km of an active well and within 1 week from the start date of the hydraulic fracturing activity, it was considered here to be an induced seismicity event and the corresponding well as being associated with induced seismicity. This resulted in a subset of 507 events identified as induced seismicity events. Next, we analyzed the relation between induced seismicity in the Montney region and the fluid injection volumes used during hydraulic fracturing operations. Figure 3.2 shows the raw temporal distribution of hydraulic fracturing injection volumes over the three-year period being considered, superimposed with the induced seismicity events recorded. The latter are plotted according to their event magnitudes. Figure 3.3 plots the cumulative induced seismicity count against the cumulative injection volume for those wells assessed as being associated with induced seismicity (Figure 3.3). The results show a good linear regression fit, indicating that injection volume is an important contributing factor for induced seismicity in NEBC. 70 Figure 3.1: Spatial distribution of earthquakes (red circles) and hydraulic fracturing wells (black circles) included in the respective databases compiled for the Montney play (black polygon) in northeastern British Columbia. 71 Figure 3.2: Timing and magnitude of induced seismicity events (red circles) superimposed over the daily fluid injection volumes for hydraulic fracturing wells in the Montney (black bars). Figure 3.3: Number of induced seismicity events versus cumulative injection volume for hydraulic fracturing wells associated with induced seismicity in the Montney. Compared are the data (red points) and the linear regression fit (blue line). The goodness-of-fit is reported with respect to the R2 value (i.e. coefficient of determination). 72 3.3.2 Statistical Analysis Using the subset of seismogenic wells determined, we then carried out a statistical analysis to compare the distributions of the hydraulic fracturing injection parameters associated with these wells to those for the parent distributions involving all wells in the Montney. The statistical distributions were analyzed using the Kolmogorov-Smirnov (KS) test (Massey, 1951) to compare whether the respective distributions are different. In these tests, we used the standard significance level of 0.05. The distributions for the fluid average injection rate, injection volume and maximum injection pressure per well, together with the computed p values, are reported in Figure 3.4. These indicate a rejection of the null hypothesis in both cases of injection rate and volume (PKS<0.05), meaning that the distribution of injection rates and volumes for the seismogenic wells are not a random subset and have a different distribution than the parent set comprised of all wells in the study area. This was not the case for maximum injection pressure for which the calculated p value was higher than 0.05. Next, Mann-Whitney (MW) U test (Gibbons & Subhabrata, 2011) was applied to determine whether the subset and parent distributions have statistically different median values. The results again point to a rejection of the null hypothesis (PMW<0.05) for injection rate and volume, indicating that in both cases, the respective medians for the distribution of injection rates and injection volumes are different between the subset of seismogenic wells and parent distribution encompassing all wells. Again, for maximum injection pressure the MW test fails to reject the null hypothesis and the results of these two tests indicate that maximum injection pressure is not associated with induced seismicity in this region. 73 Figure 3.4: Histograms showing the distribution of injection rates (left), injection volumes (middle) and injection pressure (right) for those wells associated with induced seismicity events (shaded red) compared to the parent data set encompassing all hydraulic fracturing wells in the Montney (shaded grey). Figure 3.5 provides the results from a similar analysis, in this case comparing the histograms for the subset of seismogenic wells against the remaining subset of non-seismogenic wells. The results of the statistical analysis show that the seismogenic wells have different medians for the distribution of injection rates and injection volumes than the non-seismogenic wells. In contrast, the results show that the injection well pressure is not statistically associated (PKS,MW>0.05) with seismogenic wells. In other words, the distribution of injection pressures is the same for wells that have experienced induced seismicity as those that have not, thus suggesting that it is not a significant factor in the triggering of induced seismicity. Thus, these statistical tests show that 74 induced seismicity in the Montney is associated with fluid injection rates and/or injection volumes but not injection pressure. Figure 3.5: Histograms showing the distribution of injection rates (left) and injection volumes (middle) and maximum injection pressure (right) in the Montney for those wells associated with induced seismicity events (shaded red) compared to those that are non-seismogenic (shaded grey). The dashed lines indicate the median for each distribution. To further refine these results, we next investigated induced seismicity in the Montney based on which part of this formation is being targeted: The Upper, Middle or Lower Montney. Based on an initial inspection of the data, the Lower Montney was excluded as it was targeted by very few wells and instead the vast majority of unconventional gas wells target the Upper and Middle Montney (BC Oil & Gas Commission, 2018). Table 3.1 shows the breakdown of the wells 75 associated with induced seismicity in the Montney separated by the Upper versus Middle Montney. The number of total wells that target the Upper Montney is almost twice that for the Middle Montney; however, the percentage of wells in the Middle Montney associated with induced seismicity (18%) is more than that of the Upper Montney (12%). A similar trend is observed when considering only large events that are M>3: a higher number of these events are associated with hydraulic fracturing operations targeting the Middle Montney (2.6% of Middle Montney wells) compared to those targeting the Upper Montney (0.7% of Upper Montney wells). These results are in agreement with current industry and regulator understanding that the Middle Montney is more susceptible to induced seismicity. Table 3.1: Associated wells with Induced Seismicity Separated by Hydraulic Fracturing Target Formation Total Wells Wells Associated with IS Percentage Wells Associated with M>3 events Percentage Upper Montney 892 106 12% 6 0.7% Middle Montney 458 83 18% 12 2.6% Overall 1350 189 14% 18 1.3% Figure 3.6 presents the results of the correlations between injection rates and volumes, on induced seismicity as a function of which part of the Montney formation was being targeted These show that injection volume is statistically associated with induced seismicity (PMW<0.05) for the Middle Montney but not Upper Montney (Figure 3.6a, b). Interesting are the results for injection rate. The analysis indicates that injection rate is not statistically associated with induced seismicity for wells targeting the Middle Montney (Figure 3.6d). However, it is associated for wells targeting the Upper Montney (Figure 3.6c). 76 Figure 3.6: Histogram of operational parameters, including distribution of injection volumes and rates, for the Upper and Lower Montney. Plotted are the data for the seismogenic wells (shaded red) and all wells (shaded grey). The P values for comparing these two distributions are inset [KS, Kolmogorov-Smirnov; MW, Mann Whitney]. The dashed lines show the respective average values for the seismogenic wells (red) and for all wells (black). The observed differences in the association of injection volume and rate to induced seismicity in the Upper and Middle Montney can be the result of several factors ranging from differences in the geology to specifics related to their geomechanical characteristics and properties. Of note are several studies that have shown the Montney formation to be a fractured and faulted reservoir (Davies et al., 2018; James, 2011; Ouenes et al., 2013; Rogers & McLellan, 2014). It is well known that one of the requirements for the activation mechanisms of fault slip and induced 77 seismicity is the presence of a pathway for fluid pressure diffusion and/or stress transfer from the injection source to the fault (Allen et al., 2019). Here, we propose that the observation in differences in association of operational parameters with induced seismicity between Upper and Middle Montney can be a result of differences in natural fracture network densities, with Upper Montney having a lower fracture density (or absence of natural fractures) than Middle Montney. No similar studies report on the presence or characteristics of natural fractures in the Upper Montney, although it is expected that these would be somewhat different from the Middle Montney given differences in the rock types present (i.e. stratigraphy), dolomitization, clay/mica contents, and turbidites, and thus the corresponding physical properties (see Davies et al., 2018). However, this is acknowledged to be an unknown. To further analyze how variations in the natural fracture network may influence the association of induced seismicity with fluid injection rates and volumes, we compare the lower and upper bounds of the injection rates and volumes for the seismogenic wells to those for all wells in the Upper and Middle Montney. Figure 3.7 shows the box plots for each, where the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively, and the central mark indicates the median. The whiskers extend to the most extreme data points not considered to be outliers. Examining Figure 3.7a, it can be seen that the value of the median injection volume for the seismogenic wells in the Middle Montney is higher compared to the Upper Montney. One explanation might be that the characteristics of the fracture network in the Middle Montney is more inclined to promote leak-off and a larger invaded zone around the hydraulic fracture, thus requiring a higher volume of injected fluids. Inspecting the plot for the injection rates (Figure 3.7b), it can be seen that the median value is slightly higher for the Upper Montney. It could be 78 argued that with the absence or lower intensity of more persistent fractures or a less interconnected fracture network, a higher injection rate is needed to achieve the required stimulated volume. At the very least, the fact that the narrow range of injection rates for the seismogenic wells is effectively the same as that for the population of all wells suggests that injection rate alone is not associated with induced seismicity in this formation and that the interplay with geology must also be considered. Figure 3.7: Box plot of injection volumes and rates for all hydraulic fracturing wells in the Montney as well as the subsets of seismogenic wells in both the Upper and Middle Montney. For each box, the bottom and top edges indicate the 25th and 75th percentiles, respectively, and the central mark indicates the median. The whiskers extend to the most extreme data points not considered outliers. 79 Thus, the presence, absence or differences in the extent and characteristics of the natural fracture networks in the Upper and Middle Montney can be used to infer a mechanistic explanation for differences observed in the results of the statistical analyses regarding the association of induced seismicity with fluid injection volumes and rates. In the absence or limited presence of a natural fracture network, the hydraulic fracturing treatments would generate a more limited invaded zone and restricted connectivity pathways for fluid flow and pore pressure increases between the injection source and any susceptible faults in the area. In this case, the triggering of an induced seismicity event would more likely require that the hydraulic fracture directly intersect with a critically stressed fault. However, with the presence of a more extensive or interconnected natural fracture network, the hydraulic fracturing treatment would generate a much larger invaded zone with more connectivity and fluid flow pathways for pore pressure diffusion. The increased leak-off would reduce the influence of injection rate on induced seismicity (as shown in Figure 3.6d), and require larger injection volumes to achieve the required stimulated volume (as shown in Figure 3.7a). Conversely, the larger invaded zone increases the likelihood of interacting with a critically stressed fault and does not require the hydraulic fracture to intersect the fault for it to see localized increases in pore pressures. Faults that are proximal will also see increases in pore pressures through the fracture network. Given that the Lower Montney is cited to be similarly fractured as the Middle Montney (Rogers & McLellan, 2014), this would explain observations that the majority of the induced seismicity in the Montney region occurs in the Middle and Lower Montney. In the case of the Lower Montney, where there is little hydraulic fracturing activity, the connectivity of natural fracture pathways with the Middle Montney, where there is significant hydraulic fracturing activity, would presumably make it possible for 80 the downward migration of injected fluids and pore pressure diffusion into a higher stress environment, increasing the likelihood of encountering a critically stressed fault. Other factors such as magnitude and orientation of the in situ stress field (Amini & Eberhardt, 2019; Brudzinski & Kozłowska, 2019; Holmgren et al., 2019; Mahani et al., 2019) and mechanical properties of the reservoir rock (Eaton & Schultz, 2018; Pawley et al., 2018; Schultz et al., 2016; Van der Baan & Calixto, 2017) can also have a controlling influence on the induced seismicity responses being compared. Since there is no information on the fault attributes that produced these events (orientation, size, distance from injection point, etc.), it is not possible to draw conclusions regarding specific events or to confirm that responses used to establish the trends identified are solely due to injection rates and volumes. Nevertheless, the number of data points used in the empirical analysis is significant, limiting bias and censoring, and providing reasonable explanations for the observed trends. 3.4 Numerical Modelling Study 3.4.1 Model Setup To support the findings from the empirical study and to further investigate the mechanistic influence of the injection rate and volume on induced seismicity in the presence of a fracture network, we performed a series of numerical simulations using the 3-D distinct element code 3DECTM (Itasca Consulting Group Inc., 1999). The distinct-element method allows for a rock mass to be modeled as an assemblage of deformable blocks defined by a network of interface contacts of variable orientation, spacing and persistence. These interfaces can be assigned behavior models and representative properties to simulate natural fractures that may undergo 81 shear slip and/or aperture opening in response to changes in the modeled effective stresses. Similarly, the interface contacts can be bonded to represent a potential fracture pathway (referred to here as an incipient fractures) that may initiate, propagate and open in response to an injected fluid pressure. The treatment of these interfaces in 3DEC is hydro-mechanically coupled. Fluid flow and pore pressures along the fracture network are calculated via a cubic law aperture relationship. This relates mechanical deformation occurring in the form of normal closure or opening of the contact with changes in its hydraulic conductivity and the subsequent distribution of fluid pressures; conversely, changing fluid pressures along an interface will result in a corresponding change in mechanical aperture, as well as in the effective stresses acting along the discontinuity thereby creating the potential for slip. Thus, the representation of the rock mass can include a more persistent cross-cutting interface representing a continuous fault together with a network of non-persistent interfaces representing either smaller-scale natural fractures and/or incipient fractures along which a simulated hydraulic fracture may develop. This makes the distinct element method and the 3DEC code highly suitable for modeling hydraulic fracturing (Zangeneh et al., 2015) and induced seismicity (Zangeneh et al., 2013). The 3DEC analysis performed for this study was carried out to simulate fault slip under different fluid injection rates and volumes, while keeping all other parameters constant, for model geometries with and without natural fractures. The size of the model was determined through preliminary modeling to enable the simulation of induced seismicity events with similar magnitudes to those in the Montney database. The resulting model used, consists of a block with the dimensions of 1600 m x 1600 m x 800 m edge lengths, cut by an inclined fault and network of parallel and orthogonal natural and incipient fractures (Figure 3.8). The rock blocks were 82 modeled as elastic with a density of 2700 kg/m3, Young’s modulus of 30 GPa and Poisson’s ratio of 0.25. These were based on average properties reported for the Montney based on geophysical logs and laboratory testing (Vishkai et al., 2017). The incipient fractures were modeled using intact rock strengths representative of weak sedimentary rock, applying peak values that subsequently dropped to residual values if the strength criterion is exceeded and a fracture initiated. These were further modeled such that fluid flow was not permitted along the incipient fracture until the contact failed in tension or shear (i.e., a fracture was formed in response to fluid injection). This was not the case for the representation of natural fractures in the model, when included. Natural fractures were modeled with the same properties as the incipient fractures but with zero cohesion and zero tensile strength, and fluid flow was enabled from the start. Lastly, the fault was modeled as a discrete, weak interface, also with peak and residual values to account for irregularities and asperities along the surface. This, together with the incipient and natural fractures, were modeled assuming a Coulomb slip constitutive relationship. The corresponding properties are provided in Table 3.2. Table 3.2: Discontinuity strength and deformation properties used in numerical simulations Interface Properties Incipient Fractures Natural Fractures Fault Friction angle (deg.) 40 40 30 Cohesion (MPa) 10 0 0 Tensile strength (MPa) 0.5 0 0 Residual friction angle (deg.) 25 25 15 Residual cohesion (MPa) 0.1 0 0 83 The boundary conditions for the model were selected to simulate a thrust fault stress regime, which is generally considered to be the dominant stress regime in the northern Montney NEBC region (Amini & Eberhardt, 2019; Mahani et al., 2019). The vertical stress (σV) was assumed to be aligned with the minimum principal stress (σ3), the maximum principal stress (σ1) was assumed to be horizontal in the fault dip direction with a stress ratio of 1.6σV, and the intermediate principal stress (σ2) was assumed to be horizontal in the fault strike direction with a stress ratio of 1.3σV. The vertical stress gradient was assigned assuming a fluid injection point for a horizontal well located at 2400 m depth. This point is 90 m above the fault. A series of different injection rates were modeled, each for a duration of 30 minutes. Monitoring points were included across the fault surface to track stress drop in response to slip, which together with the area of slip, were used to calculate the event magnitude associated with fault slip. 84 Figure 3.8: 3-D distinct element model developed for the numerical simulations, representing a rock mass with edge lengths of 1600 m, cut by a persistent thrust fault and a network of orthogonal incipient and natural fractures (lower diagram). The upper diagram shows a cut away highlighting the through going dipping fault (shaded red) and the location of the fluid injection point small (yellow sphere), which is located 90 meters above and normal to the fault surface. 85 3.4.2 Model Results Five different fluid injection rates were modeled for the 3DEC numerical hydraulic fracturing simulations: 0.6, 3, 6, 12 and 60 m3/min. The range of 3 to 12 m3/min corresponds with the lower and upper bounds of the injection rates typically used in the Montney (Figure 3.7). The added cases of 0.6 and 60 m3/min were included to test extremes. For each injection rate, two sets of models were analyzed both including a network of incipient fractures, but one interspersed with natural fractures, representative of the Middle Montney, and one without natural fractures, inferred to be representative of the Upper Montney. Figure 3.9 shows the results for the injection rate of 6 m3/min applied for a 30 minutes period. Both show the hydraulic fracture extending in the σ1 direction, with an invaded zone developing around it. The invaded zone in this case refers to the rock mass volume that experiences a pore pressure increase through fracture connectivity and conductivity, whether involving generated or natural fractures, and corresponding leak-off (this is discussed in more detail in the Discussion). For the case without fractures (Figure 3.9a), the fluid injection rate translates into a high localized zone of pore pressures developing around the injection point and within the primary hydraulic fracture as there is less opportunity for pressure leak-off. The leak-off that does occur involves adjoining secondary induced fractures generated along the incipient fractures in response to the build-up of high pore pressures in the primary hydraulic fracture. The generation of these fractures include the opening of incipient fractures in the σ3 and σ2 directions, with both propagating in the σ1 direction. In contrast, when a natural fracture network is included (Figure 3.9b), the same injection rate results in a longer hydraulic fracture enabled by the open natural 86 fractures along its path. The presence of the natural fracture network adjacent to the injection point and hydraulic fracture also makes it possible for more leak-off into the invaded zone. Figure 3.9: 3-D distinct-element results showing the simulated hydraulic fracture and invaded zone for the two modeled reservoir geometries: (a) without natural fractures and incipient fractures only (Upper Montney type case), and (b) with both natural and incipient fractures (Middle Montney type case). The fluid injection rate for both models was 6 m3/min, which is approximately the average injection rate used in the Montney for hydraulic fracturing operations. (c-d) Conceptual illustration of hydraulic fracturing, leak-off and secondary fracturing resulting in an invaded zone of elevated pore pressures (shaded grey) for scenarios cases (c) without and (d) with the presence of an interconnected natural fracture network. 87 The result is that the case without a natural fracture network has a smaller invaded zone (Figure 3.9a,c), which translates into a much smaller area of interaction with the dipping fault where the hydraulic fracture intersects it (Figure 3.10a). In response, both the area of slip and average magnitude of slip is quite small (Figure 3.10b). The magnitude of this event was calculated using the seismic moment which can be calculated using the measured slip area and average shear displacement across that slip area. The model results indicate that the slip event in this case produced a Mw = 2.5 event. In comparison, the model where a natural fracture network is present produced a much larger invaded zone, through leak-off into these fractures, for the same fluid injection rate and duration (Figure 3.9b,d). This larger invaded zone extended both horizontally with the hydraulic fracture to intersect the dipping fault but also vertically downward along the fluid flow pathways created by the natural fractures to intersect elsewhere with the fault. The net result is a much larger area on the fault surface seeing a pore pressure increase (Figure 3.10c). The corresponding slip event (Figure 3.10d) generated a Mw = 3.9 event. Thus, the effect of an interconnected natural fracture network can be seen to create multiple fluid flow pathways to reach the fault, which has two effects: 1) pore pressure magnitudes reaching the fault were higher; and 2) a larger invaded zone increases both the likelihood of interacting with a critically stressed fault and a larger area on the fault being affected by a pore pressure increase. The result is a larger induced seismicity event. Whereas when there the natural fracture network is more limited, the smaller invaded zone reduces the likelihood of intersecting a critically stressed fault and the area on the fault affected by a pore pressure increase is smaller resulting in a smaller induced seismicity event being triggered. 88 Figure 3.10: (a) Pore pressure distribution on fault plane and (b) corresponding shear slip response to a fluid injection rate of 6 m3/min for 30 minutes for the model scenario without natural fractures (i.e., incipient fractures only). (c) Pore pressure distribution on fault plane and (d) corresponding shear slip response to the same fluid injection rate and duration for the model scenario with natural fractures. These numerical simulations were repeated for a range of different injection rates. It is noted that the injection rates used in the numerical simulations cannot be directly related to those in the empirical database, at least with respect to accounting for differences in fluid viscosity. The numerical simulations assume a viscosity representative of water. Mechanistically, however, they represent comparative responses to low and high injection rates. The results are shown in Figure 3.11. These first suggest that there is an apparent upper bound for the maximum induced 89 seismicity event that can be generated for the assumed model geometry, fault attributes, and input parameters. For the case where a natural fracture network is present, each fluid injection rate produced a larger invaded zone (relative to the models without a natural fracture network) resulting in a very similar large magnitude induced seismicity event. These show a very slight upward trend of increasing event magnitude with increasing injection rate. In contrast, without the presence of a natural fracture network, higher injection rates are required to create the necessary fractures and fluid flow pathways to reach and interact with the fault. Smaller injection rates result in a more limited extent of the hydraulic fracture generated, and therefore a smaller area on the fault being affected by pore pressure increases. Accordingly, the induced seismicity event magnitude is seen to sharply increase with fluid injection rate. These results should be conditioned by noting that the distance between the injection point and the fault (in this case 90 m) is a highly sensitive parameter, while also recognizing that the location of adversely oriented faults relative to a hydraulic fracturing injection point is difficult to establish. 90 Figure 3.11: Magnitude distribution of events calculated from numerical simulations of fault slip for the different fluid injection rates modeled. 3.5 Discussion The results of the empirical and numerical analyses point to the importance of reservoir geology, including the intensity of the natural fracture network relative to the presence of a critically stressed fault, as well as other factors such as the magnitudes and orientation of the principal stresses, when investigating the effect of fluid injection rate and volume on triggering induced seismicity. The implications of these analyses are not limited to hydraulic fracturing activities and can be applied to other practices such as wastewater disposal and enhanced geothermal systems. Specific to the Montney in northeastern British Columbia, the numerical modeling results presented help to provide a mechanistic explanation to statistical relationships observed in the 91 empirical analysis. The Middle Montney is deeper than the Upper Montney so there is already an expectation that the Middle Montney would have a higher susceptibility to higher magnitude events. In addition, though, the empirical data indicates that the Upper Montney is more sensitive to injection rate than the Middle Montney. The numerical models explain this with respect to the relative size of the invaded zone around a hydraulic fracture, which was shown to be much more limited for reservoirs with a lower intensity of natural fractures. A smaller invaded zone results in a lower likelihood of interacting with a critically stressed fault and thus a lower susceptibility to induced seismicity. The invaded zone in this case refers to the extent of pore pressure increase, which is facilitated by the zone of enhanced fracture conductivity generated during hydraulic fracturing. Cipolla et al. (2010) describes this as a “network-fracture conductivity” that incorporates the primary hydraulic fracture as well as a complex network of propped and unpropped secondary fractures constituting the stimulated reservoir volume. Dusseault et al. (2011) further describe this with respect to the presence of natural fractures as a dilated zone that surrounds a multi-stage hydraulic fracture, where natural fractures have been opened permanently by wedging and through asperity propping in response to shear displacements and offset. In this context, a low intensity of natural fractures creates a setting more sensitive to injection rate, with higher injection rates generating more secondary tensile fractures branching off of the primary hydraulic fracture and thus more leak-off and a larger invaded zone. Even so, the extent of this invaded zone is still smaller than that where a more interconnected natural fracture network is present. 92 The presence of an interconnected natural fracture network would result in a larger invaded zone via leak-off into the fractures, increasing the likelihood of interacting with a critically stressed fault. In this case, a higher susceptibility to induced seismicity can be expected together with larger magnitude events, both because of the greater depths involved and therefore higher stresses but also the potential for a larger zone of interaction with a critically stressed fault. Although the depth difference between the Upper and Middle Montney would result in an expectation of higher magnitude events for the latter, the numerical modeling results demonstrated that even when the stresses are kept equal between models without and with the natural fracture network, the latter showed a higher susceptibility to larger magnitude events due to the larger invaded zone interacting with a critically stressed fault 90 m away from the injection point. Thus, another area of interest that builds on the susceptibility of triggering an induced seismicity event is the influence of injection volume and rate on the magnitude distribution of these events. Figure 3.12 and Figure 3.13 respectively show the relation between injection volume and rate with the number of events associated with each seismogenic well, as identified in the empirical analysis. As a first pass, this shows that many wells experience multiple induced seismicity events, and a few are associated with a high number of events (>50). It should be noted that most hydraulic fracturing operations related to unconventional gas development, especially those in the Montney, involve multiple horizontal wellbores drilled from a single well pad, and each horizontal wellbore involves multiple hydraulic fracturing stages along the well axis (i.e., multi-stage hydraulic fracturing). Thus, the potential for multiple induced seismicity events is not unexpected where susceptible conditions have been encountered. For each of these, the marker is 93 color coded to the maximum event magnitude associated with that well. Further divided in each figure are separate plots for the smaller events less likely to be felt (M3). With respect to injection volume, the plot in Figure 3.12 for small events (M3) shows that these are primarily (70%) associated with wells with injection volumes above the average injection volume of 1.4e4 m3. This is consistent with conclusions reached in previous studies in Canada (Farahbod et al., 2015; Schultz et al., 2018; Schultz et al., 2017) and the central and eastern United States (Van der Baan & Calixto, 2017). Furthermore, the data in Figure 3.12 shows that the larger magnitude events (>M3) tend to correspond with wells with a smaller number of associated earthquakes. We interpret this as potentially pointing to differences in the fault attributes and rock strength where stronger faults allow strain energy to build up, thus resulting in fewer but larger events associated with the well. Many of these also have high injection rates plotting above the average for the region (Figure 3.13; see plot for >M3). Again, a higher injection rate potentially points to stronger rock and the need to overcome higher breakdown pressures. For wells associated with higher numbers of earthquakes (>50), the largest of the multiple events are relatively low magnitude events (i.e., M3) induced seismicity events. Similarly, where lower injection rates and volumes result in a high number of lower magnitude events, this might point to a weaker fault system less likely to generate larger magnitude events. Caution is of course still required as these results only speak to general trends in the Montney and the specific conditions for a given well may include a combination of factors that respond in a contrary manner. These results are also relative to a lower magnitude of completeness in the database of ML 1.8. Figure 3.12: Injection volume analysis of wells associated with induced seismicity in the Montney for events with M<3, M>3 and considering all events. The dashed line marks the average injection volume of all wells in the database. 95 Figure 3.13: Injection rate analysis of wells associated with induced seismicity in the Montney for events with M<3, M>3 and considering all events. The dashed line marks the average injection volume of all wells in the database. 3.6 Conclusions Determining the dominant factor in triggering injection induced seismicity has proven to be difficult. This is because susceptibility is often dependent on multiple factors; certain geological and operational conditions must be met for a fluid injection operation to trigger seismicity. However, operational factors such as fluid injection volume, rate and pressure are of special interest as these can be engineered and controlled. An empirical study supported by 3-D numerical modeling was performed to investigate the influence of these operational factors on the magnitude and magnitude distribution of induced seismicity events related to hydraulic fracturing activities in the Montney play of northeastern British Columbia. The results of the 96 empirical analysis found that both injection volume and rate are statistically associated with induced seismicity in this region, but injection pressure is not. Breaking the analysis down further with respect to differences between the Upper and Middle Montney, injection rate was only seen to only be associated with induced seismicity for wells that targeted the Upper Montney. In comparison, injection volume was associated with induced seismicity for wells that targeted Middle Montney. Given the known presence of an interconnected natural fracture network in the Middle Montney, whereas this is unknown for the Upper Montney, we suggest that these differences can be explained mechanistically by the inferred differences in the intensity or characteristics of the natural fracture networks present. This hypothesis was supported by the results derived from a series of 3-D numerical simulations using the distinct-element method and commercial software 3DEC. The results showed that under similar stress conditions and the same fault attributes, properties, and distance from the injection source, the inclusion of a natural fracture network had a significant influence on the susceptibility to fault slip and magnitude of the induced seismicity events. This was related to the size of the invaded zone that developed with the hydraulic fracture generated for a given fluid injection rate. With the presence of an interconnected natural fracture network, the larger invaded zone served to increase both the likelihood of interacting with a critically stressed fault and the area on the fault experiencing a higher pore pressure, resulting in a larger area of slip on the fault and a larger magnitude event. For the scenario without an interconnected fracture network, where the invaded zone was not aided by natural fractures and was restricted to the primary hydraulic fracture and secondary pressure-induced fractures generated, it was seen that the susceptibility to induced seismicity was sensitive to 97 injection rate. In this case, higher injection rates resulted in more secondary pressure-induced fractures and thus a larger invaded zone, which in turn resulted in larger magnitude events. This behavior was not observed in the models that included the fracture network as the natural fractures, more so than the pressure-induced secondary fractures, influenced the extent of the invaded zone. Thus, the modeling results demonstrate that the presence of natural fractures play an important role in influencing the susceptibility of a well to induced seismicity, together with the magnitudes of the events, by increasing the pathways for injection fluid to reach the fault in which a larger area on the fault experiences a higher pore pressure increase. Together, the empirical and numerical results point to the potential of applying injection rate and volume thresholds for mitigating induced seismicity hazards related to hydraulic fracturing operations. Caution is of course required in applying these results to future hydraulic fracturing operations, as the results only consider injection rate and volume, and speak to general trends in the Montney and it must be recognized that conditions for a given well may include a combination of factors that respond in a different or unique manner. 98 Chapter 4: Machine learning analysis of factors influencing induced seismicity susceptibility in the Montney Formation of northeastern British Columbia 4.1 Introduction Unconventional gas resources represent an emerging low-cost, clean burning energy source, the export of which is touted by proponents as presenting a greener transition option to replace more carbon intensive fossil fuels like coal. In northeast British Columbia (NEBC) alone, new discoveries and advancements in extraction technologies have led to estimates of gas-in-place of 1500 Tcf, enough to support development and export operations for more than 150 years. However, with the development of these new resources come new challenges. Amongst these are public, First Nations and regulator concerns regarding induced seismicity associated with hydraulic fracturing and wastewater injection operations. Both activities involve the injection of large volumes of fluids into deep formations, which serve to create localized increases in pore pressures and stress changes acting on critically stressed faults, resulting in fault slip and induced seismicity. Notable events in NEBC include one M4.4 and two M4.6 events between 2014 and 2018. In response to these events, and other environmental concerns, the British Columbia government appointed a scientific panel to review hydraulic fracturing practices and their impacts (Allen et al. 2019). In their review, a key knowledge gap identified (in relation to induced seismicity susceptibility) was that the effects of different geological and operational factors on the spatial 99 and temporal distribution of events are not well understood and vary in importance for different unconventional gas plays. Although separating these from one another is a complex task, it is also recognized that a massive amount of geological, operational, and seismic data are being collected from hydraulic fracturing activities for which robust analysis methods are needed. The rapid development of multivariate statistics and machine learning techniques to analyze large data sets is especially attractive and conducive to this problem, although there is not a lot of experience yet in applying these to induced seismicity hazard assessments (likelihood, severity, etc.), especially in analyzing both operational and geological parameters together. Distinguishing between these is of interest as the influence of geological factors on induced seismicity susceptibility of a targeted formation cannot be controlled or manipulated (outside of avoidance), whereas operational factors (i.e., completion related) can be controlled offering a means to potentially mitigate induced seismicity hazards for a susceptible formation. This chapter presents the results of applying different machine learning algorithms to determine the relative importance of several geological and operational parameters in relation to the triggering of induced seismicity (i.e., feature importance). This is done for data compiled for the Montney Formation in NEBC. The algorithms applied and compared include Decision Tree, Random Forest and Gradient Boost methods. In addition to testing the robustness of these algorithms through a comparative analysis, guidance is provided in the use of machine learning to identify influencing factors as a step towards developing induced seismicity susceptibility maps. 100 4.2 Background A significant increase in the seismicity rate in western Canada in recent years has been associated with the development of unconventional oil and gas resources, including the related activities of hydraulic fracturing (Bao & Eaton, 2016) and wastewater disposal (Schultz et al., 2014). There are numerous operators conducting these activities, each using different operational parameters (e.g., fluid injection volumes and rates) tailored for the local geological setting and targeted formation, as well as further shaped by in-house objectives, experiences and optimization efforts. This raises the question of what are the effects of these parameters on induced seismicity susceptibility and magnitude distribution? The question of susceptibility addresses the likelihood that a particular well will generate induced seismicity; this can be viewed as a classification problem (seismogenic or not seismogenic). The question of magnitude distribution addresses the potential severity. It is unlikely that a single causative factor is solely responsible for an induced seismicity event. Instead, multiple parameters, both geological and operational, can play an influencing role. With respect to operational parameters, it has been argued that the moment release attributable to induced earthquakes is related to the net volume of the injected fluid, with empirical trends established that link an upper limit for the moment magnitude to injection volume (Hallo et al., 2014; McGarr, 2014). The data analyzed in these studies included a mix of hydraulic fracturing and waste water disposal activities in sedimentary rocks and enhanced geothermal projects in crystalline rocks, both involving different regions (Europe, United States, Australia). Weingarten et al. (2015) carried out a similar study combining information from public databases on injection wells with available earthquake catalogues but instead concluded that injection rate is 101 the most important operational parameter affecting induced seismicity. Their study focussed on data from hydraulic fracturing and wastewater disposal activities in the eastern and central United States. For the Western Canadian Sedimentary Basin (WCSB), where the Montney formation which is the focus of this study is situated, different studies generally conclude that injection volume is the most important operational factor in triggering induced seismicity. Schultz et al. (2018) investigated the relationship between injection parameters and induced seismicity in the Duvernay shale play in Alberta and concluded that events are associated with completions that used larger injection volumes and that seismic productivity scales linearly with injection volume. Their analysis further showed that injection pressure and rate have an insignificant association with seismic response, and that geological factors account for the variability in induced seismicity susceptibility observed in the region. The results of previous chapter show that both injection volume and rate influence triggering induced seismicity in this region. With respect to geological parameters, Gobel (2015) compared several fluid injection operations in California and Oklahoma and examined the temporal and spatial variations in their induced seismicity responses. His results suggest that operational parameters for fluid injection are likely of secondary importance and that the primary controls on injection induced seismicity are the site-specific geology and geological setting. Van der Baan and Calixto (2017) compared current and historic seismicity rates in six U.S. states and three Canadian provinces to past and present oil and gas production. Their study showed that in addition to injection volumes, local and regional-scale geology and tectonics influence increased earthquake hazard susceptibility. Amini and Eberhardt (2019) similarly compared induced seismicity and well data for several key North 102 American unconventional gas plays, with focus on magnitude distribution relative to differences in the tectonic in situ stress regime. They found that stress regime has a significant influence on event magnitude with thrust fault stress regimes, as exists in parts of the Montney, being more susceptible to large magnitude events compared to a strike-slip or normal fault stress regimes. In Oklahoma, Shah and Keller (2017) combined geophysical and drill hole data to map subsurface geologic features in the crystalline basement and found that most induced seismicity events are located where the crystalline basement is likely composed of fractured intrusive or metamorphic rock; areas with extrusive rock or thick sedimentary cover (>4 km) exhibited little induced seismicity. They concluded that the differences in seismicity may be due to variations in permeability structure – within intrusive rocks, fluids can become narrowly focused in fractures and faults, causing a concentrated increase in local pore fluid pressures, whereas more distributed pore space in sedimentary and extrusive rocks may relax pore fluid pressures. Hincks et al. (2018) developed an advanced Bayesian Network to model joint conditional dependencies between spatial, operational, and seismicity parameters in Oklahoma. They found that injection depth relative to crystalline basement most strongly correlates with seismic moment release and that the combined effects of depth and volume are critical, as injection rate becomes more influential near the basement interface. Similar findings were reported by Skoumal et al. (2015) and Currie et al. (2018) for hydraulic fracturing operations in Ohio. The latter showed that seismicity occurred along faults below the injection interval in the Precambrian crystalline basement. From seismic reflection lines, they showed that these fault systems intersected the injection interval targeted by the well, providing permeability pathways for fluid pressure increases leading to fault slip. 103 Specific to the geology of the WCSB, Schultz et al. (2016) found that hypocenters of induced seismicity clusters in Alberta coincided with the margins of the Devonian carbonate reefs and interpreted this spatial correspondence as the result of geographically biased activation potential, possibly as a consequence of reef nucleation preference to paleobathymetric highs associated with Precambrian basement tectonics. Their work provided evidence that in some areas Paleozoic and Precambrian strata are likely in hydraulic communication which points to the important role of regional- and local-scale geological factors in the nature of induced seismicity. Eaton and Schultz (2018) also suggested natural processes involving the transformation of organic material (kerogen) into hydrocarbons and cracking to produce gas, can cause fluid overpressures resulting in an increased susceptibility to induced seismicity. They presented two examples from the WCSB where induced seismicity attributed to hydraulic fracturing are strongly clustered within areas characterized by high pore-pressure gradients. The above examples highlight the importance of different operational and geological factors on induced seismicity, but do so from the perspective of studying the influence of a single factor. It is unlikely, however, that a single causative factor is solely responsible for an induced seismicity event. Instead, multiple factors, both geological and operational, can play an influencing role. Therefore, it is important to consider and understand the cause and effect of different geological and operational factors on the spatial and temporal distribution of induced seismicity events. However, this is not a simple task and requires probing a wide variety of linear and non-linear associations and interaction terms between factors affecting induced seismicity without assuming a priori knowledge on the nature of the relationships between these factors. 104 The use of machine learning and data analytics are quickly evolving as a means to identifying hidden patterns and extracting information from large data sets. In the geosciences and rock engineering, it has been applied to predicting rockburst potential in deep mines (Pu et al., 2018; Ribeiro e Sousa et al., 2017) and squeezing behavior in deep tunnels (Sun et al., 2018), as well as developing geological maps using remote sensing data (Cracknell & Reading, 2014) and analyzing laboratory data from rock testing (Millar & Clarici, 2002) and blasting (Liu & Liu, 2017). In the context of earthquake seismology, machine learning has been applied to a variety of problems such as earthquake identification (Rouet-Leduc et al., 2017) and forecasting (Panakkat & Adeli, 2009). Building on these studies, it is recognized that a massive amount of geological, operational, and seismic data is being collected with hydraulic fracturing activities, and the size and complexity of these datasets have made traditional empirical and statistical analyses inefficient and ineffective. This has led to recent studies by Pawley et al. (2018) who combined tectonic, geomechanical and hydrological data with induced seismicity data related to hydraulic fracturing operations in the Duvernay play in Alberta to train a logistic regression algorithm to map and develop an induced seismicity potential map. Their results suggest that the proximity to basement, formation overpressure, minimum horizontal stress, proximity to reef margins, lithium concentrations, and natural seismicity rate are the dominant controlling factors to triggering induced seismicity within the study area. Zhang et al. (2020) used machine learning on real-time induced seismicity location problems for small events in Oklahoma by accessing seismic waveform data from a regional network. They designed a fully convolutional network (FCN), to predict a 3-D image of the earthquake location probability from a volume of input data recorded at multiple network stations. Their results showed that the designed system is capable of locating small events of ML ≥ 2.0 with a mean epicenter error of 4-6 km. 105 4.3 Data Compilation and Preparation A database of 16,945 hydraulic fracturing stages within the Montney Formation was compiled and analyzed using multiple sources reporting well activities in British Columbia (BC Oil & Gas Commission, 2018; geoLOGIC systems Ltd., 2019). This was combined with a second database involving a comprehensive earthquake catalogue compiled for NEBC and western Alberta (Visser et al., 2017). This catalogue was produced specifically to study induced seismicity in this region and consists of 4916 events for the period of January 2014 to December 2016 with a magnitude of completeness of 1.8 (ML). To prepare the data for analysis using a supervised machine-learning algorithm (as discussed in the next section), it was necessary to determine the output labels. Here, induced seismicity was considered as a binary-classification problem with respect to the observed seismic activity. Wells were classified as being either ‘aseismic’ or ‘seismic’ based on spatial and temporal correlations with hydraulic-fracturing operations. This was done by cross-correlating the earthquake catalogue with the well database and applying a series of spatial and temporal filters to identify the subset of earthquake events that are likely induced seismicity events related to hydraulic fracturing and waste water disposal activities. The first step was to clip the data to only include earthquakes located within the boundaries of the Montney study area and to filter out events spatially associated with non-oil and gas related anthropogenic activities, such as those from mining and construction (i.e., blasting). This step reduced the total number of events being considered from 4916 to 2867. Next, a spatial filter was applied to search for all event locations that were within a 5 km radius of an active hydraulic fracturing well. The 5 km radius represents 106 the uncertainty in the event location accuracy reported for the earthquake catalogue. To this, a 3-month temporal filter was applied (e.g., Atkinson et al., 2016). Thus, if an earthquake event was recorded as occurring within 5 km of an active well and within 3 months from the start date of the hydraulic fracturing activity, it was considered here to be an induced seismicity event and the corresponding well was classified as being seismogenic. This resulted in a subset of 543 events identified as induced seismicity events. The input parameters (from now on will be referred to as features as it’s a more common name for parameters when discussing machine learning) for the machine learning algorithms were selected from available geological and operational data in the Montney. The values of these features were either directly available at each well or were interpolated using the average of the three closest data points. The operational features were treated differently as only wells that had complete data throughout all features were included. After data processing, data from 11,415 stages (out of 16,945) were used for the machine learning analysis. The input features are described in Table 4.1. Due to data sparsity and the limited number of measurements, the pressure-gradient value on each well was estimated using the average of the three closest data points. The resulting pressure gradients ranged from 5.4 to approximately 18 kPa/m. In addition to the reservoir pore pressure, information regarding the maximum horizontal-stress (𝑆𝐻𝑚𝑎𝑥) direction was included. This was calculated for each well based on 𝑆𝐻𝑚𝑎𝑥 azimuths extracted from the World Stress Map database (Heidbach et al., 2018) and interpolated for each well. Two different values were investigated: 1) the difference between the local 𝑆𝐻𝑚𝑎𝑥 and the horizontal well azimuths; and 2) the difference 107 between the local and regional 𝑆𝐻𝑚𝑎𝑥 azimuths, with N 45° E considered as the regional 𝑆𝐻𝑚𝑎𝑥 direction in the Montney. Injection depth was also considered as a proxy for the magnitude of stresses in this region since no direct measurement was available. Specific to the local geology, the vertical distance between the injection depth and the top of the Montney and Debolt formations were considered, together with that to the Precambrian basement. For these, a negative value indicates above the formation top and a positive value refers to below the formation top. It should be noted that there is a high degree of uncertainty in the interpolated values for the top of the basement due to a lack of measurements. The final geological data included in the machine learning analysis was the distance from the injection point to mapped faults. For this parameter, the shortest distance between the well head and closest fault was considered. 108 Table 4.1: Input features used for machine learning analysis. Feature name Abbreviation Source Number of Available Data Points Feature Category Lower limit Upper limit Median Pore Pressure Gradient (KPa/m) PP_grad BCOGC1 2252 Geological 5.4 18 12.4 Local and Regional SHmax Azimuth Difference (deg) Az_L_R World Stress Map2 601 Geological 0.3 12 5.6 Local SHmax and Well Direction Azimuth Difference (deg) Az_L_W BCOGC1 16945 Operational 0.04 88 41.7 Distance to Faults (m) Dist_F Geoscience BC3 16945 Geological 2.2 70396 1772.8 Average Injection Rate (m3/min) Rate geoLogic4 15472 Operational 0.16 22 8.8 Injection Depth (m) Inj. Depth geoLogic4 16945 Operational 1312 3416 2052 Stage Injection Volume (m3) Volume geoLogic4 16603 Operational 0.99 5645.22 592 Maximum Injection Pressure (MPa) Max_P geoLogic4 13889 Operational 4.3 99.48 56.1 Distance to top of Montney (m) D_Mont BCOGC1, geoLogic4 8232 Geological -100 427.6 78.6 Distance to top of Debolt (m) D_Deb BCOGC1, geoLogic4 2183 Geological -478 8.6 -247 Distance to Basement (m) D_Base BCOGC1, geoLogic4 28 Geological -2449.2 -129.9 -1661.7 Well Completion Length (m) Comp_Len geoLogic4 16945 Operational 91.59 3835.44 1614.1 1 BC Oil & Gas Commission (2018), 2 Heidbach et al. (2018), 3 GeoscienceBC (2018), 4 geoLOGIC systems Ltd. (2019)109 4.4 Machine Learning Algorithm Development Machine learning can be undertaken using supervised or unsupervised algorithms. Supervised learning is where the input features and an output result are given, and an algorithm is used to learn the mapping function between these. The goal is to approximate the mapping function so well that for any new input data, the output can be predicted for that specific data. This contrasts with unsupervised learning where only the input data is known, and no corresponding output variables are given. The goal for unsupervised learning is to model the underlying structure or distribution in the data in order to learn more about the data. For this study, supervised learning was used as the initial data analysis (as described in the previous section) included identifying which wells were associated with induced seismicity and which were not. These represent the correct answers to the classification problem for training the mapping function; the corresponding data associated with each set of wells is referred to as the training data. Three different supervised machine learning algorithms were used that are generally considered to be robust for classification problems: Decision Tree, Random Forest, and Gradient Boost. These methods were chosen because of the ease of interpretability of their results and also because they are not sensitive to the scale of input data. The objective of the algorithm is to iteratively make predictions on the training data and to correct these until the algorithm achieves an acceptable level of performance. Decision trees are a non-parametric supervised learning method used for classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. Decision trees have two advantages: the resulting model can easily be visualized and understood by non-experts, and the algorithms are completely invariant to scaling of the data. As each 110 feature is processed separately, and the possible splits of the data do not depend on scaling, no preprocessing of features is needed for decision tree algorithms. The main limitation of decision trees is that they tend to over fit the data and provide poor generalization performance. A random forest is essentially a collection of decision trees, where each tree is slightly different from the others. The idea behind random forests is that each tree might do a relatively good job of predicting but will likely over fit part of the data. If many trees are built, all of which work well and over fit the data in different ways, the amount of overfitting can be reduced by averaging their results. The gradient boost regression tree is another ensemble method that combines multiple decision trees to create a more powerful model. These can be used for both regression and classification. In contrast to the random forest approach, gradient boosting works by building trees in a serial manner, where each tree tries to correct the mistakes of the previous one. The main idea is to combine many simple models (known as weak learners) that can provide good predictions on parts of the data, and so more and more trees are added to iteratively improve performance. All three machine learning models were built using Scikit-learn, a Python library for machine learning (Pedregosa et al., 2011). The dataset was divided into training and validation sets, accounting for 75% and 25% of the full dataset, respectively. The training dataset was further divided into training and test datasets which were used to train the algorithm using 50 cross validations, with the training set accounting for 98% of the training dataset and the test set for 2% at each cross validation run (Figure 4.1). 111 Figure 4.1: Splitting of the compiled dataset into training and testing sub-sets, with the testing data to be used for validation. The training data was used to train and evaluate the optimum tree depths of the three algorithms using 50-fold cross validations. Figure 4.2 shows the results of 50 cross validations for each algorithm. At each run, the accuracy of the model for a specific tree depth is calculated. In these figures the orange line represents the accuracy of the training set. The blue line shows the mean cross validation accuracy and the shaded area represents the confidence interval (±2 standard deviations) for the calculated means. For these plots, an accuracy of 1 represents 100% accuracy. This determines if the training dataset is over fitted, and alongside this, the confidence interval was used to determine the optimum tree depth. Based on the cross-validation results, the depths of 16, 12 and 8 were chosen for the Decision Tree, Random Forest and Gradient Boost algorithms, respectively. 112 Figure 4.2: Results of the 50-fold cross-validations to determine the optimum tree depth for each algorithm, showing from top to bottom: Decision Tree, Random Forest, and Gradient Boost. 113 4.5 Results 4.5.1 Feature Importance The results of the cross-validations using each algorithm were further analyzed to investigate the importance of each feature on the classification outcome. Figure 4.3 shows the feature importance plots identified by each of the three classification algorithms. The bars are color coded to differentiate features that relate to the geology from those that are operational. The importance of each feature is indicated as a coefficient; these coefficients do not have a physical meaning and are compared based on the relative values and not the absolute values. Table 4.2 similarly compares the relative ranking of the features considered as determined by each algorithm. Based on the result of the decision-tree algorithm, the high importance features were determined to be pore pressure gradient, well completion length, azimuth difference between the local and regional maximum horizontal stress, distance to basement and distance to faults. Four out of these five features also form the top five ranked features from the random forest analysis, although in a slightly different order and with azimuth difference between the local and regional maximum horizontal stress being replaced by injection depth which is ranked as high importance. For the gradient boost analysis, again four of these features were ranked in the top five, the exception being that this model showed a higher sensitivity to the well direction than the distance to faults. The gradient boost model also showed very high sensitivity to azimuth difference between the local and regional maximum horizontal stress and azimuth difference between local maximum horizontal stress and well direction. Overall, in all models, pore 114 pressure gradient, completion length and distance to basement were among the high importance features. In all models no change in the grouping of operational features were observed; injection rate and maximum injection pressure were ranked lowest in importance, and injection volume ranked in the middle. This is an interesting result because injection rate and volume are often cited as operational features that have a significant influence on induced seismicity (McGarr, 2014; R. Schultz et al., 2018). This appears to hold partly true in the case of injection volume, but operational features such as well completion length and injection depth, which have not been thoroughly studied, appear to have a stronger correlation with induced seismicity. 115 Figure 4.3: -Feature importance calculated using three different supervised machine learning algorithms: Decision Tree, Random Forest, and Gradient Boost. 116 Table 4.2: Relative ranking of feature importance based on the different machine learning results. Rank Decision Tree Random Forest Gradient Boost 1 Pore Pressure Gradient Well Completion Length SHmax Azimuth 2 Well Completion Length Distance to Basement Well Direction Azimuth 3 SHmax Azimuth Distance to Faults Distance to Basement 4 Distance to Basement Injection Depth Pore Pressure Gradient 5 Distance to Faults Pore Pressure Gradient Well Completion Length 6 Well Direction Azimuth Well Direction Azimuth Injection Volume 7 Injection Depth SHmax Azimuth Distance to Faults 8 Injection Volume Injection Volume Injection Depth 9 Distance to top of Debolt Distance to top of Debolt Distance to top of Debolt 10 Distance to top of Montney Distance to top of Montney Injection Rate 11 Injection Pressure Injection Rate Distance to top of Montney 12 Injection Rate Injection Pressure Injection Pressure 4.5.2 Model Validation To validate the trained algorithms, they were next applied to the test data that were set aside to evaluate each model’s performance (see Figure 4.1). The test data comprised 25% of the full dataset and was not previously used to train the algorithms. To evaluate the performance of each algorithm, a confusion matrix was calculated. Also known as an error matrix, the confusion matrix allows visualization of the performance of a supervised machine learning algorithm by reporting the number of true and false positives and true and false negatives (Figure 4.4). These are based on predictions using the test data relative to the mapping functions determined from 117 the training data. In this case, a true positive would be a correct prediction that a well is associated with induced seismicity and a true negative would be a correct prediction that the well is not. Similarly a false positive would be the prediction of a well being associated with induced seismicity where there was none, and a false negative would be a prediction of a well not having experienced induced seismicity when it had. The results from calculating a confusion matrix for each algorithm are shown in Figure 4.5. Comparing these, the random forest and gradient boost methods performed slightly better than the decision tree classifier. However, all three models performed with very high accuracy (97-98%). Figure 4.4: Calculation and use of a confusion matrix to compare test data against correlations derived from machine learning using training data. 118 Figure 4.5: Confusion matrixes comparing predicted versus actual results for the trained decision tree, random forest and gradient boost models against the validation test data. IS: induced seismicity events; No-IS: no induced seismicity events. To further interpret the results, a SHapley Additive exPlanations (SHAP) analysis was run. SHAP is a game theory approach used to help interpret predictions from complex models, for example the output from machine learning models. SHAP assigns each feature an importance value for a particular prediction and shows there is a unique solution for each class of additive feature importance that adheres to desirable properties (Lundberg & Lee, 2017). The SHAP Tree Explainer tool is a subcategory of SHAP that is specifically built for interpreting tree models such as decision trees and random forests. The SHAP value plot can show the positive and negative relationships of the predictors with the target variable. The analysis presented here is for the random forest results as these performed slightly better than the decision tree model. Figure 4.6 presents the summary plot from the SHAP analysis, which combines feature importance with feature impact. Based on this plot, the following information can be gained. First, each feature is ordered according to its importance (starting with the most important). Note that the SHAP plot is calculated for one instance of the random forest model, whereas the 119 ranking of feature importance in Figure 4.3 and Table 4.2 is based on an averaging of 50 cross validation runs. Thus, the order of feature importance between the two might be slightly different. Next, points are plotted to show the distribution of the SHAP values using color to represent the feature value and stacking of overlapping points in the y-axis direction to give a sense of the distribution of the SHAP values. From this the impact (both positive and negative) is shown through the horizontal location of stacking, which shows whether the effect of that value is associated with a higher or lower prediction. This can be compared to whether the value for that variable/observation is high (red) or low (blue). Thus, for example, it can be seen that high values of completion length have a high and positive impact on the quality rating. The high values related to this feature are indicated by the red color of the points, and the high, positive impact is shown by its extent on the X-axis. The results of the random forest SHAP feature analysis helps to validate the meaningfulness of the algorithm results. Inspecting both Figure 4.3 and Figure 4.6, it can be seen that key influencing features such as completion length, pore pressure gradient and injection depth have a positive correlation with induced seismicity. The influence of pore pressures in the Montney has been studied by Eaton and Schultz (2018), who demonstrated a positive relationship between over-pressured areas and induced seismicity. The positive correlation of completion length is also valid as higher completion lengths correspond with larger stimulated volumes and therefore a higher probability of adversely interacting with a critically stressed fault. Features such as distance from the basement or distance to a known fault have negative correlations, meaning shorter distances between the injection point and basement or fault increase the likelihood of triggering an induced seismicity event. 120 Figure 4.6: Results for a SHAP feature importance analysis using the Random Forest trained model. 4.6 Discussion Machine-learning models are highly dependent on the quality and quantity of the input data. For this analysis, where data for a feature was either limited or the spatial distribution and/or coverage of the data was sparse relative to the distribution of the wells, this was compensated for using linear interpolation. However, large distances between points can reduce the accuracy of interpolation, as can the interpolation method itself (e.g., assigning linear vs non-linear weightings). With time, as new data becomes available, including that for other features, the induced seismicity susceptibility model can be updated to improve its predictive capabilities. 121 Based on the results obtained, an interesting observation is the correlation of injection volume and SHAP values. As can be seen in Figure 4.6, high injection volumes have a negative correlation with triggering induced seismicity. This might be interpreted as high injection volumes reduce the risk of induced seismicity, which is counter to general experience. Thus, empirical analyses and machine learning data correlations for feature analysis do have their limitations and should be constrained by an understanding of the physics of fault slip and induced seismicity mechanisms. When comparing the ranking of feature importance, the decision tree and random forest models provided similar results. In these models, high importance features included the geological parameters pore pressure gradient, distance to basement and distance to faults. The operational features well completion length and injection depth were also highly ranked by the random forest model. For the gradient boost method, the ranking was slightly different with the importance of the stress field showing the greatest influence, both with respect to the geological feature 𝑆𝐻𝑚𝑎𝑥 and related operational feature of the well direction relative to 𝑆𝐻𝑚𝑎𝑥 azimuth difference. This difference must be related to differences in how these models work: gradient boost models are based on shallow trees (high bias, low variance) and they reduce error mainly by reducing bias. On the other hand, decision tree and random forest models use fully grown trees (low bias, high variance) and they reduce model’s error by reducing variance. Bias are the simplifying assumptions made by a model to make the target function easier to learn, and variance is the amount that the estimate of the target function will change if different training data was used. The source of bias in this problem is the number of features that are used for classification and by including both operational and geological features the overall bias tends to be less than just 122 considering a single category of features. On the other hand, the variance of the data can be calculated and as seen in Table 4.3 it is high for parameters such as volume and distance to faults whereas for stress features it is low. This difference in variance of the data can help to explain the observed difference in the results of these algorithms. The values of variances reported in Table 4.3 are divided by mean for each feature in order to make them unitless and therefore comparable. Table 4.3: Variance of some of the features used in this study Feature Az_L_R* Az_L_W† Dist_F‡ Volume⁋ Variance/Mean 17.8 2.2 19113.1 370.4 *Local and Regional SHmax Azimuth Difference †Local SHmax and Well Direction Azimuth Difference ‡ Distance to Faults ⁋ Injection Volume The comparison of feature importance between geological and operational parameters show that overall geological parameters generally ranked higher in importance. In all models, the operational parameters ‘average injection rate’ and ‘average injection pressure’ consistently were ranked as the least influential. The only high importance operational features were the completion length and the horizontal well direction relative to the 𝑆𝐻𝑚𝑎𝑥 azimuth. As was shown in Figure 4.6, completion length had a positive correlation with seismicity and can be thought of as the impacted volume being stimulated by hydraulic fracturing. The larger the stimulated volume, the higher the probability of intersecting a fault (directly or indirectly). 123 4.7 Conclusions The application of machine learning was investigated for the purpose of ranking the influence of geological and operational parameters on the classification problem of induced seismicity susceptibility (i.e., distinguishing between wells that are associated with induced seismicity and those that are not). Three different algorithms: Decision Tree, Random Forest and Gradient Boost, were tested using data related to hydraulic fracturing activities in the Montney Formation in northeastern British Columbia. All models were initially trained on a subset of 75% of the total data compiled using a 50-fold cross validation analysis. The remaining 25% of the data formed a validation set to test the trained models. The validation results showed a high accuracy of successful predictions (97-98%) for all three models. The classification results were used to calculate the relative importance of all features on whether a well had or had not been associated with induced seismicity. Geological features were differentiated from operational features as the latter are of particular interest as they can be controlled or manipulated to mitigate induced seismicity hazards. However, it was the geological features that generally rated higher with respect to correlation with wells associated with induced seismicity. In all models, pore pressure gradient (hydrostatic versus overpressured) rated highly as having a major influence. For the decision tree and random forest trained models, distance to basement and distance to known faults also ranked highly, whereas for the gradient boost, the 𝑆𝐻𝑚𝑎𝑥 azimuth was the geological feature that ranked highly. For the operational features, the completion length was ranked as high importance. 124 Overall, the results of these analysis agree with the current understanding of the parameters that influence induced seismicity susceptibility such as reservoir over pressure (Eaton & Schultz, 2018), influence of stress regime (Amini & Eberhardt, 2019), and stimulated injection volume (Schultz et al., 2018). These point to the importance of understanding the geology of the Montney formation including the 3-D seismic mapping of faults and making in situ stress measurements. The machine learning algorithms investigated here can be used to better understand induced seismicity by determining and ranking the factors that influence induced seismicity susceptibility and therefore further improve industry practices and regulator oversight. 125 Chapter 5: Conclusions 5.1 Summary Induced seismicity related to hydraulic fracturing and wastewater disposal is a key subject of public, industry, regulator and government concern. This thesis presents the results from a series of detailed and advanced data analysis and numerical modelling investigations into the different geological and operation factors that influence the likelihood of induced seismicity and the corresponding magnitude distribution. Three research questions were identified, and hypotheses posed, which formed the central objectives of the thesis: • In Chapter 2, the first objective addressed the question of how do tectonic stress regimes influence the magnitude and magnitude distribution of induced seismicity events. For this, the dominant stress regime for several major shale gas basins across North America was determined and the corresponding b-values of induced seismicity events in these basins calculated. In parallel, results from a series of 3-D numerical simulations were used to provide mechanistic insights to the data analysis results. • In Chapter 3, the second research objective investigated the sensitivity of fault slip and induced seismicity to operational factors, with specific focus on data and experiences from northeastern British Columbia. Operational parameters were posed to be of special interest as they offer a means to mitigate induced seismicity hazards through their control or manipulation. An empirical analysis was performed to investigate the association of injection rate and volume with induced seismicity in the Montney Formation, with distinctions being made relative to different depths within the Montney. Again, results 126 from a series of 3-D numerical simulations were used to provide mechanistic understanding into the findings from the empirical analysis. • In Chapter 4, the third research objective applied three machine learning algorithms and compared their results to investigate the relative importance of different geological and operational parameters on well susceptibility to induced seismicity. Focus was again placed on data and experiences from the Montney region of northeastern British Columbia. Algorithms were trained to classify the data and predict ‘Seismic’ versus ‘Aseismic’ wells. Together, the findings of the research presented in this thesis help to improve our understanding of the mechanisms by which induced seismicity are triggered, their sensitivities relative to local site conditions, and the factors that influence and control the corresponding range of event magnitude distributions. These findings can be used to help improve forecasting and mitigation strategies to allow designers and operators to better assess the likelihood and therefore risk of large induced seismicity events. Furthermore, the methods and procedures developed, as well as the key findings obtained, serve to guide future researchers investigating knowledge gaps associated with our understanding of induced seismicity. 5.1.1 Influence of Tectonic Stress regime The influence of different tectonic stress regimes on the magnitude and magnitude distribution of induced seismicity events related to hydraulic fracturing activities were investigated by performing an empirical analysis supported by a series of 3-D numerical modelling simulations. 127 For this, a comprehensive database of existing in situ stress indicators, earthquake catalogues and hydraulic fracturing well data for key North American shale gas plays was compiled. The stress data showed that the stress regime for the Montney play in northeastern British Columbia is predominantly a thrust fault regime (i.e., high horizontal stress relative to the vertical stress with depth). In comparison, the stress regime for the Woodford play in Oklahoma was seen to be that of a strike-slip fault regime (i.e., a horizontal stress similar to the vertical stress with depth). For the seismicity and well data, a spatial and temporal filter was applied to identify those wells associated with induced seismicity events for each play. The resulting frequency-magnitude plots showed a lower b-value (i.e., events skewed towards larger magnitudes) for induced seismicity occurring in a thrust fault regime (Montney) compared to a strike-slip regime (Woodford). These findings were supported by the results obtained for the 3-D numerical simulations using the distinct element method. The numerical modelling results showed that under similar geological and operational conditions, the tectonic in situ stress regime has a significant influence on the amount of strain energy that can be stored on the fault and therefore released through injection induced slip. The results showed that both the magnitude and area of slip are greater in the case of a thrust fault stress regime, compared to the other fault/stress regimes if all other factors are kept equal. This is due to the higher horizontal stresses for the thrust fault case and how the stress field resolves, relative to the fault dip angle, to act on the fault. This combination acts to promote the stored strain energy being released as a single large event. In contrast, numerical simulations for the other faulting types showed these to be more susceptible to either smaller slip magnitudes or 128 smaller slip areas with the stored strain energy being released as multiple smaller slip events. Accordingly, the modelling results show a higher upper bound can be expected in the distribution of induced seismicity events when operating in a thrust fault stress regime compared to a strike-slip or normal fault stress regime. 5.1.2 Influence of Operational Parameters An empirical study supported by 3-D numerical modelling was performed to investigate the influence of several key operational factors (fluid injection volume, rate and pressure) on the magnitude and magnitude distribution of induced seismicity events related to hydraulic fracturing activities in the Montney play of northeastern British Columbia. These parameters are of special interest as these can be engineered and controlled to potentially mitigate induced seismicity hazards. The results showed that both injection volume and rate are statistically associated with induced seismicity in this region, but injection pressure is not. Furthermore, breaking the analysis down to the different depths within the Montney formation being targeted by the hydraulic fracturing operations, showed that injection rate was only associated with induced seismicity for wells in the upper parts of the Montney. In comparison, injection volume was associated with induced seismicity for wells that targeted middle parts of the Montney. Given the known presence of an interconnected natural fracture network in the middle Montney, whereas this is unknown for the upper Montney, a mechanistic model was postulated suggesting that these differences can arise due to differences in the intensity or characteristics of the natural fracture networks present. 129 This hypothesis was supported by the results derived from a series of 3-D numerical simulations using the distinct-element method. Simulation results showed that under similar stress conditions, fault attributes, material properties, and distance from the injection source, the inclusion of a natural fracture network had a significant influence on the susceptibility to fault slip and magnitude of the induced seismicity events. This was related to the size of the invaded zone that developed with the hydraulic fracture generated for a given fluid injection rate. With the presence of an interconnected natural fracture network, the larger invaded zone served to increase both the likelihood of interacting with a critically stressed fault and the area on the fault experiencing a higher pore pressure, resulting in a larger area of slip on the fault and a larger magnitude event. For the scenario without an interconnected fracture network, where the invaded zone was not aided by natural fractures and was restricted to the primary hydraulic fracture and secondary pressure-induced fractures generated, it was seen that the susceptibility to induced seismicity was sensitive to injection rate. In this case, higher injection rates resulted in more secondary pressure-induced fractures and thus a larger invaded zone, which in turn resulted in larger magnitude events. This behavior was not observed in the models that included the fracture network as the natural fractures, more so than the pressure-induced secondary fractures, influenced the extent of the invaded zone. Thus, the modeling results demonstrate that the presence of natural fractures play an important role in influencing the susceptibility of a well to induced seismicity, together with the magnitudes of the events, by increasing the pathways for injection fluid to reach the fault in which a larger area on the fault experiences a higher pore pressure increase. 130 Together, the empirical and numerical results point to the potential of applying injection rate and volume thresholds for mitigating induced seismicity hazards related to hydraulic fracturing operations. Caution is of course required in applying these results to future hydraulic fracturing operations, as the results only speak to general trends in the Montney and it must be recognized that conditions for a given well may include a combination of factors that respond in a different or unique manner. 5.1.3 Susceptibility to Induced Seismicity The application of machine learning was investigated for the purpose of ranking the influence of geological and operational parameters on the classification problem of induced seismicity susceptibility (i.e., distinguishing between wells that are associated with induced seismicity and those that are not). Three different machine learning classification algorithms: Decision Tree, Random Forest and Gradient Boost, were tested using data related to hydraulic fracturing activities in the Montney play in northeastern British Columbia. All three models showed high accuracy (97-98%) in correctly classifying whether wells were associated with seismicity or not. In terms of computing time, the Decision Tree classifier had the fastest run time followed by Random Forest and Gradient Boost Classifiers. Overall, geological factors were seen to generally rate higher than operational factors with respect to correlation with wells associated with induced seismicity. In all models, pore pressure gradient (hydrostatic versus overpressured) rated highly as having a major influence. For the decision tree and random forest trained models, distance to basement and distance to known 131 faults also ranked highly, whereas for the gradient boost, the 𝑆𝐻𝑚𝑎𝑥 azimuth was the geological feature that ranked highly. For the operational factors, completion length and well direction relative to the local 𝑆𝐻𝑚𝑎𝑥 azimuth difference were ranked as high importance features. These findings point to the importance of understanding the geology of the Montney formation including the 3-D seismic mapping of faults and making in situ stress measurements. The machine learning algorithms investigated can be used to better understand induced seismicity by determining and ranking the factors that influence induced seismicity susceptibility and therefore further improve industry practices and regulator oversight. 5.2 Key Scientific Contributions The key scientific contributions made in this thesis that contribute to the body of research and knowledge on induced seismicity are as follows: i) Compilation and analysis of a comprehensive database related to hydraulic fracturing activities in northeastern British Columbia, with focus on the Montney, comprising public-domain geology, in situ stress, and operational data. This database alongside with the developed filtering techniques were used to show that: wells in regions with high horizontal stresses are more susceptible to larger magnitude events; both injection volume and rate are associated with the triggering of induced seismicity; and geological 132 factors generally rank higher than operational parameters with respect to association with induced seismicity susceptibility and magnitude distribution. ii) Demonstration of the utility and importance of using large, reservoir-scale, 3-D numerical modelling, specifically the distinct-element method, to simulate hydraulic fracturing and injection induced seismicity. The procedures developed were shown to be particularly effective in providing mechanistic understanding of data trends and associations observed from empirical analysis. For example, it was shown that the presence of a natural fracture network can significantly influence the susceptibility of a well to fault slip and generating an induced seismicity event. This relates to the natural fractures increasing the connectivity for fluid flow leading to larger invaded zones and therefore increased likelihood of higher pore pressures intersecting and interacting with a critically stressed fault. iii) Demonstration and guidance in the use of machine learning algorithms to uncover associations in large data sets comprised of geological and operational data integrated with induced seismicity data derived from hydraulic fracturing operations. This included the use of different plots of machine learning output to analyze and visualize the relative importance of various geological and operational features. A first-of-its-kind analysis and ranking of feature importance specific to the Montney Formation in northeastern British Columbia was provided. It was shown that geological factors generally rate higher than operational factors with presence of overpressures, distance to basement, distance to known faults and SHmax azimuth being the geological features that ranked highly. For 133 operational factors, completion length and well orientation were also ranked highly as important features. 5.3 Future Works Several limitations and opportunities for future work have been identified based on the results presented in this thesis. These include: • It was demonstrated that the quantity and quality of data is of central importance in any empirical analysis, especially in the training of machine learning algorithms. As experiences in the Montney, northeastern British Columbia and Western Canadian Sedimentary Basin continue, the collection of additional geology, stress state and operational data will create opportunities to further improve the results of the data analysis and machine learning investigations demonstrated in this thesis. Of specific value would be the addition of 3-D seismic data to that made available to researchers, or given the sensitivity of this data to the commercial interests of those who collect the data, the sharing of interpreted faults derived from 3-D seismic data. • The continued collection of data will also allow for statistical analyses to consider longer time intervals. The analyses conducted in this thesis were limited to available data and earthquake catalogues for the period of January 2014 to December 2016 (with a data freeze being required for practical purposes and to maintain continuity between Chapters 3 and 4). As more comprehensive earthquake catalogues become available it is recommended to extend the time window in which the statistical analysis is performed. 134 • Induced seismicity is not limited to hydraulic fracturing activities and experiences from other practices such as enhanced geothermal systems and CO2 sequestration can be used to improve understanding of induced seismicity. Experiences from enhanced geothermal systems are particularly of interest as they predominantly involve crystalline rock, providing a notable geological contrast to experiences in shale gas which involve sedimentary rocks but for which the distance to the crystalline basement appears to be an important factor. • The numerical modelling presented in this thesis were largely based on conceptual realizations dictated by the relationships observed in the empirical analyses. Further extension of this work would benefit from using site-specific data to properly constrain and calibrate the models. The use of microseismic data, for example, represents one avenue for further investigation of hydraulic fracturing behavior and model calibration, by simulating synthetic seismicity to compare to specific hydraulic fracturing performance and rock mass response in an actual shale gas reservoir. • Numerical modelling work could also be expanded to explore the use of newer 3-D hybrid FEM-DEM brittle fracturing techniques using commercial software like ELFEN or IRAZU. Research using these programs have been largely restricted to 2-D analyses, whereas 3-D considerations were required for the studies carried out in this thesis. However, as computing power continues to increase and the 3-D versions of these 135 software improve, they will offer a means to analyze the role of intact rock progressive failure and interactions between hydraulic fracturing and natural fracture networks. Other new 3-D hydraulic fracturing simulators, for example the lattice-model program XSite offer additional purpose-built options. • Further research is needed to understand the role of natural fracture networks and discontinuities on fluid flow and pressure transfer from well to fault. Again, new numerical modelling programs undergoing verification and validation offer one means to investigate this, together with the use of microseismic data for field-based insights and added validation. • Further research is required on how to discriminate between natural earthquakes and induced seismicity events. The outcome of this research has direct influence on the results and accuracy of statistical analyses. At present, the spatial coverage and surface-based installation of seismic monitoring networks limit the accuracy of detected event locations to epicentres on surface with significant errors being associated with the depth component. • Further research is needed to validate the results of machine learning analyses with state-of-the-art laboratory experiments and advanced 3-D discontinuum-based numerical modelling techniques to provide mechanistic understanding into the induced seismicity associations determined. This will help contribute to the application of machine learning algorithms as a tool for developing induced seismicity susceptibility maps. 136 References Akber Hassan, W. A., & Jiang, X. (2012). Upscaling and its application in numerical simulation of long-term CO2 storage. Greenhouse Gases: Science and Technology, 2(6), 408–418. https://doi.org/10.1002/ghg Allen, D. ., Eberhardt, E., & Bustin, A. (2019). Scientific Review of Hydraulic Fracturing in British Columbia (Issue February). Amini, A., & Eberhardt, E. (2019). Influence of tectonic stress regime on the magnitude distribution of induced seismicity events related to hydraulic fracturing. Journal of Petroleum Science and Engineering, 182(June 2019). https://doi.org/10.1016/j.petrol.2019.106284 Amitrano, D. (2003). Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b value. Journal of Geophysical Research, 108, 1–15. https://doi.org/10.1029/2001JB000680 Arthur, D. J., & Langhus, B. P. G. (2008). An overview of modern shale gas development in the United States. ALL Consulting, 1–21. http://www.all-llc.com/publicdownloads/ALLShaleOverviewFINAL.pdf ASTM International. (2004). Standard Test Method for Determination of the In-Situ Stress in Rock Using the Hydraulic Fracturing Method 1. Changes, i, 1–8. https://doi.org/10.1520/D4645-08 Atkinson, G. M., Eaton, D. W., Ghofrani, H., Walker, D., Cheadle, B., Schultz, R., Shcherbakov, R., Tiampo, K., Gu, J., Harrington, R. M., Liu, Y., Van Der Baan, M., & Kao, H. (2016). Hydraulic fracturing and seismicity in the western Canada sedimentary basin. Seismological Research Letters, 87(3), 631–647. https://doi.org/10.1785/0220150263 137 Atkinson, G. M., Eaton, D. W., & Igonin, N. (2020). Developments in understanding seismicity triggered by hydraulic fracturing. Nature Reviews Earth & Environment, 1(5), 264–277. https://doi.org/10.1038/s43017-020-0049-7 Bao, X., & Eaton, D. W. (2016). Fault activation by hydraulic fracturing in western Canada. Science, 354(6318), 1406–1409. https://doi.org/10.1126/science.aag2583 Barclay, J. E., Krause, F. F., Campbell, R. I., & Utting, J. (1990). Dynamic casting and growth faulting: Dawson Creek Graben Complex, Carboniferous-Permian Peace River Embayment, western Canada. Bulletin of Canadian Petroleum Geology, 38 A(40790), 115–145. BC Oil & Gas Commission. (2018). BC Oil & Gas Commission online library. http://www.bcogc.ca/ BC Oil and Gas Commision. (2012). Investigation of Observed Seismicity in the Horn River Basin. August, 29. BC Oil and Gas Commission. (2012). Montney Formation Play Atlas for NEBC. October, 36. BC Oil and Gas Commission. (2014). Investigation of Observed Seismicity in the Montney Trend. December, 32. Berger, Z. (1994). Satellite Hydrocarbon Exploration. Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-78587-0 Brudzinski, M. R., & Kozłowska, M. (2019). Seismicity induced by hydraulic fracturing and wastewater disposal in the Appalachian Basin, USA: a review. Acta Geophysica, 67(1), 351–364. https://doi.org/10.1007/s11600-019-00249-7 Cappa, F., & Rutqvist, J. (2011). Impact of CO2 geological sequestration on the nucleation of earthquakes. Geophysical Research Letters, 38(17), 2–7. https://doi.org/10.1029/2011GL048487 138 Cesca, S., Rohr, A., & Dahm, T. (2012). Discrimination of induced seismicity by full moment tensor inversion and decomposition. Journal of Seismology, 17(1), 147–163. https://doi.org/10.1007/s10950-012-9305-8 Cipolla, C. L., Lolon, E. P., Erdle, J. C., & Rubin, B. (2010). Reservoir modeling in shale-gas reservoirs. SPE Reservoir Evaluation and Engineering, 13(4), 638–653. https://doi.org/10.2118/125530-PA Cloos, H. (1928). Experimenten zur inneren Tektonik. Centralblatt fur Mineralogie and Paleontologie 1928B, 609. Cracknell, M. J., & Reading, A. M. (2014). Geological mapping using remote sensing data: A comparison of five machine learning algorithms, their response to variations in the spatial distribution of training data and the use of explicit spatial information. Computers and Geosciences, 63, 22–33. https://doi.org/10.1016/j.cageo.2013.10.008 Cundall, P. a., & Hart, R. D. (1992). Numerical Modelling of Discontinua. Engineering Computations, 9(2), 101–113. https://doi.org/10.1108/eb023851 Currie, B. S., Free, J. C., Brudzinski, M. R., Leveridge, M., & Skoumal, R. J. (2018). Seismicity Induced by Wastewater Injection in Washington County, Ohio: Influence of Preexisting Structure, Regional Stress Regime, and Well Operations. Journal of Geophysical Research: Solid Earth, 123(5), 4123–4140. https://doi.org/10.1002/2017JB015297 Davies, G. R., Watson, N., Moslow, T. F., & MacEachern, J. A. (2018). Regional subdivisions, sequences, correlations and facies relationships of the Lower Triassic Montney Formation, west-central Alberta to northeastern British Columbia, Canada — with emphasis on role of paleostructure. Bulletin of Canadian Petroleum Geology, 66(1), 23–92. Davies, R., Foulger, G., Bindley, A., & Styles, P. (2013). Induced seismicity and hydraulic 139 fracturing for the recovery of hydrocarbons. Marine and Petroleum Geology, 45, 171–185. https://doi.org/10.1016/j.marpetgeo.2013.03.016 Dixon, J. (2009). The Lower Triassic Shale member of the Montney Formation in the subsurface of northeast British Columbia. https://doi.org/10.4095/248223 Dixon, James. (2011). A review of the character and interpreted origins of thick, mudstone-encased sandstone bodies in the middle triassic doig formation of western canada. Bulletin of Canadian Petroleum Geology, 59(3), 261–276. https://doi.org/10.2113/gscpgbull.59.3.261 DOE/EIA. (2016). Annual Energy Outlook 2016. Dusseault, M. B., Mclennan, J., & Shu, J. (2011). Massive multi-stage hydraulic fracturing for oil and gas recovery from low mobility reservoirs in China. Petroleum Drilling Techniques, January. Earthquakes Canada, G. (2018). Earthquakes Canada, GSC,Earthquake Search (On-line Bulletin). http://earthquakescanada.nrcan.gc.ca/stndon/NEDB-BNDS/bull-eng.php Eaton, D. W., & Schultz, R. (2018). Increased likelihood of induced seismicity in highly overpressured shale formations. Geophysical Journal International, 214(1), 751–757. https://doi.org/10.1093/gji/ggy167 Eberhardt, E., & Amini, A. (2018). Hydraulic Fracturing. In Encyclopedia of Engineering Geology. Encyclopedia of Earth Sciences Series (pp. 1–6). Springer. https://doi.org/10.1007/978-3-319-12127-7_159-1 Edwards, Barclay, Gibson, Kvill, & Halton. (1994). Triassic Strata of the Western Canada Sedimentary Basin. Canadian Society of Petroleum Geologists and Alberta Research Council. https://ags.aer.ca/document/Atlas/chapter_16.pdf 140 El-Isa, Z. H., & Eaton, D. W. (2014). Spatiotemporal variations in the b-value of earthquake magnitude-frequency distributions: Classification and causes. Tectonophysics, 615–616, 1–11. https://doi.org/10.1016/j.tecto.2013.12.001 Ellsworth, W. (2013). Injection-Induced Earthquakes. Science. https://doi.org/10.1126/science.1225942 Ellsworth, W. L. (2013). Injection-Induced Earthquakes. Science, 341(July), 1–8. https://doi.org/10.1785/gssrl.83.2.250 Eyre, T. S., Eaton, D. W., Garagash, D. I., Zecevic, M., Venieri, M., Weir, R., & Lawton, D. C. (2019). The role of aseismic slip in hydraulic fracturing–induced seismicity. Science Advances, 5(8), 1–11. https://doi.org/10.1126/sciadv.aav7172 Farahbod, A. M., Kao, H., Cassidy, J. F., & Walker, D. (2015). How did hydraulic-fracturing operations in the Horn River Basin change seismicity patterns in northeastern British Columbia, Canada? The Leading Edge, 34(6), 658–663. https://doi.org/10.1190/tle34060658.1 Farahbod, A. M., Kao, H., Walker, D. M., Cassidy, J. F., & Calvert, A. (2015). Investigation of regional seismicity before and after hydraulic fracturing in the Horn River Basin, northeast British Columbia. Canadian Journal of Earth Sciences, 52(2), 112–122. https://doi.org/10.1139/cjes-2014-0162 Feng, Y. J., & Shi, X. W. (2012). Hydraulic Fracturing Process: Roles of In Situ Stress and Rock Strength. Advanced Materials Research, 616–618, 435–440. https://doi.org/10.4028/www.scientific.net/AMR.616-618.435 Fracfocous. (2017). National hydraulic fracturing chemical registry. http://fracfocus.org/ Galis, M., Ampuero, J. P., Mai, P. M., & Cappa, F. (2017). Induced seismicity provides insight 141 into why earthquake ruptures stop. Science Advances, 3(12). https://doi.org/10.1126/sciadv.aap7528 geoLOGIC systems ltd. (2019). geoSCOUT version 8.12; geo- LOGIC systems ltd., mapping, data management and analysis software. Gibbons, J. D., & Subhabrata, C. (2011). Nonparametric Statistical Inference. Gibson, & Barclay. (1989). Middle Absaroka Sequence, the Triassic stable Craton. In: Ricketts, B.D., Ed., Western Canada Sedimentary Basin: A Case History. Canadian Society of Petroleum Geologists. Göbel, T. (2015). A comparison of seismicity rates and fluid-injection operations in Oklahoma and California: Implications for crustal stresses. The Leading Edge, 34(6), 640–648. https://doi.org/10.1190/tle34060640.1 Guglielmi, Y., Cappa, F., Avouac, J.-P., Henry, P., & Elsworth, D. (2015). Seismicity triggered by fluid injection – induced aseismic slip. Science, 348(6240), 1224–1227. https://doi.org/10.1126/science.aab0476 Gulia, L., & Wiemer, S. (2010). The influence of tectonic regimes on the earthquake size distribution: A case study for Italy. Geophysical Research Letters, 37(10), 1–6. https://doi.org/10.1029/2010GL043066 Gupta, H. K., & Chadha, R. K. (1995). Induced Seismicity (H. K. Gupta & R. K. Chadha (Eds.)). Birkhäuser Basel. https://doi.org/10.1007/978-3-0348-9238-4 Gutenberg, B., & Richter, C. F. (1945). Seismicity of the earth. Bulletin of the Geological Society of America, 56(6), 603–667. https://doi.org/10.1130/0016-7606(1945)56[603:SOTE]2.0.CO;2 Hallo, M., Oprsal, I., Eisner, L., & Ali, M. Y. (2014). Prediction of magnitude of the largest 142 potentially induced seismic event. Journal of Seismology, 18(3), 421–431. https://doi.org/10.1007/s10950-014-9417-4 Hanks, T. C., & Kanamori, H. (1979). A moment magnitude scale. Journal of Geophysical Research B: Solid Earth, 84(B5), 2348–2350. https://doi.org/10.1029/JB084iB05p02348 Hasegawa, H. S., Wetmiller, R. J., & Gendzwill, D. J. (1989). Induced seismicity in mines in Canada-An overview. Pure and Applied Geophysics PAGEOPH, 129(3–4), 423–453. https://doi.org/10.1007/BF00874518 Healy, J. H., Rubey, W. W., Griggs, D. T., & Raleigh, C. B. (1968). The Denver Earthquakes. Science, 7(1), 79–92. Heidbach, O., Rajabi, M., Cui, X., Fuchs, K., Müller, B., Reinecker, J., Reiter, K., Tingay, M., Wenzel, F., Xie, F., Ziegler, M. O., Zoback, M. Lou, & Zoback, M. (2018). The World Stress Map database release 2016: Crustal stress pattern across scales. Tectonophysics, 744(July), 484–498. https://doi.org/10.1016/j.tecto.2018.07.007 Heidbach, O., Tingay, M., Barth, A., Reinecker, J., Kurfeß, D., & Müller, B. (2010). Global crustal stress pattern based on the World Stress Map database release 2008. Tectonophysics, 482(1–4), 3–15. https://doi.org/10.1016/j.tecto.2009.07.023 Herrmann, R. B., Benz, H., & Ammon, C. J. (2011). Monitoring the Earthquake source process in North America. Bulletin of the Seismological Society of America, 101(6), 2609–2625. https://doi.org/10.1785/0120110095 Hincks, T., Aspinall, W., Cooke, R., & Gernon, T. (2018). Oklahoma’s induced seismicity strongly linked to wastewater injection depth. Science, 359(6381), 1251–1255. https://doi.org/10.1126/science.aap7911 Holmgren, J. (2015). Induced Seismicity in the Dannemora Mine , Sweden Inducerad seismicitet 143 vid. 341. Holmgren, J. M., Atkinson, G. M., & Ghofrani, H. (2019). Stress drops and directivity of induced earthquakes in the western Canada sedimentary basin. Bulletin of the Seismological Society of America, 109(5), 1635–1652. https://doi.org/10.1785/0120190035 Horton, S. (2012a). Disposal of Hydrofracking Waste Fluid by Injection into Subsurface Aquiders Triggers Earthquake Swarm in Central Arkansas with Potential for Damaging Earthquakes. Seismological Research Letters. Horton, S. (2012b). Disposal of Hydrofracking Waste Fluid by Injection into Subsurface Aquifers Triggers Earthquake Swarm in Central Arkansas with Potential for Damaging Earthquake. Seismological Research Letters, 83(2), 250–260. https://doi.org/10.1785/gssrl.83.2.250 Hubbert, M. K., & Willis, D. G. (1957). Mechanics Of Hydraulic Fracturing. Transactions of the AIME, 210(01), 153–168. https://doi.org/10.2118/686-G Hummel, N., & Shapiro, S. a. (2013). Case History Nonlinear diffusion-based interpretation of induced microseismicity : A Barnett Shale hydraulic fracturing case study. Geophysics, 78(5), B211–B226. https://doi.org/10.1190/GEO2012-0242.1 Itasca Consulting Group Inc., 1999. (1999). UDEC/3DEC (Universal Distinct Element Code in 2/3 Dimensions), Version 6.0. Jeffrey, R. G., Bunger, a. P., Lecampion, B., Zhang, X., Chen, Z. R., van As, a., Allison, D. P., de Beer, W., Dudley, J. W., Siebrits, E., Thiercelin, M., & Mainguy, M. (2009). Measuring Hydraulic Fracture Growth in Naturally Fractured Rock. SPE Annual Technical Conference and Exhibition, 1–19. https://doi.org/10.2118/124919-MS Johnson, E. G., & Johnson, L. A. (2012). Hydraulic fracture water usage in northeast British 144 Columbia: locations, volumes and trends. In Geoscience Reports. Kao, H, Shan, S.-J., Bent, A., Woodgold, C., Rogers, G., Cassidy, J. F., & Ristau, J. (2012). Regional centroid-moment-tensor analysis for earthquakes in Canada and adjacent regions: an update. Seismological Research Letters, 83(3), 505–515. https://doi.org/10.1785/gssrl.83.3.505 Kao, Honn, Hyndman, R., Jiang, Y., Visser, R., Smith, B., Babaie Mahani, A., Leonard, L., Ghofrani, H., & He, J. (2018). Induced Seismicity in Western Canada Linked to Tectonic Strain Rate: Implications for Regional Seismic Hazard. Geophysical Research Letters, 45(20), 11,104-11,115. https://doi.org/10.1029/2018GL079288 Keranen, K. M., Weingarten, M., Abers, G. A., Bekins, B. A., & Ge, S. (2014). Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science, 345(6195). Kim, W. Y. (2013). Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio. Journal of Geophysical Research: Solid Earth, 118(7), 3506–3518. https://doi.org/10.1002/jgrb.50247 King, G. E. (2012). Hydraulic Fracturing 101: What Every Representative, Environmentalist, Regulator, Reporter, Investor, University Researcher, Neighbor and Engineer Should Know About Estimating Frac Risk and Improving Frac Performance in Unconventional Gas and Oil Wells. SPE Hydraulic Fracturing Technology Conference, 1–80. https://doi.org/10.2118/152596-MS Kraft, T., & Deichmann, N. (2014). High-precision relocation and focal mechanism of the injection-induced seismicity at the Basel EGS. Geothermics, 52, 59–73. https://doi.org/10.1016/j.geothermics.2014.05.014 145 Kwiatek, G., Martínez-Garzõn, P., Dresen, G., Bohnhoff, M., Sone, H., & Hartline, C. (2015). Effects of long-term fluid injection on induced seismicity parameters and maximum magnitude in northwestern part of the Geysers geothermal field. Journal of Geophysical Research: Solid Earth, 120(10), 7085–7101. https://doi.org/10.1002/2015JB012362 Langenbruch, C., & Zoback, A. M. D. (2016). How will induced seismicity in Oklahoma respond to decreased saltwater injection rates? Science Advances, 2(11), 1–10. https://doi.org/10.1126/sciadv.1601542 Lisjak, A., & Grasselli, G. (2014). A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. Journal of Rock Mechanics and Geotechnical Engineering, 6(4), 301–314. https://doi.org/10.1016/j.jrmge.2013.12.007 Liu, K., & Liu, B. (2017). Optimization of smooth blasting parameters for mountain tunnel construction with specified control indices based on a GA and ISVR coupling algorithm. Tunnelling and Underground Space Technology, 70, 363–374. https://doi.org/10.1016/j.tust.2017.09.007 Lund Snee, J.-E., & Zoback, M. D. (2016). State of stress in Texas: Implications for induced seismicity. Geophysical Research Letters, 43, 10,208–10,214. https://doi.org/10.1002/2016GL070974 Lundberg, S., & Lee, S.-I. (2017). A Unified Approach to Interpreting Model Predictions Scott. Nips. Mahani, A. B., Esfahani, F., Kao, H., Gaucher, M., Hayes, M., Visser, R., & Venables, S. (2019). A systematic study of earthquake source mechanism and regional stress field in the southern Montney unconventional play of northeast British Columbia, Canada. Seismological Research Letters, 91(1), 195–206. https://doi.org/10.1785/0220190230 146 Mahani, A. B., Schultz, R., Kao, H., Walker, D., Johnson, J., & Salas, C. (2017). Fluid injection and seismic activity in the Northern Montney Play, British Columbia, Canada, with special reference to the 17 August 2015 Mw 4.6 induced earthquake. Bulletin of the Seismological Society of America, 107(2), 542–552. https://doi.org/10.1785/0120160175 Majer, E. L., Baria, R., Stark, M., Oates, S., Bommer, J., Smith, B., & Asanuma, H. (2007). Induced seismicity associated with Enhanced Geothermal Systems. Geothermics, 36(3), 185–222. https://doi.org/10.1016/j.geothermics.2007.03.003 Massey, F. J. (1951). The Kolmogorov-Smirnov Test for Goodness of Fit. Journal of the American Statistical Association, 46(253), 68. https://doi.org/10.2307/2280095 Maxwell, S. C., Urbancic, T. I., Steinsberger, N., Energy, D., & Zinno, R. (2002). Microseismic Imaging of Hydraulic Fracture Complexity in the Barnett Shale. SPE International, 1–9. https://doi.org/10.2523/77440-MS Mazzoldi, A., Rinaldi, A. P., Borgia, A., & Rutqvist, J. (2012). Induced seismicity within geological carbon sequestration projects: Maximum earthquake magnitude and leakage potential from undetected faults. International Journal of Greenhouse Gas Control, 10, 434–442. https://doi.org/10.1016/j.ijggc.2012.07.012 McGarr, A. (2014). Maximum magnitude earthquakes induced by fluid injection. Journal of Geophysical Research: Solid Earth, 1–12. https://doi.org/10.1002/2013JB010597.Received McNamara, D. E., Benz, H. M., Herrmann, R. B., Bergman, E. A., Earle, P., Holland, A., Baldwin, R., & Gassner, A. (2015). Earthquake hypocenters and focal mechanisms in central Oklahoma reveal a complex system of reactivated subsurface strike-slip faulting. Geophysical Research Letters, 42(8), 2742–2749. https://doi.org/10.1002/2014GL062730 Millar, D., & Clarici, E. (2002). Investigation of back-propagation artificial neural networks in 147 modelling the stress-strain behaviour of sandstone rock. Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94), 5, 3326–3331. https://doi.org/10.1109/ICNN.1994.374770 Mori, J., & Abercrombie, R. E. (1997). Depth dependence of earthquake frequency-magnitude distributions in California: Implications for rupture initiation. Journal of Geophysical Research-Solid Earth, 102(B7), 15081–15090. https://doi.org/Doi 10.1029/97jb01356 OCC. (2018). Oklahoma Corporation Commission. Ouenes, A., Aimene, Y. E., & Nairn, J. (2013). Interpretation of Microseismic Using Geomechanical Modeling of Multiple Hydraulic Fractures Interacting with Natural Fractures – Application to Montney Shale. Canadian Society of Exploration Geophysists, 46–54. Panakkat, A., & Adeli, H. (2009). Recurrent Neural Network for Approximate Earthquake Time and Location Prediction Using Multiple Seismicity Indicators. Computer-Aided Civil and Infrastructure Engineering, 24(4), 280–292. https://doi.org/10.1111/j.1467-8667.2009.00595.x Pawley, S., Schultz, R., Playter, T., C., H., S., T., L., & Hauck, T. (2018). The geological susceptibility of induced earthquakes in the Duvernay play. Geophysical Research Letters, 45, 1786–1793. https://doi.org/https://doi.org/10.1002/ 2017GL076100 Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., & Duchesnay, ´Edouard. (2011). Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12, 2825–2830. https://doi.org/10.1145/2786984.2786995 148 Pirayehgar, A., & Dusseault, M. B. (2015). NUMERICAL INVESTIGATION OF SEISMIC EVENTS ASSOCIATED WITH HYDRAULIC FRACTURING. ISRM International Symposium. Montreal, Canada. Preisig, G., Eberhardt, E., Gischig, V., Roche, V., van der Baan, M., Valley, B., Kaiser, P. K., Duff, D., & Lowther, R. (2015). Development of connected permeability in massive crystalline rocks through hydraulic fracture propagation and shearing accompanying fluid injection. Geofluids, 15(1–2), 321–337. https://doi.org/10.1111/gfl.12097 Pu, Y., Apel, D. B., & Lingga, B. (2018). Rockburst prediction in kimberlite using decision tree with incomplete data. Journal of Sustainable Mining, 17(3), 158–165. https://doi.org/10.1016/j.jsm.2018.07.004 Ribeiro e Sousa, L., Miranda, T., Leal e Sousa, R., & Tinoco, J. (2017). The Use of Data Mining Techniques in Rockburst Risk Assessment. Engineering, 3(4), 552–558. https://doi.org/10.1016/J.ENG.2017.04.002 Riedel, W. (1929). Zur mechanik geologischer brucherscheinungen. Centralblatt fur Minerologie, Geologie, und Paleontologie 1929B, 354. Ristau, J., Rogers, G. C., & Cassidy, J. F. (2007). Stress in western Canada from regional moment tensor analysis. Canadian Journal of Earth Sciences, 44(2), 127–148. https://doi.org/10.1139/e06-057 Rogers, S., & McLellan, P. (2014). Investigation of the Effects of Natural Fractures and Faults on Hydraulic Fracturing in the Montney Formation , Farrell Creek Gas Field , British Columbia. International Discrete Fracture Network Engineering Conference, October, 1–13. Rouet-Leduc, B., Hulbert, C., Lubbers, N., Barros, K., Humphreys, C. J., & Johnson, P. A. 149 (2017). Machine Learning Predicts Laboratory Earthquakes. Geophysical Research Letters, 44(18), 9276–9282. https://doi.org/10.1002/2017GL074677 Rutqvist, J., Rinaldi, A. P., Cappa, F., & Moridis, G. J. (2013). Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs. Journal of Petroleum Science and Engineering, 107, 31–44. https://doi.org/10.1016/j.petrol.2013.04.023 Rutqvist, J., Rinaldi, A. P., Cappa, F., & Moridis, G. J. (2015). Modeling of fault activation and seismicity by injection directly into a fault zone associated with hydraulic fracturing of shale-gas reservoirs. Journal of Petroleum Science and Engineering, 127(January 2011), 377–386. https://doi.org/10.1016/j.petrol.2015.01.019 Sainoki, A., & Mitri, H. (2014). Estimation of Mining-induced Fault-slip Potential in Underground Mines. October, 1491–1500. Scholz, C. H. (1968). The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. Bulletin of the Seismological Society of America, 58(1), 399–415. http://www.bssaonline.org/cgi/content/abstract/58/1/399 Schorlemmer, D., Wiemer, S., & Wyss, M. (2005). Variations in earthquake-size distribution across different stress regimes. Nature, 437(7058), 539–542. https://doi.org/10.1038/nature04094 Schultz, R., Atkinson, G., Eaton, D. W., Gu, Y. J., & Kao, H. (2018). Hydraulic fracturing volume is associated with induced earthquake productivity in the duvernay play. Science, 359(6373), 304–308. https://doi.org/10.1126/science.aao0159 Schultz, Ryan, Corlett, H., Haug, K., Kocon, K., Maccormack, K., Stern, V., & Shipman, T. (2016). Linking fossil reefs with earthquakes: Geologic insight to where induced seismicity 150 occurs in Alberta. Geophysical Research Letters, 43, 2534–2542. https://doi.org/10.1002/2015GL067514.Received Schultz, Ryan, Skoumal, R. J., Brudzinski, M. R., Eaton, D., Baptie, B., & Ellsworth, W. (2020). Hydraulic fracturing‐induced seismicity. Reviews of Geophysics, 58(3), 1–43. https://doi.org/10.1029/2019RG000695 Schultz, Ryan, Stern, V., & Gu, Y. J. (2014). An investigation of seismicity clustered near the Cordel Field, west central Alberta, and its relation to a nearby disposal well. Journal of Geophysical Research: Solid Earth, 119(4), 3410–3423. https://doi.org/10.1002/2013JB010836 Schultz, Ryan, Wang, R., Gu, Y. J., Haug, K., & Atkinson, G. (2017). A seismological overview of the induced earthquakes in the Duvernay play near Fox Creek, Alberta. Journal of Geophysical Research: Solid Earth, 122(1), 492–505. https://doi.org/10.1002/2016JB013570 Shah, A. K., & Keller, G. R. (2017). Geologic influence on induced seismicity: Constraints from potential field data in Oklahoma. Geophysical Research Letters, 44(1), 152–161. https://doi.org/10.1002/2016GL071808 Shapiro, S. A., & Dinske, C. (2009). Fluid-induced seismicity: Pressure diffusion and hydraulic fracturing. Geophysical Prospecting, 57(2), 301–310. https://doi.org/10.1111/j.1365-2478.2008.00770.x Skoumal, R. J., Brudzinski, M. R., & Currie, B. S. (2015). Earthquakes induced by hydraulic fracturing in Poland township, Ohio. Bulletin of the Seismological Society of America, 105(1), 189–197. https://doi.org/10.1785/0120140168 Smith, M.B., Montgomety, C. T. (2015). Hydraulic Fracturing. In Reference Module in 151 Chemistry, Molecular Sciences and Chemical Engineering. https://doi.org/10.1016/B978-0-12-409547-2.11013-3 Soeder, D. J. (1988). Porosity and Permeability of Eastern Devonian Gas Shale. SPE Formation Evaluation, 3(01), 116–124. https://doi.org/10.2118/15213-PA Spada, M., Tormann, T., Wiemer, S., & Enescu, B. (2013). Generic dependence of the frequency-size distribution of earthquakes on depth and its relation to the strength profile of the crust. Geophysical Research Letters, 40(4), 709–714. https://doi.org/10.1029/2012GL054198 Sumy, D. F., Cochran, E. S., Keranen, K. M., Wei, M., & Abers, G. a. (2014). Observations of static Coulomb stress triggering of the November 2011 M 5.7 Oklahoma earthquake sequence. Journal of Geophysical Research: Solid Earth, 119(3), 1904–1923. https://doi.org/10.1002/2013JB010612.Received Sun, Y., Feng, X., & Yang, L. (2018). Predicting Tunnel Squeezing Using Multiclass Support Vector Machines. Advances in Civil Engineering, 2018, 17–20. https://doi.org/10.1155/2018/4543984 Tagliabue, J. (2013). Parts of Low Country Are Now Quake Country. The New York Times, 2–4. The Ultimate Potential for Unconventional Petroleum from the Montney Formation of British Columbia and Alberta - Energy Briefing Note. (2013). Uddameri, V., Morse, A., & Tindle, K. J. (2015). Hydraulic Fracturing Impacts and Technologies: A Multidisciplinary Perspective. https://doi.org/10.1016/j.jngse.2015.12.020 USGS. (2018). USGS Earthquake Cataloge. https://earthquake.usgs.gov/earthquakes/search/ van As, A., & Jeffrey, R. G. (2000). Hydraulic fracturing as a cave inducement technique at northparkes mines. MassMin 2000 Proceedings. 152 Van der Baan, M., & Calixto, F. J. (2017). Human-induced seismicity and large-scale hydrocarbon production in the USA and Canada. Geochemistry, Geophysics, Geosystems, 18(7), 2467–2485. https://doi.org/10.1002/2017GC006915 Van Der Baan, M., Eaton, D., & Dusseault, M. (2013). Microseismic Monitoring Developments in Hydraulic Fracture Stimulation. Effective and Sustainable Hydraulic Fracturing, 439–466. https://doi.org/10.5772/45724 van der Elst, N. J., Page, M. T., Weiser, D. A., Goebel, T. H. W., & Hosseini, S. M. (2016). Induced earthquake magnitudes are as large as (statistically) expected. Journal of Geophysical Research: Solid Earth, 121(6), 4575–4590. https://doi.org/10.1002/2016JB012818 Vermylen, J., & Zoback, M. (2011). Hydraulic fracturing, microseismic magnitudes, and stress evolution in the Barnett Shale, Texas, USA. SPE Hydraulic Fracturing Technology …, 2009, SPE 140507. https://doi.org/10.2118/140507-MS Vishkai, M., Wang, J., Wong, R. C. K., Clarkson, C. R., & Gates, I. D. (2017). Modeling geomechanical properties in the montney formation, Alberta, Canada. International Journal of Rock Mechanics and Mining Sciences, 96(May), 94–105. https://doi.org/10.1016/j.ijrmms.2017.04.001 Visser, R., Smith, B., Kao, H., Babaie Mahani, A., Hutchinson, J., & McKay, J. E. (2017). A comprehensive earthquake catalogue for northeastern British Columbia and western Alberta, 2014-2016. Geological Survey of Canada, Open File, 8335(February 2018), 1–28. https://doi.org/10.4095/306292 Walsh, F. R., & Zoback, M. D. (2015). Oklahoma’s recent earthquakes and saltwater disposal. Science Advances, 1(5), e1500195–e1500195. https://doi.org/10.1126/sciadv.1500195 153 Warpinski, N. R., Du, J., & Zimmer, U. (2012). Measurements of Hydraulic-Fracture-Induced Seismicity in Gas Shales. SPE Hydraulic Fracturing Technology Conference, February, 6–8. https://doi.org/10.2118/151597-PA WEC. (2013). World Energy Resources. In World Energy Council. Weidler, R., Gérard, A., Baria, R., Baumgärtner, J., & Jung, R. (2002). Hydraulic and micro-seismic results of a massive stimulation test at 5 km depth at the European Hot-Dry-Rock test site, Soultz, France. Twenty-Seventh Workshop on Geothermal Reservoir Engineering, 27, 95–100. Weingarten, M., Ge, S., Godt, J. W., Bekins, B. a, & Rubinstein, J. L. (2015). High-rate injection is associated with the increase in U.S. mid-continent seismicity. Science, 348(6241), 1336–1340. https://doi.org/10.1126/science.aab1345 Westerhaus, M., Wyss, M., Yilmaz, R., & Zschau, J. (2002). Correlating variations of b values and crustal deformations during the 1990s may have pinpointed the rupture initiation of the Mw = 7.4 Izmit earthquake of 1999 August 17. Geophysical Journal International, 148(1), 139–152. https://doi.org/10.1046/j.0956-540x.2001.01554.x Wiemer, S., & Wyss, M. (2000). Minimum Magnitude of Completeness in Earthquake Catalogs: Examples from Alaska, the Western United States, and Japan. Bulletin of the Seismological Society of America, 90(4), 859–869. https://doi.org/10.1785/0119990114 Williamson, W. H., & Woolley, D. R. (1980). Hydraulic fracturing to improve the yield of bores in fractured rock. Australian Government Publishing Service, Canberra. Woessner, J., & Wiemer, S. (2005). Assessing the quality of earthquake catalogues: Estimating the magnitude of completeness and its uncertainty. Bulletin of the Seismological Society of America, 95(2), 684–698. https://doi.org/10.1785/0120040007 154 World Energy Scenarios. (2016). In World Energy Council. https://doi.org/10.2307/j.ctv6cfqsc.15 Zang, A., Oye, V., Jousset, P., Deichmann, N., Gritto, R., McGarr, A., Majer, E., & Bruhn, D. (2014). Analysis of induced seismicity in geothermal reservoirs - An overview. Geothermics, 52, 6–21. https://doi.org/10.1016/j.geothermics.2014.06.005 Zang, A., Yoon, J. S., Stephansson, O., & Heidbach, O. (2013). Fatigue hydraulic fracturing by cyclic reservoir treatment enhances permeability and reduces induced seismicity. Geophysical Journal International, 195(2), 1282–1287. https://doi.org/10.1093/gji/ggt301 Zangeneh, N., Eberhardt, E., & Bustin, R. M. (2015). Investigation of the influence of stress shadows on horizontal hydraulic fractures from adjacent lateral wells. Journal of Unconventional Oil and Gas Resources, 9, 54–64. https://doi.org/10.1016/j.juogr.2014.11.001 Zangeneh, Neda, Eberhardt, E., & Bustin, M. (2013). A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing. Effective and Sustainable Hydraulic Fracturing, 12. https://doi.org/http://dx.doi.org/10.5772/56191 Zangeneh, Neda, Eberhardt, E., & Bustin, R. M. (2015). Investigation of the influence of natural fractures and in situ stress on hydraulic fracture propagation using a distinct-element approach. Canadian Geotechnical Journal, 52(7), 926–946. https://doi.org/10.1139/cgj-2013-0366 Zhang, X., Zhang, J., Yuan, C., Liu, S., Chen, Z., & Li, W. (2020). Locating induced earthquakes with a network of seismic stations in Oklahoma via a deep learning method. Scientific Reports, 10(1), 1–12. https://doi.org/10.1038/s41598-020-58908-5 Zoback Mark D., & Arent, D. J. (2014). Shale Gas Development Opportunities and Challenges. 155 3–7. 156 Appendix A: Data and Resources The data used in Chapter 3 has been made available on the University of British Columbia’s digital repository (cIRcle). 157 Appendix B: Feature Definitions The following are the descriptions and details of the operational and geological features used for the machine learning analyses presented in Chapter 4. Each feature is listed with a short description, the type of feature (categorical or numerical), and the number of missing data entries where applicable. Operational Parameters Average Injection Rate (m3/min) (numerical) Average injection rate for each stage. 1473 missing entries for which the stages are removed from analysis. Injection Depth (m) (numerical) True vertical depth of each hydraulic fracturing stage. Stage Injection Volume (m3) Injection volume of each stage. 342 missing entries for which the stages are removed from analysis. Maximum Injection Pressure (MPa) (numerical) Maximum injection pressure of each stage. 3057 missing entries for which the stages are removed from analysis. 158 Local SHmax and Well Direction Azimuth Difference (deg) (numerical) The azimuth difference between the orientation of horizontal section of the well and maximum horizontal stress orientation at well. Well Completion Length (m) (numerical) The distance between first and last hydraulic fracturing stages at each well. Geological Parameters Pore Pressure Gradient (KPa/m) (numerical) Pore pressure gradient of Montney at each well calculated from interpolation of 2252 available data points in BCOGC’s database. Local and Regional SHmax Azimuth Difference (deg) (numerical) Azimuth difference between regional (assumed 45 degrees here) and local (interpolated at each well) maximum horizontal stress at each well. Distance to Faults (m) (numerical) The distance between each well and the closest mapped fault in Montney. No cut off value is considered. Distance to top of Montney (m) (numerical) Vertical distance between each stage injection depth and the interpolated top of Montney (TVD). 159 Distance to top of Debolt (m) (numerical) Vertical distance between each stage injection depth and the interpolated top of Debolt (TVD). Distance to Basement (m) (numerical) Vertical distance between each stage injection depth and the interpolated top of Basement (TVD).