Chapter 6

Simulating hydraulic fracturing using Finite-Discrete Element Method (FDEM): Effects of pre-existing joints and induced seismicity Murad AbuAisha a,b , David W. Eaton a , Jeffrey Priest b and Ron Wong b a Dept. of Geoscience, Univ. of Calgary, Calgary, AB, T2N 1N4, Canada. b Dept. of Civil Engineering, Univ. of Calgary, Calgary, AB, T2N 1N4, Canada. E: murad.abuaisha@ucalgary.ca

Summary

This numerical study investigates near-wellbore behaviour of Hydraulic Fracture (HF) growth using a Finite-Discrete Element Modeling (FDEM) approach. Our results rarely exhibit a simple bi-wing fracture geometry, but are consistent with analytical solutions and exhibit many features that have been documented in the field. In the case of a homogeneous medium, our simulations show that isotropic far-field stresses (Shmax = Shmin) promote the development of complex fractures, whereas anisotropic far-field stresses (Shmax Shmin) promote branching and curving fracture growth in the general direction of Shmax. A set of pre-existing randomly space distributed joints around the wellbore indicate that fractures prefer to initiate at the joints' tips once intersected by a fluid-driven fracture. They also reveal that possible seismic events due to formation deformation are induced because of shear slippage of critically stressed joints. The presence of multiple joints around the wellbore also increases the extent to which fluid-driven fractures can grow for the same injection energy. Introducing single isolated joints that do not intersect the wellbore, usually created by previous fracturing tests, imposes a lateral stress gradient. The presence of such a gradient leads to asymmetric fracture initiation, generally away from the pre-existing joint. The asymmetric fracture growth is deemed to cause a symmetric micro-seismicity around the wellbore, which explains some of the seismic phenomena observed in the fields.

Introduction

Hydraulic fracturing (HF) is a wellbore stimulation method used extensively in the development of unconventional oil and gas reservoirs to enhance connectivity of low-permeability (< 1 mD) oil and gas-bearing rock formations [START_REF] Nagel | Discrete element hydraulic fracture modeling -Evaluationg changes in natural fracture aperture and transmissivity: the Canadian Unconventional Resources Conference[END_REF]; [START_REF] King | Thirty years of gas shale fracturing: What have we learnt?: the SPE Annual Technical Conference and Exhibition[END_REF]; [START_REF] Hubbert | Mechanics of Hydraulic fracturing[END_REF]; [START_REF] Fjaer | Petroleum related rock mechanics[END_REF]. Injection of fluid under high pressure overcomes the rock tensile strength to create and/or extend fractures AbuAisha (2014); [START_REF] Bruel | Modelling heat extraction from forced fluid flow through simulated fractured rock masses: Evaluation of the Soultz-Sous-Forêts site potential[END_REF]; Barree et al. (2002). Through the creation of high-permeability pathways, HF serves to increase reservoir drainage by effectively increasing wellbore radius and thus reducing the hydrocarbon flow distance. HF has also been a controversial issue, especially in the case of naturally fractured formations for which concerns have been expressed about the potential for groundwater contamination [START_REF] Osborn | Modelling heat extraction from forced fluid flow through simulated fractured rock masses: Evaluation of the Soultz-Sous-Forêts site potential: Proceedings of the[END_REF]. Consequently, robust methods for simulating HF growth are critical.

Simple solutions of stress distribution around a borehole provide a starting point for analysis of HF operations [START_REF] Hubbert | Mechanics of Hydraulic fracturing[END_REF]; [START_REF] Fjaer | Petroleum related rock mechanics[END_REF]. In the simplest cases, these solutions are based upon assumptions of elastic, isotropic and homogeneous media [START_REF] Adachia | Computer simulation of hydraulic fractures[END_REF]Dusseault (2013). However, the presence of pre-existing fractures and weak bedding planes render such assumptions questionable, motivating the use of a more general numerical simulation approach. In this context a new technique that combines both Discrete Element Method (DEM) [START_REF] Munjiza | Combined single and smeared crack model in combined finite-discrete element analysis: International[END_REF] and Finite Element Method (FEM) [START_REF] Hughes | The finite element method; linear static and dynamic finite element analysis[END_REF] is well suited to this problem. FDEM methods have capabilities for simulating fluid-driven fractures in jointed/pre-fractured rock masses. This approach enables investigation of the transition from continuum to discontinuum behaviour through the use of FEM to simulate continuum behaviour of the intact material and DEM to capture interactions between fractured blocks [START_REF] Munjiza | The Combined Finite-Discrete Element Method[END_REF]; [START_REF] Mahabadi | Y-Geo: a new combined finite-discrete element numerical code for geomechanical applications[END_REF]; Lisjak et al. (2013).

This paper is organized as follows. First, two HF tests are simulated in an elastic continuous medium, considering both isotropic (Shmax = Shmin) and anisotropic (Shmax Shmin) far-field stress states. Next, tests are performed while considering a heavily fractured formation close to the wellbore. Thirdly, the presence of lateral stress gradients is investigated by introducing single isolated joints in the medium that are offset from the wellbore. Finally, the FDEM-HF approach is used to simulate a real DFIT test performed at the Montney formation in British Columbia/Canada.

The FDEM: Overview

