Properties assigned to the modeled planes of weakness and fault.

The study of fault slip in response to fluid injection offers a means to understand the complex hydromechanical behavior of shale gas and oil reservoir systems during hydraulic fracturing operations, together with the induced seismicity, and corresponding mitigation measures, arising from such events. In this paper, a series of numerical simulations are performed to investigate the relationship between hydraulic fracturing (i.e. fluid injection) and the response of a naturally fractured rock mass to transient fluid pressures. The analysis is carried out using the discontinuum-based distinct-element program UDEC assuming a fracture flow system. The conceptual reservoir model consists of a critically stressed fault plane and the surrounding rock mass containing planes of weakness, for which a hydraulic fracture is numerically simulated and the response modeled using a transient, coupled hydro-mechanical solution. The results demonstrate the influence of fluid diffusion generated by the fracing fluid after shut-in on the triggering of fault slip. The simulation is then used to interpret the associated seismic events and their relationship to the injections and shut-in pressures, and to estimate the maximum magnitude of the induced seismic event.

## 1. Introduction

Hydraulic fracture operations involve the injection of fluids from a wellbore into a formation to maximize the extraction of oil and gas resources from the reservoir. These stimulations allow the development of previously uneconomical reserves, for example shales gas. However, it has been well established that these fluid injections also induce seismicity (e.g., Hollister and Weimer [1], Ohtake [2], Fletcher and Sykes [3], Pearson, [4], Talwani and Acree, [5], Simpson et al. [6], Zoback and Harjes [7]).

The occurrence of induced seismicity, whether through hydraulic fracturing, enhanced geothermal systems, or carbon dioxide (CO_{2}) sequestration, represents an important consideration in the risk profile of an injection well. Zoback [8] suggested that because of the critically stressed nature of the crust, fluid injection in deep wells can trigger earthquakes through increases in pore pressure (and decreases in effective stress) in the vicinity of the pre-existing faults. The increased pore pressure reduces the shear resistance to fault slip, allowing elastic energy already stored in the surrounding rocks to be released (Healy et al. [9]).

This paper describes a numerical experiment investigating the influence of pore pressure diffusion on the modeling of seismicity after a hydraulic fracture treatment. One objective is to evaluate the utility of 2-D distinct element techniques to model both pore pressure build-up and diffusion along existing rock weakness planes in response to fluid injection and subsequent changes to the effective stress field. A second objective is to estimate the corresponding maximum seismic magnitude that is likely to occur in response to the hydraulic fracture treatment.

## 2. Effect of fluid injection on induced seismicity

The occurrence of injection-induced seismicity is usually confined in both space and time, with pressure build-up and diffusion controlling the spatial and temporal pattern of the seismicity after a hydrofrac treatment. Spatially, the problem can be conceptualized as a scenario in which a wellbore in the vicinity of a critically-stressed fault is pressurized and the injected fluids diffuse towards the fault, increasing pore pressures and reducing the effective normal stress until slip is triggered along a portion of the fault. The transient nature of the fluid front as it radiates outward from the injection well allows for the triggering of events after injection and shut-in, referred to here as post-injection seismicity.

The importance of fluid pressure in generating induced seismicity was demonstrated in 1962 in Denver, Colorado, when it was observed that a series of seismic events were being generated shortly after a chemical weapons plant had begun injecting contaminated wastewater down a deep well. Detailed studies by Hollister and Weimer [1] and Healy et al. [9] confirmed the causal observations of Evans [10] that the seismic events were induced by the waste fluid injection. Similar relations between local seismicity and fluid injection were likewise observed at the Rangely Oil Field, Colorado by Pakiser et al. [11], Healy et al. [12] and Gibbs et al. [13], and in the Los Angeles basin by Teng et al. [14].

More recently, post-injection seismicity associated with gas shale and geothermal projects have entered the public spotlight with heightened sensitivity directed towards hydraulic fracturing practices and induced seismicity. Between 2000 and 2005, numerous seismic events were recorded at the European geothermal project in Soultz-sous-Forêts, France, with events as large as M=2.5 being recorded during injection and M=2.6 several days after shut-in (Charlety et al. [15]). For stimulations carried out in 2000, 2003 and 2004 the largest seismic events occurred just after or after shut-in. Similar experiences were encountered at Basel, Switzerland where induced seismicity reached a Richter Magnitude of 3.4 despite precautionary reductions of the injection rate, leading to suspension of its hot dry rock enhanced geothermal systems project (Haring et al. [16]).