The FDEM approach was first suggested by [START_REF] Munjiza | A combined finite-discrete element method in transient dynamics of fracturing solids[END_REF] and [START_REF] Munjiza | Penalty Pre-existing joints and HF function method for combined finite-discrete element systems comprising large number of separate bodies[END_REF]. It is a hybrid technique that combines the advantages of the FEM and DEM methods. While the medium is undergoing elastic deformation, the behaviour of intact material is explicitly modelled by the FEM method. When fractures are initiated, the interaction between discontinuous blocks is then captured by the DEM. The FDEM is capable of tracking crack initiation and propagation by applying the principles of non-linear elastic fracture mechanics [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses[END_REF][START_REF] Barenblatt | The mathematical theory of equlibrium cracks in brittle fracture[END_REF]. In this study we use an implementation of the FDEM method, known as Y-Geo [START_REF] Mahabadi | Y-Geo: a new combined finite-discrete element numerical code for geomechanical applications[END_REF]; [START_REF] Lisjak | Continuum-discontinuum analysis of failure mechanisms around unsupported circular excavations in anisotropic clay shales[END_REF], which consists of three main computational modules that exchange information at every time step (Fig. 1):

1. a geomechanics solver based on the FDEM approach, which captures the mechanical response of the rock mass (i.e. deformation and fracturing); 2. a cavity volume calculator, which dynamically tracks the evolution of wet fractures within the model and computes variations in cavity volume due to elastic deformation, fracturing, and fluid compressibility;

3. a pump model, which computes the fluid pressure based on the injection flow rate and cavity volume.

Solid solver (FDEM)

Elastic deformation Rock fractures Rock joints

Cavity volume calculator

Wet boundary tracking F r a c t u r e t o p o l o g y

Fluid compressibility model

Cavity volume, V F l u i d p r e s s u r e , p The Y-Geo code uses a FEM mesh of elastic triangular Delaunary elements1 with nondimensional rectangular cohesive fracture elements embedded between the triangular ones. While the medium is undergoing elastic loading, the fracture elements are assigned large penalty coefficients to eliminate them from the elasticity matrix. As soon as the tensile and/or the shear strengths of the material are reached, the material starts undergoing inelastic deformation which is expressed in terms of fracture energy dissipation (Fig. 2(a)). Once the fracture energies, G Ic of mode I and G IIc of mode II, are dissipated, the fracture elements are removed and fractures are initiated. At this point the positions of the separated blocks are tracked by the DEM.

Pumping boundary conditions

Solid solver of the Y-Geo code

The Y-Geo simulator is capable of tracking tensile (mode I) as well as shearing (mode II) of fracture growth. Hybrid fracturing can also be simulated as a combination of both modes (mode I-II). The fracturing mode is triggered based on the relative displacement of adjacent elastic triangular elements containing the fracture element. Similar to the cohesive model originally proposed for concrete [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF], a mode I crack initiates when the crack tip opening (δ) reaches a critical value (δ p ) which is related to the cohesive tensile strength of the rock (f t ) (Fig. 2(a)). As the fracture propagates and the crack tip opening increases, the normal bonding stress (σ) is assumed to decrease until a residual opening value (δ c ) is reached and a traction-free surface is created. Mode II fracturing is simulated by a slipweakening model Ida (1976). The shear bonding stress (τ ) is a function of the amount of slip (s) and the normal stress on the fracture (σ n ) (Fig. 2(a)). The critical slip (s p ) occurs when the shear bonding stress exceeds the cohesive shear strength of the rock (f s ) defined as,

f s = c + σ n tan(φ i ), (6.1)
where c is the internal cohesion, φ i is the material internal friction angle. Upon undergoing critical slip (s p ), the shear strength is gradually reduced to a residual value (f r ) which is a pure frictional resistance,

f r = σ n tan(φ f ), (6.2)
where φ f is the fracture friction angle.

In the mixed mode I-II of fracturing (Fig. 2(b)), the coupling between crack opening and slip is defined by an elliptical relationship where the failure envelope is defined as,

δ -δ p δ c -δ p 2 + s -s p s r -s p 2 ≥ 1. (6.3)
Pre-existing defects/rock joints are modelled as non-dimensional planes that can sustain two types of repulsive and shear forces. Repulsive forces are calculated using the normal penalty coefficient p n (Fig. 2(a)). These forces prevent element penetrability and can affect the stress conditions in the area close to the rock joint. Shear forces are simply calculated based on Eq. 6.2. If at any point the rock joint is intersected by a fluid-driven fracture, the fluid pressure percolation is distributed evenly over the entire discontinuity. In this study, rock joints are applied in the model initially, i.e. before geostatic stresses are applied.

Pressure-driven fracturing

Fluid injection and pressure-driven fracture propagation are captured by a simplified approach based on the principle of mass conservation for a compressible fluid injected into a deformable solid. The model is hydromechanically coupled in the sense that variations in cavity volume, due to either rock elastic deformation or fracturing, affect the pressure of the compressible fluid, which in turn, affects rock deformation and failure (Fig. 1). Fluid leakoff into the rock matrix is neglected, as the rocks of interest have exceptionally low matrix permeability and the simulations only capture the initial onset of fracture growth. The fluid-driven fracturing in the Y-Geo code is simulated as follows (Fig. 3): 1-the fluid injection induced pressure is assigned to the elements edges in contact with the borehole (initial wet boundaries); 2-as injection continues, the increasing pressure causes new fractures to be initiated; 3-the topology of the newly created fractures as well as connectivity with the previous wet boundaries are calculated by the geomechanics/solid solver; 4-by knowing the new topology and connectivity, cavity volume can be computed; 5-a fluid compressibility model is used to calculate fluid pressure on the new wet boundaries while considering any incoming fluid mass.

Mathematically, Green's theorem is used to calculate cavity volume by evaluating the following integral over all wet boundaries,

V = 1 2 x dy -y dx. (6.4)
Numerically, this integration is calculated as,

V = 1 2 i x i y i+1 -y i x i+1 , (6.5)
with i indexing over the nodes of the wet boundaries with coordinates (x i , y i ).

Once the cavity volume is computed, the fluid compressibility model is implemented to calculate fluid pressure on all the newly created wet boundaries,

K f = -V dp dρ f , (6.6)
with K f being fluid bulk modulus and ρ f being fluid density. While considering the amount of injected fluid at each time step, the fluid mass (m) is integrated from the flow rate,

p = p 0 + K f log m V ρ f 0 , (6.7)
where p 0 and ρ f 0 are the reference fluid pressure and density respectively. Based on the length and orientation of each wet element edge, the fluid pressure is converted into equivalent nodal forces and substituted for in the geomechanics solver.

Simulations of HF in elastic homogeneous media

In this section we consider fluid-driven fracturing in a homogeneous elastic medium subject to both isotropic and anisotropic far-field stress states. Fracture patterns are analysed and the values of fracturing pressures are compared to analytical solutions from Fjaer et al. [START_REF] Fjaer | Petroleum related rock mechanics[END_REF]. Similar works could have been done by other researchers, however, we just use them as an introduction for proceeding sections.

Geometry, meshing and material properties

The Boundary Value Problem (BVP) considered here is representative of a small-scale fracture test, such as a Diagnostic Fracture Injection Test (DFIT) [START_REF] Hawkes | Interpretation of closure pressure in the unconventional montney using PTA techniques: the SPE Hydraulic Fracturing Technology Conference[END_REF]. It uses a vertical wellbore of 0.1 m diameter located in 8 m × 8 m rock formation. To capture deformation process changes near the wellbore, the mesh is intensively refined with 0.003 m elements in this area. Away from the zone of the intensive refinement, element size is gradually increased to 0.3 m (Fig. 4).

Zone The material properties pertaining to rock formation and needed in the Y-Geo code are listed in Table 1. The input data as well as the in situ stress measurements correspond to a granitic formation of 2500 m depth Evans et al. (1972). Although the geomechanical characteristics of granite differ somewhat from shale formations, the main focus of this work is to in-vestigate fluid-driven fracture topology due to pre-existing Joints. Therefore, the specific details of the particular unit are less important. The stress values displayed later represent the Terzaghi's effective stresses, i.e. stress values above the formation fluid pressure. The sign convention adopted in this research attributes negative sign to compressive stresses.

Boundary and initial conditions

To study the effect of the far-field stress state, HF simulations are performed while assuming: 1-isotropic far-field stresses of σ h = σ H = -4.85 MPa and; 2-anisotropic far-field stresses of σ h = -4.9 MPa and σ H = -7.05 MPa. As effective stresses are used directly in the calculations, the values of fluid pressure represent the effective fluid pressure (the increment above the formation fluid pressure). A geostatic/elastic step is initially simulated without the wellbore being present by applying the initial far field stresses on the boundary of the model and allowing the forces to come into equilibrium throughout the model (zero velocity). The boundary conditions are then set to zero-displacement to allow for borehole excavation. The excavation of the borehole is implemented in two steps: the first step included reducing the elastic modulus of the borehole material by 10%; second step included removing the elements forming the borehole material. The final stage of the simulation included injecting fluid into the borehole at constant flow rate of Q = 65 l/s.

Simulations and discussion

Figure 5 shows the numerical treatment pressure profiles for the cases of isotropic and anisotropic stress states. The deviatoric part of the stress tensor facilitates fracturing by combining mode I and mode II [START_REF] Atkinson | Fracture mechanics of rock[END_REF]; [START_REF] Bieniawski | Mechanism of brittle fracture of rock: Part I-theory of the frac-ture process[END_REF]. The threshold of fracturing pressure is reduced by ∼10% for stress anisotropy of ∼27%. The analytical solution for stress distribution around a wellbore can be used to calculate a mathematical expression for the threshold of fracturing pressure. In terms of total stresses, Fjaer et al. [START_REF] Fjaer | Petroleum related rock mechanics[END_REF] found the following expression for borehole pressure at failure in the presence of mud cake2 :

p HF = -[(3σ h -σ H ) + p -f t ],
(6.8)

with p being the far-field formation pressure.

In our simulations the fracturing pressure for the cases of isotropic and anisotropic far-field stresses are 12.74 and 11.47 MPa respectively, while the analytical values derived using Eq. 6.8 are 12.7 and 10.5 MPa respectively. The discrepancy between the numerical and analytical values is related to the fact that in the numerical solution both mode I and mode II fractures are considered with a fracture energy term applied (G Ic and G IIc ) that accounts for yielding of the rock once the strength limit of the rock is exceeded (plastic zone after failure). In the analytical solution only Mode I fractures are considered and no fracture energy term is applied. This discrepancy is more apparent for the case of anisotropic far-field stresses which appear to induce more mode II fractures3 . The numerical response of treatment pressure also corresponds to typical loggings observed in HF field tests [START_REF] Hawkes | Interpretation of closure pressure in the unconventional montney using PTA techniques: the SPE Hydraulic Fracturing Technology Conference[END_REF]. The initial linear portion (Fig. 5) corresponds to the elastic behaviour of the intact rock. Once the fracturing threshold is reached, the injection pressure reduces due to the formation of cavities created by fractures. The injection pressure eventually reaches a steady state of fracture propagation pressure (∼5 MPa) slightly greater than the far-field minimum stress adopted in our simulations. Figure 6 shows the fracture trajectories for the two cases of isotropic and anisotropic stress states. These trajectories are plotted in the field of velocity magnitude as an indication of the elastic wave spreading in the medium due to fracturing. The comparison is shown for the same discrete time steps for each simulation. For this simple case of a homogeneous medium, initial fracture growth shows a complexity that deviates from the simplified model of bi-wing tensile fractures that are often used to represent hydraulic fractures Nordgren (1972). Similar observations of complex fracturing, as shown for the isotropic stress field have been observed in the Barnett shales in north Texas [START_REF] Vermylen | Hydraulic Fracturing, Microseismic Magnitudes, and Stress Evolution in the Barnett Shale[END_REF]. For the anisotropic stress field, fractures grow, as expected, in the direction of maximum stress, with curving fractures linked to mode II failure being more evident [START_REF] Zhou | Micromechanical Modelling of the Complete Stress-Strain Relationship for Crack Weakened Rock Subjected to Compressive Loading[END_REF].

Simulations of HF in jointed rock masses

Rock formations often contain pre-existing natural joints and fractures that act as zones of weakness and/or stress concentrators [START_REF] Gale | Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments[END_REF]. In some cases such joints may be created due to stresses caused by the drilling process [START_REF] Brudya | Drillinginduced tensile wall-fractures: Implications for determination of in-situ stress orientation and magnitude[END_REF], or through perforations in the well casing prior to fluid injection. In this section we first consider the effect of numerous space randomly distributed joints around the wellbore to investigate fracture interaction. We then consider the effect of lateral stresses from different joints close to the wellbore on HF and fracture trajectory as a key to understand HF induced asymmetric microseismicity Shapiro and Dinske (2009).

Heavily fractured rock formations

This section investigates fluid-driven fractures interaction with pre-existing joints in heavily fractured rock formations. It also explores possible micro-seismic events due to fracture shear slips incited by fracturing elastic waves. A set of space randomly distributed joints, (Fig. 7), with average length of 3.1 cm and fracture density of 320 (1/m 2 ) is created in the zone of intensive refinement show in Fig. 4. A free intact space is left around the borehole to allow for pure fluid-driven fracture growth before interaction with joints.

x-coordinate (m) (Longitudinal) The rock joints are implemented in the simulation before the application of the geostatic step (see section 6.2.1) and therefore lead to stress perturbations during the initial injection of fluid prior to fracture formation. This perturbation to initial (pre-fracturing) stress field affects the threshold of fracturing, typically increasing it by 0.25 MPa for the case of isotropic far-field stresses and by 0.45 MPa for the case of anisotropic far-field stresses. This will be further discussed in the Sect. (6.4.2).

Figure (8) shows the fracture patterns for the two cases of the far-field stresses at time 1.36 ms since injection. The joints around the borehole alternate the fluid-driven fracture patterns, compare with Fig. 6. Once intersected by a growing fracture, the joints obviously exploit fluid pressure and lead it to the tips where fracture initiation is the easiest. The joints also increase the extent to which fluid-driven fractures can reach for the same injected energy. Possible fracture shear slippage due to fracturing elastic waves is observed as well in both far-field states at some joints; some events are marked in the isotropic stress state. These slips are merely created by the fracturing elastic waves as they are not affected by fluid pressure. Such slips will lead to micro-seismic events.

Figure ( 9) demonstrates the micro-seismic events (MEs) locations at time steps of 1.28 and 1.4 ms for the case of isotropic far-field stress state. Each broken element represents a ME and the strain energy released is determined and used to calculate the event's magnitude Lisjak et al. (2013). The events' locations delineate the fracture patterns except for some scattered incidents (Fig. 9(a)) driven by the deformation of the medium (joint shear slippage). Figure (10) shows the magnitudes of the MEs which range between -4 and -5. Such values are not high and the events will not be felt on the surface, however they are in the range of induced HF microseismicity Schoenball et al. (2013). It is worthwhile to mention that joint shear slippage can also be incited by fluid diffusion altering critically stressed joints. However, this is to be implemented in our code and will be discussed in future works.

The space distribution of the pre-existing joints (Fig. 7) is generated using a Matlab Mfile, they are randomly distributed in space and their lengths are normally distributed around 3.1 cm. The results of the previous simulations can drastically change depending on the space distribution and lengths of these joints. Yet, this section helps to conclude the following remarks:

1. pre-existing joints create a pre-fracturing stress state that affects the patterns of the very initiated fractures, compare the isotropic states of Figs ( 6) and ( 8); 2. rock joints increase the extent to which fluid-driven fractures can reach; 3. once a joint is intersected by a fluid-driven fracture, it becomes easier for new fractures to initiate at the joint's tips; 4. MEs due to joint shear slippage can be created by rock deformation. Such events can go as far as fluid-driven fractures can reach as long as there are critically stressed joints that can be activated, Fig. 9(b).

Single isolated joints: Lateral stress gradient

In this section the effect of single isolated preexisting joints that do not intersect the wellbore is investigated in relation to fracture growth topology and fracturing pressure. We consider a single rock joint of length L (where L = 6D, D is the diameter of the wellbore) located at discrete distances (0.3L, 0. The presence of such joints around the wellbore as well as their orientations are expected to affect fracture growth. Indeed, HF induced micro-seismicity is related to fracture evolution and tends to concentrate in the regions where mode II (shear) fractures are prevalent [START_REF] Shapiro | Fluidinduced seismicity: Pressure diffusion and hydraulic fracturing[END_REF]. Fluid diffusion is also dominant in the regions of fracture evolution and therefore the long-term shearslippage micro-seismic/seismic events [START_REF] Shapiro | Fluidinduced seismicity: Pressure diffusion and hydraulic fracturing[END_REF].

The rock joints are implemented in the simulation before the application of the geostatic step (see section 6.2.1) and therefore lead to stress perturbations during the initial injection of fluid prior to fracture formation. All the simulations that will be presented in this section are performed under an anisotropic far-field stress state. Figure 12 shows contours of stress magnitude 4 for the reference case (no rock joints) and for the three cases of different joint orientation at a distance of 0.3L from the wellbore immediately prior to the onset of fracturing. It can be seen that for cases 1 and 2 appreciable perturbations in the stress field are clear when compared to the reference case, although for case 3 no appreciable perturbations are evident. These differences in stress perturbations lead to an increase in the onset of fracturing pressure/time step for when fracturing occurs. This is demonstrated by identical and earlier time step for the reference case and the case 3 of transversal joints. The existence of rock joints at this distance from the wellbore lead to stress concentrations around the rock joint, this is especially evident for case 2 where localized zones of shear failure/slip can been seen at the joint tips.

Figure 13 shows fracture trajectories for the three rock joint orientations described in Fig. 11 at a distance of 0.3L from the wellbore. It can be inferred that perturbations in the stress field resulting from the presence of the rock joints lead to significantly different fracture evolution/topology, compare case 1 at t = 1.26 ms in Fig. 13 with that of Fig. 6 at the same time step. For cases 1 and 2 with longitudinal and oblique rock joints, it is clear that fractures initiate in a direction that is oriented away from the zone of stress perturbation created by the preexisting joints. This conclusion is more evident in case 2 for an oblique joint. For case 3 with transversal joint, i.e. normal to the direction of maximum stress field, the perturbation to initial stress field does not appear to have a significant effect on the initiation or the resultant preferred direction of fracture growth. Once the fracture is initiated however, the stress field resulting from the rock joint seems to suppress any potential branching of the fracture.

The influence of the rock joint orientation can be also observed in the fracturing pressure profiles plotted in Fig. 14(a). For case 3, the threshold of fracturing pressure is not significantly altered and closely matches that observed for the reference case with no joints. However, the presence of longitudinal or oblique pre-existing joints, that leads to stress perturbations, has increased the fracturing threshold by ∼5.2%.

Figure 15 shows fracture topology for the three cases of the fracture joint at a distance 0.6L from wellbore. As expected, the effect of stress perturbation due to pre-existing joints becomes smaller as the joint is located further away from wellbore. Cases 1 and 2 still exhibit preferential growth of fractures towards left/asymmetric fracture growth, at least for early time frame of stimulation (t < 1.38 ms). However, fractures in case 3 exhibit the same fracture growth as if the rock joint did not exist (compare of instance case 3 at 1.26 ms in Fig. 15 with that in Fig. 6.

Fracturing pressure profiles of the three cases of different joint orientations at distance 0.6L from wellbore are shown in Fig. 14(b). The figure highlights that the breakdown pressure matches that for the case of no joints. In considering Fig. 14, it can be seen that the post-4 Stress magnitude is calculated based on the maximum principal stresses as σ 2 1 + σ 2 2 .

Anisotropic stress state

Case 1 (L) -0.3L Case 2 (I) -0.3L Stress magnitude (MPa) peak behaviour of the pressure is independent of the distance of the joint away from the wellbore.

For the rock joint at a distance L from the wellbore no appreciable differences in fracture trajectories can be qualitatively observed as shown in Fig. 16 when compared with those in Fig. 6 at 1.26 ms. The excessive branching on the right hand side of the wellbore for case 3 is assumed to be related to the relatively small time step that was used as the model was showing instability. This instability is related to the intensive refinement of the mesh in the direction if the joint, i.e. direction of maximum far-field stress.

FDEM-HF verification against field DFIT

The purpose of this section is to validate our approach of FDEM hydraulic fracturing by verifying the simulations against field data. Figure ( 12) dissects a comparison between the analytical response of the FDEM-HF and the field data of Encana DFIT at Montney formation. Clearly the numerical response captures the elastic response as well as the threshold of fracturing ∼76 MPa. The post-peak behaviour is showing a bit of discrepancy, i.e. the numerical response is showing higher fluid pressure values. This could be related to the fact that fluid diffusion as well as fluid leakoff are not implemented in our code so far. It is due to mention that only the elastic and the pre-shutdown behaviour is simulated in Fig. (12). The Shutdown and post shutdown behaviours require fully coupled poroelastic model which will be discussed in a future work.

Conclusions and future work

The FDEM method that combines the advantages of FEM and DEM is used in this study to simulate HF growth, with particular application to diagnostic fracture injection test (DFIT) analysis. The FDEM is capable of tracking fluiddriven fractures topology in homogeneous and pre-fractured media. Consistent with HF observations, the study shows that an isotropic stress state favours complex fracture growth compared with anisotropic stress state. Fluid-driven fracture patterns are considerably affected if the rock medium is heavily fractured. Pre-existing joints are deemed to increase the extent of fracture growth. Critically stressed joints may also slip leading to seismic events just because of formation deformation. In addition, the study shows that once a pre-existing joint is intersected by a fluid-driven fracture, fractures prefer to initiate at the joints' tips. Rock joints that are offset from the wellbore give rise to stress gradients that lead to asymmetric fracture growth with a tendency for fractures to initiate in a direction away from the rock joint. After intersecting a pre-existing joint, the fluid induced fracture system exploits and extends natural fractures.

The FDEM method reproduces features predicted by analytical models and empirical observations and thus is expected to provide a powerful tool to explore more complex fracturing scenarios. These capabilities are deemed important for simulating HF tests and predicting the direction and extent of created fractures for well monitoring real HF field tests. Future research will focus on simulating HF by FDEM is a fully coupled poroelastic frame work that can account for the effect of fluid diffusion and leakoff.

FFigure 1 :

 1 Figure 1: Interaction between computational modules of the Y-Geo code.

Figure 2 :

 2 Figure 2: (a) Constitutive behaviour of cohesive fracture elements. The curves represent a relation between normal and tangential bonding stress, σ and τ , versus crack relative displacement, δ (opening) and s (sliding). GIc and GIIc represent the amount of energy dissipated per unit length of fracture for mode I and mode II fractures respectively. (b) Elliptical coupling relationship between crack opening (δ) and crack slip (s) for mixed-mode fracturing (Eq. 6.3).

Figure 3 :

 3 Figure 3: Simulation of tensile fracture growth using Y-Geo. Only FEM elements adjacent to the fracture are shown, for simplicity. Once the strength limit of the material (ft) is reached the fracturing energy (GIc) starts dissipating and the fracture elements start yielding which is represented by the green area. Also note the change of tensile strength over the yielded elements.

Figure 4 :

 4 Figure4: The geometry and boundary conditions of BVP chosen for the HF simulations. The diameter of the vertical wellbore is 0.1 m. Two cases of far-field stress loading: 1-isotropic when the longitudinal stress is σ h ; 2-anisotropic when the longitudinal stress is σH .

Figure 5 :

 5 Figure 5: Numerical fracturing pressure profiles for the cases of isotropic and anisotropic far-field stress states. It is evident that HF is easier in media subjected to anisotropic stress fields.

Figure 6 :

 6 Figure 6: Fracture topology/trajectory for the cases of isotropic and anisotropic stress states and for two different time steps. The trajectories are plotted in the field of velocity magnitude.

Figure 7 :

 7 Figure7: A set of space randomly distributed joints is generated around the borehole to study fracture interaction and shear-slip micro-seismicity. Joint length is normally distributed around 3.1 cm and fracture density is 320/m 2 .

Figure 8 :

 8 Figure 8: Fracture patterns for the cases of isotropic and anisotropic stress states in a heavily fractured rock formation at time of 1.36 ms since injection. The trajectories are plotted in the field of velocity magnitude.

  Figure9: Micro-seismic events locations due to fracturing with an isotropic far-field stress state. The events delineate the fracture topology except for some scattered locations incited by fracturing elastic waves/formation deformation.

Figure 12 :

 12 Figure12: Effect of pre-existing joints on the field of stress magnitude just before the initiation of HF. The threshold of fracturing for case 3 is the same as the no joints case, i.e. t = 1.08 ms. This can be related to the reduced magnitude of perturbation on stress field created by the transversal joint.

Figure 13 :

 13 Figure 13: Fracture trajectories for the three cases of joint configurations described in Fig. 11 at time steps of 1.26, 1.38 and 1.68 ms and at distance 0.3L from wellbore. The trajectories are plotted in the field of velocity magnitude.

Figure 17 :

 17 Figure 17: Encana DFIT fluid pressure profile June/2015, Montney formation. The red dots represent the field data and the solid line exemplifies the numerical response by FDEM-HF.

Table 1 :

 1 Rock properties for a HF test of the BVP defined previously.

	Nature	Parameter	Value	Unit
	Elasticity	Drained Young's modulus, E	35	GPa
		Drained Poisson's ratio, ν	0.27	-
	Fracture	Tensile strength, f t	3	MPa
		Cohesion, c	24	MPa
		Mode I fracture energy, G Ic	5	N/m
		Mode II fracture energy, G IIc	80	N/m
		Material internal friction angle, φ i	38	( • )
	Computational Damping coefficient, µ	5.6 × 10 5 kg/m/s
		Normal contact penalty, p n	250	GPa.m
		Shear contact penalty, p t	30	GPa.m
		Fracture penalty, p f	120	GPa.m

Time since injection, t (ms)

Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid skinny triangles.

The presence of a mud cake creates a relatively impermeable barrier.

Reader is referred to Fjaer et al.[START_REF] Fjaer | Petroleum related rock mechanics[END_REF] for more details about deriving Eq. 6.8.

Acknowledgment

The authors would like to thank the sponsors of the Microseismic Industry Consortium and NSERC for their financial support.


David W Eaton
email: eatond@ucalgary.ca

Mason Mackay

Anton Biryukov

Murad Abuaisha

Identifying critically stressed faults without triggering slip : A discussion

Identifying and mapping critically stressed faults is a key first step for assessing hazards from induced seismicity. The classical paradigm of mapping faults with seismic or well log data may not be universally successful, however, as faults are sometimes cryptic and not easily identified. This is especially true for near-vertical strike-slip faults, for which vertical reflection offsets and/or diffractions may not be conspicuous, or in the case of low-angle thrust faults such as an intracutaneous wedge within a thick shale sequence. We propose alternative approaches for identifying critically stressed faults that involve: 1) modeling and mapping of stress perturbations caused by the presence of a critically stressed fault system, 2) mapping localized zones of anomalous pore pressure in a target shale formation, and 3) dynamic triggering of microseismicity due to high-amplitude surface waves resulting from large teleseismic earthquakes. The goal of this work is to promote further discussion.

Introduction

The identification and mapping of potential seismogenic faults without activating them is an important element for risk assessment for injection-induced seismicity [START_REF] Walters | Characterizing and Responding to Seismic Risk Associated with Earthquakes Potentially Triggered by Fluid Disposal and Hydraulic Fracturing[END_REF]. Traditionally, the identification of faults is carried out using well logs [START_REF] Rafiq | Reservoir Characterization using Microseismic Facies Analysis Integrated with Surface Seismic Attributes[END_REF] or by mapping reflection offsets in seismic profiles (Figure 1). However, this approach is not always viable, especially in the case of near-vertical strike-slip faults or sub-horizontal faults. Strike-slip faults commonly have a nearvertical fault orientation and may not result in any obvious vertical offset across the fault. Moreover, thrust-fault systems commonly contain low-angle (flat) segments. If such segments occur within a thick shale sequence they may be difficult to identify. In this paper, we propose alternative approaches to identification of critically stressed faults. These alternative approaches consist of:

1. Mapping perturbations in stress conditions near a fault, including changes in stress magnitude and direction (Figure 2). Such variations in stress could be identified using DFIT, wellbore breakouts [START_REF] Bell | Attempts to detect open fractures and nonsealing faults with dipmeter logs[END_REF] and cross-dipole sonic logs. Under favourable conditions it may be possible to map stress variations in detail using seismic data that is calibrated using these independent sources. 2. Mapping pore pressure compartments within a reservoir. We propose that faults may act as barriers to stress equalization within a reservoir. 3. Identifying critically stressed faults by mapping microearthquakes that are dynamically triggered by the passage of large-amplitude surface waves from distant earthquakes. Implementation of this approach would require deployment of a seismograph array in advance of hydraulic fracturing.