Experiences with hydraulic fracturing treatments in gas shales have generated similar seismicity, generally viewed as low risk due to the remote nature of the locations. In the Bowland shale, Lancashire, UK, two notable events were recorded with Richter Magnitudes of 1.5 and 2.3. These occurred almost 10 hours after shut-in (de Pater and Baisch [17]). de Pater and Baisch [17] suggested that reducing the treatment volume is one means to mitigate induced seismicity. In a detailed study in the Horn River Basin in northeastern British Columbia, Canada, approximately 38 events ranging between Richter Magnitudes 2.2 and 3.8 were recorded between 2009 and 201l, with only one being felt on surface in this remote region (BCOG, [18]). In their 2012 report of the findings from this investigation, the B.C. Oil and Gas Commission found that several of the Horn River events greater than M 2.0 were located along faults intersecting the wellbores. They also note that there were other instances where faults intersected wellbores without anomalous events being detected (BCOG, [18]).

## 3. Methodology

The Distinct Element Method (DEM) applies a Lagrangian formulation to compute the motion and interaction between a series of discrete deformable blocks, representing the problem domain, via compliant contacts and Newton’s equation of motion (Cundall and Hart, [19]). This enables the problem domain to be divided through by one or more discontinuity sets of variable orientation, spacing and persistence. One fundamental advantage of the DEM is that pre-existing (natural) joints in the rock mass can be modeled explicitly and allow for joints to undergo large deformations in shear (slip) or opening (dilation). The 2-D commercial code UDEC (Universal Distinct Element Code; Itasca Consulting Group, 1999) is used here to model the response of a jointed rock mass subjected to static loading and hydraulic injection.

UDEC is capable of modeling the behavior of weak jointed rock masses in which both the deformation and yielding of weak rock and slip along pre-existing discontinuities are important controlling factors. Progressive failure associated with crack propagation and fault slip can be simulated by the breaking of pre-existing contacts between the pre-defined joint bounded blocks, which although deformable, remain intact.

Key for simulating hydrofracturing, UDEC has the capability to model fluid flow through the defined fracture network. A fully coupled hydro-mechanical analysis can be performed, in which fracture conductivity is dependent on mechanical deformation of joint apertures and, conversely, joint water pressures can affect the mechanical computations of joint aperture. The blocks in this assemblage are treated as being impermeable, and fracture flow is calculated using a cubic law relationship for joint aperture:

where, *k* is a joint conductivity factor (dependent on the fluid dynamic viscosity), *a* is the contact hydraulic aperture, Δ*P* is the pressure difference between the two adjacent domains, and *l* is the length assigned to the contact between the domains. Since the UDEC formulation is restricted to the modeling of fracture flow, it should be noted that leak-off along the fractures diffusing into the rock matrix is assumed to be negligible (only leak-off into other fractures is considered). Furthermore, the cubic law flow assumption limits tortuosity. When a joint contact is broken, the fluid flows into the joint.

## 4. Simulation setup and input parameters

The rock mass modeled in this study is represented by two persistent orthogonal planes of weakness (Figure 1). These serve as incipient planes along which hydrofrac propagation is restricted. The simulation of induced seismicity is executed through the inclusion of a fault, which extends across the model. An injection well is located such that a significant portion of the fluid injected diffuses towards the fault and eventually penetrates the fault following shut-in. The fluid pressure decays slowly after the injection and the disturbed pressure front diffuses through the surrounding rock mass. Although the fluid pressures decrease with time and distance, there is still sufficient pressure to trigger fault slip. The fault model is based on interactions between neighboring fault segments allowing the model to simulate slippage on a single contact together with the subsequent interactions and responses of its neighboring contacts.

The input material parameters include both those for the rock matrix and incipient planes of weakness and fault. The rock matrix was modeled as being elastic, assuming typical values for shale (density=2500 kg/m^{3}, Young’s modulus=30 GPa, Poisson’s ratio=0.25). The incipient planes of weakness and fault were modeled assuming a Coulomb-slip constitutive model with both peak and post-peak properties. These are given in Table 1. The depth of the injection and horizontal plane represented by the model is 1000 m. The maximum and minimum horizontal stresses were assumed to be 30 and 20 MPa, respectively.

Discontinuity property | Incipient fractures | Fault | Units |

Friction angle | 30 | 20 | degrees |

Residual friction angle | 25 | 20 | MPa |

Cohesion | 1.0 | 0.0 | MPa |

Residual cohesion | 0.0 | 0.0 | MPa/m |

Tensile strength | 0.5 | 0.0 | MPa/m |

Residual tensile strength | 0.0 | 0.0 | degrees |

Normal stiffness | 1e4 | 10 | MPa |

Shear stiffness | 1e3 | 1 | MPa |

## 5. Post-injection seismicity simulations

Numerical simulations were performed using the model developed to investigate the influence of hydraulic fracturing on fluid pressure changes around the neighboring fault and any subsequent shear slippage along the fault. Both fluid pressure build-up during the injection until the time of shut-in as well as fluid pressure diffusion after the shut-in were considered.

### 5.1. Fluid pressure build-up during injection and diffusion after shut-in