These three proposed methods are examined in the following sections.

Stress Modeling

In order to calculate how the stress field around a fault is perturbed both in magnitude and direction, the finite-difference code FLAC3D was used (Itasca, 2012). The finite-difference approach allows for materials to be represented as cubic polyhedral elements that are free to deform following Hooke's elastic constitutive law (Itasca, 2012),

∆σ ij = 2G∆ ij + K - 2 3 G ∆ kk δ ij , (35.1)
where Einstein notation is used implying the change in stress is a function of the Kroenecker delta δ ij , bulk modulus, K, and shear modulus G. An elastic material formulation without a failure criterion was used in order to represent the complete mechanics of stress distribution around a fault without the complexity of failure mechanisms that one would expect to be associated with dynamic processes.

The model approximates a stress condition that is characterized by both strike slip and thrust fault regimes based on the critically stressed crust theorem, wherein the ratio of maximum and minimum principal stresses, σ 1 to σ 3 is dependent on a coefficient of friction, µ, of 0.6 [START_REF] Zoback | Reservoir geomechanics[END_REF]:

σ 1 σ 3 = µ 2 + 1 + µ 2 2 . (35.2)
The convention used by the numerical code is that negative stresses are compressive. The stress state is at a modelled depth of 3 km assuming an average crustal density of 2700 kg/m 3 and a hydrostatic pore pressure gradient. The fault is approximated as a circle with a diameter of 3 km and the model size is 12 × 12 × 12 km . The 25 m computational grid is gradationally meshed in order to balance computation time and numerical accuracy with mesh elements at the fault face. The material properties and model parameters are summarized in Table 1. The effects of gravity were neglected and the model was initialized with an in situ stress state. Boundary conditions were applied with fixed corners and constrained to only allow for shear displacement at the model boundary. The model was then stepped to equilibrium and the fault region material was changed to a weaker material by assigning a fault material with half the bulk and shear modulus of the surrounding rockmass. The model was then stepped to a new equilibrium and analysis of the stress perturbation was made. The results of stress modelling are summarized in Figures 345.

Pore pressure compartments

In sedimentary basins, regions of pore overpressure can be compartmentalized by permeability barriers [START_REF] Hunt | Generation and Migration of Petroleum from Abnormally Pressured Fluid Compartments[END_REF]. Breaches in integrity of seals can occur due to fracturing, induced by generation of hydrocarbons within the compartments and/or the thermal expansion of pore fluids. An episodic process of resealing and breakout cycles in intervals of thousands of years has been proposed by [START_REF] Hunt | Generation and Migration of Petroleum from Abnormally Pressured Fluid Compartments[END_REF]. Figure 6 shows an example of compartmentalized pore overpressure within the Montney region of Alberta and British Columbia and its possible relationship to anomalous induced seismicity (Eaton et al., this volume). Table 1: FLAC3D fault zone modelling parameters. The following parameters are chosen to represent the rockmass and fault properties of an assumed fault zone that may be a candidate for induced seismicity.

Dynamic Triggering

Dynamic triggering happens as surface seismic waves from an initial earthquake propagate through the Earth's crust and perturb the stresses at local critically stressed faults, triggering secondary earthquakes. Once the seismic wave train has passed and ground shaking ends, the crust in the local area returns to its previous stress state modified according to the cumulative stress drops associated with locally triggered earthquakes. Dynamic triggering of secondary earthquakes was widely accepted in the scientific community following the 1992 M7.3 Landers earthquake in southern California. Earthquake rates increased dramatically across the western United States at distances well beyond the aftershock zone in a few days after the earthquake [START_REF] Hill | Seismicity remotely triggered by the magnitude 7.3 Landers, California, earthquake[END_REF]. Furthermore, dynamic triggering has been observed across the globe in various geologic and tectonic environments. It has been shown to take place at distances from the initial rupture varying from meters [START_REF] Kilb | Triggering of earthquake aftershocks by dynamic stresses[END_REF] to over thousands of kilometers [START_REF] West | Periodically Triggered Seismicity at Mount Wrangell, Alaska, after the Sumatra Earthquake[END_REF]. The main ruptures have been noticed to trigger other earthquakes at different time scales. In many cases secondary earthquakes occurred within minutes to hours following the primary seismic surface wave arrival [START_REF] West | Periodically Triggered Seismicity at Mount Wrangell, Alaska, after the Sumatra Earthquake[END_REF][START_REF] Gomberg | Earthquake nucleation by transient deformations caused by the M = 7.9 Denali, Alaska, earthquake[END_REF].