Experience gained from mapping hundreds of hydraulic fracturing treatments with downhole geophones has shown that the occurrence of seismic event induced during a treatment is greatly influenced by the injection volume and rate used. Here, the hydraulic fracturing simulation was conducted by pressurizing the wellbore in the vicinity of a critically stressed fault (Figure 1). Figure 2 shows the joint fluid pressure distribution in the rock mass at the time that fluid injection is stopped (i.e., shut-in). It can be seen that the fluid pressure has not reached the fault at the end of the hydraulic fracture treatment and shut-in. Here the fluid pressure treatment was applied for a period of 50 minutes.

After shut-in, the injected fluids continue to diffuse and radiate towards the fault, eventually penetrating it (Figure 3). The response of the fault to the elevated fluid pressures then depends on the spatial and temporal characteristics of the diffusion pulse, with a series of slip events being produced as opposed to a single event. The strongest slip event was typically observed late in the sequence, sometimes long after injection had stopped (up to 150 minutes for the model simulations performed here).

### 5.2. Shear slip of critically stressed fault

The shear slip distribution along the fault, after 50 minutes of injection and 150 additional minutes of shut-in, is presented in Figures 4. The figure shows a distribution of slip magnitudes along the length of the fault, with a maximum fault slip of approximately 10 mm. Displacements between 5 and 10 mm were observed along 680 meters of the total 800 meter fault length, with an average fault slip of 8 mm. This average slip magnitude was then used to calculate the maximum earthquake magnitude, as presented in the following section.

## 6. Seismic moment and moment magnitude calculations

The seismic moment is a direct measurement of the energy released by a seismic event and is related to the strength characteristics of the fault. It can be calculated as follows (Kramer, [20]):

Where:

Mo is the seismic moment (dyne∙cm),

µ is the rupture strength of the rock along the fault (dyne/cm^{2}),

A is the rupture area (cm^{2}), and

Here, the rupture strength of the fault is equal to 1e10 dyne/cm^{2}, the length of the fault that slipped is equal to 6.8e4 cm (680 m), and the average amount of slip is equal to 0.8 cm. Assuming a unit depth of fault slip due to the 2-D nature of the model (i.e., rupture area = 6.8e9 cm^{2}), the resulting seismic moment is equal to 5.4e19 dyne.cm.

Seismic moment can be converted into a moment magnitude using the following relationship (Hanks and Kanomori [21]):

Where:

M_{w} is moment magnitude (dyne.cm), and

Mo is seismic moment (dyne.cm).

The moment magnitude calculated from equation 3 is 2.45. Figure 5 is used to convert the calculated moment magnitude to the Richter local magnitude. As seen in Figure 5, for earthquake magnitudes smaller than 6, the moment and local magnitudes are near equal; the local magnitude for the calculated moment magnitude is equal to 2.45.

Figure 6 uses a well-established seismological relationship that correlates earthquake magnitude to the size of the fault that slipped and seismic moment (Stein and Wysession [23]). This suggests that only faults that are at least tens of kilometers long are capable of producing large seismic events with magnitudes exceeding 6 (Zoback and Gorelick [8]). Zoback and Gorelick [8] note here that the fault size in Figure 6 is a lower bound value that refers to the size of the fault segment that slips in a given earthquake. As a geological feature, the total fault length is generally much longer than the part (segment) that slips during an individual event. As shown in Figure 6, an active fault slip segment of 680 meters, as produced in the UDEC model, is capable of producing an earthquake with a magnitude between 2.3 and 3.8 for a fault slip displacement between 1mm and 1 cm, depending on the magnitude of stress released by the event, or stress drop. For the case modeled here, the 0.8 mm of average slip produced corresponds to an event with a magnitude equal to 2.45, which is in the range shown in Figure 6.

## 7. Discussion and conclusions

A numerical modeling study has been carried out to investigate the application of the distinct element technique to the simulation of fault slip and induced seismicity resulting from a hydraulic fracture treatment via a wellbore in the vicinity of a natural fault. The model was able to predict the occurrence of post-injection seismicity in response to diffusion of the injected fluids, in a system governed by fracture permeability, long after the hydraulic fracturing treatment is finished and shut-in is initiated.

In most areas, regional-scale faults should be easily identified during geological site characterization studies. Smaller-scale faults or those that are shallow dipping (i.e., that do not daylight on surface and therefore would not be detected through surface mapping) may be more difficult to locate *a priori*. Experience gained from monitoring hundreds of hydraulic fracturing treatments with downhole geophones has shown that the occurrence of seismic events induced by the treatments is greatly influenced by the injection volume used during the operation. Improved understanding of this condition will allow designers and operators to control the amount of injection during a hydraulic fracturing treatment to minimize fluid pressure diffusion and subsequent slip of a neighboring fault.

The simulations presented here show great potential in providing a deeper understanding of the effect of natural fracture systems and pressure diffusion of fracing fluids into a neighboring fault causing shear slip and induced seismicity after a hydraulic fracture treatment. Future work will include modeling of hydrofrac fluid diffusion into a fault in the vicinity of the hydrofrac treatment site in a three-dimensional model.