In other cases, earthquakes occurring weeks to months after the initial earthquake have been interpreted as a delayed response to dynamic triggering.

Recently, van der Elst et al. ( 2013) and [START_REF] Wang | Isolated cases of remote dynamic triggering in Canada detected using cataloged earthquakes combined with a matched-filter approach[END_REF] noticed that seismic activity on a fault in the areas of wastewater injection may be related to the stress perturbations caused by strong distant earthquakes. By using the novel match-filtering approach, they have managed to extract low-magnitude events from the continuous waveform data and establish the statistical causality between the major distant earthquakes and local increase in the seismic activity in the hydraulically stimulated areas.

Having been located and shown on the map, the series of triggered event locations may be spatially clustered, thus providing an estimate of the dimensions of the fault that has undergone dynamic triggering. This suggests that such remote triggering could potentially serve as a probe for outlining the areas where critically stressed faults are present, prior to any hydrocarbon field development. This procedure could serve as a preventive manner of accidental fault activation, causing unnecessary slip, and thus reduce the risk of a production disruption.

Conclusions and Outlook

Recognition and mapping of critically stressed faults is a key component of hazard assessment for anomalous induced seismicity. Under various common circumstances, fault imaging and mapping may be problematic using traditional methods, including 3-D seismic imaging or mapping with densely spaced well log data. We propose three novel methods for indirectly inference of the existence and location of critically stressed faults: mapping perturbations to regional stress, mapping compartments in pore pressure, and clustering of naturally occurring microearthquakes that are dynamically triggered by teleseismic surface waves. Ongoing research is focused on evaluation of the potential utility of these proposed methods.

Acknowledgments

Sponsors of the Microseismic Industry Consortium are thanked for their long-standing support of this research topic. We are particularly grateful to Canadian Discovery Ltd. for permission to show the pore-pressure data, as well as the Natural Sciences and Engineering Research Council of Canada (NSERC) and Chevron Canada for their support of an Industrial Research Chair in Microseismic System Dynamics at the University of Calgary. In addition, Aquaterra Water Management, ConocoPhillips Canada, Microseismic Canada and Nanometrics are thanked for sponsorship of a collaborative research program on multi-scale monitoring of hydraulic fracturing and wastewater injection, Duvernay region.

Figure 1 :

 1 Figure 1: Seismic line in the Fox Creek region showing interpreted faulting. From Green and Mountjoy, 2005.

Figure 2 :

 2 Figure 2: Left panel: deflection of maximum horizontal stress orientation in the neighbourhood of vertical open fractures that are not oriented perpendicular to the smallest principal stress. Right panel: deflection of mean borehole breakout azimuths from 29 wells in the Scotian Shelf, offshore eastern Canada. From[START_REF] Bell | Attempts to detect open fractures and nonsealing faults with dipmeter logs[END_REF]

Figure 5 :

 5 Figure 5: Maximum principal stress rotation. Points represent trends and plunge of the primary principal stress orientation. Deviation from a perfect North to south alignment indicates a rotation of the stress field by up to 25 degrees.

Figure 3 :

 3 Figure 3: Maximum principal stress magnitudes contours around a 3 km fault. Stress perturbations indicate a lowering of stress magnitudes in the core of the fault with increase stress magnitude lobes distally from the fault edges.

Figure 4 :

 4 Figure 4: Minimum principal stress magnitude contours. Increased compressive stresses (orange) are observed distally to the fault, while zones of relaxation (green) occur in other regions.

Figure 6 :

 6 Figure 6: Heterogeneous (compartmentalized) pore-pressure gradient in the Montney trend, courtesy Canadian Discovery, Ltd. Black dots show M ≥ 2.5 earthquakes since 2008, suggesting that pore-pressure compartments may be linked to faults.

Figure 7 :

 7 Figure 7: Seismicity rate increases associated with Hector Mine and Landers mainshocks, showing that dynamically triggered events tend to align along faults. From Gomberg et al. (2001).