**3. Methodology**

Weimer [1], Ohtake [2], Fletcher and Sykes [3], Pearson, [4], Talwani and Acree, [5], Simpson

The occurrence of induced seismicity, whether through hydraulic fracturing, enhanced geothermal systems, or carbon dioxide (CO2) sequestration, represents an important consid‐ eration 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

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 buildup 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 corre‐ sponding maximum seismic magnitude that is likely to occur in response to the hydraulic

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

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],

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.

energy already stored in the surrounding rocks to be released (Healy et al. [9]).

**2. Effect of fluid injection on induced seismicity**

shut-in, referred to here as post-injection seismicity.

and in the Los Angeles basin by Teng et al. [14].

et al. [6], Zoback and Harjes [7]).

478 Effective and Sustainable Hydraulic Fracturing

fracture treatment.

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 impor‐ tant 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:

$$q = ka^3 \frac{\Delta P}{l} \tag{1}$$

**Discontinuity property Incipient fractures Fault Units** Residual tensile strength 0.0 0.0 degrees Normal stiffness 1e4 10 MPa Shear stiffness 1e3 1 MPa

A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing

http://dx.doi.org/10.5772/56191

481

**Table 1.** Properties assigned to the modeled planes of weakness and fault.

**Figure 1.** Rock mass model with two sets of weakness planes and a fault.

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.

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

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

**5. Post-injection seismicity simulations**

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 shutin. 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/m3 , 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.



**Table 1.** Properties assigned to the modeled planes of weakness and fault.

<sup>3</sup> *<sup>P</sup> q ka <sup>l</sup>*

contact is broken, the fluid flows into the joint.

480 Effective and Sustainable Hydraulic Fracturing

contacts.

for shale (density=2500 kg/m3

**4. Simulation setup and input parameters**

stresses were assumed to be 30 and 20 MPa, respectively.

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

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 shutin. 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

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

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

**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

, Young's modulus=30 GPa, Poisson's ratio=0.25). The incipient

<sup>D</sup> <sup>=</sup> (1)

**Figure 1.** Rock mass model with two sets of weakness planes and a fault.

#### **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).

**Figure 2.** Pore pressure distribution at the time of shut-in.

#### **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. Displace‐ ments 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.

**Figure 4.** Fault shear displacement 2 ½ hours after shut-in.

**Figure 3.** Pore pressure distribution 2 ½ hours after shut-in.

A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing

http://dx.doi.org/10.5772/56191

483

A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing http://dx.doi.org/10.5772/56191 483

**Figure 3.** Pore pressure distribution 2 ½ hours after shut-in.

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

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

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. Displace‐ ments 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.

treatment was applied for a period of 50 minutes.

model simulations performed here).

482 Effective and Sustainable Hydraulic Fracturing

**Figure 2.** Pore pressure distribution at the time of shut-in.

**5.2. Shear slip of critically stressed fault**

**Figure 4.** Fault shear displacement 2 ½ hours after shut-in.

#### **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]):

$$M\_o = \mu A \overline{D} \tag{2}$$

Seismic moment can be converted into a moment magnitude using the following relationship

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

Figure 6 uses a well-established seismological relationship that correlates earthquake magni‐ tude 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

**Figure 6.** Relationship among various scaling parameters for earthquakes. The larger the earthquake, the larger the fault and amount of slip, depending on the stress drop in a particular earthquake. Observational data indicate that

earthquake stress drops range between 0.1 and 10 MPa (modified after Zoback and Gorelick, [8]).

local magnitude for the calculated moment magnitude is equal to 2.45.

2 / 3log 10.7 *M M w o* = - (3)

A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing

http://dx.doi.org/10.5772/56191

485

(Hanks and Kanomori [21]):

Mw is moment magnitude (dyne.cm), and

Mo is seismic moment (dyne.cm).

Where:

Where:

Mo is the seismic moment (dyne∙cm),

μ is the rupture strength of the rock along the fault (dyne/cm2 ),

A is the rupture area (cm2 ), and

*<sup>D</sup>*¯ is the average amount of slip (cm).

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

**Figure 5.** Saturation of various magnitude scales: Mw (moment magnitude), ML (Richter local magnitude), Ms (surface wave magnitude), mb (short-period body wave magnitude), mB (long-period body wave magnitude), and MJMA (Japa‐ nese Metrological Agency magnitude). After Idriss [22] and Kramer [20].

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

$$M\_w = 2 \, / \Im \log M\_o - 10.7 \,\tag{3}$$

Where:

**6. Seismic moment and moment magnitude calculations**

μ is the rupture strength of the rock along the fault (dyne/cm2

Here, the rupture strength of the fault is equal to 1e10 dyne/cm2

), and

resulting seismic moment is equal to 5.4e19 dyne.cm.

nese Metrological Agency magnitude). After Idriss [22] and Kramer [20].

Where:

Mo is the seismic moment (dyne∙cm),

484 Effective and Sustainable Hydraulic Fracturing

*<sup>D</sup>*¯ is the average amount of slip (cm).

A is the rupture area (cm2

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]):

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 cm2

**Figure 5.** Saturation of various magnitude scales: Mw (moment magnitude), ML (Richter local magnitude), Ms (surface wave magnitude), mb (short-period body wave magnitude), mB (long-period body wave magnitude), and MJMA (Japa‐

(2)

, the length of the fault that

), the

),

*M AD <sup>o</sup>* = m

Mw 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 magni‐ tude 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

**Figure 6.** Relationship among various scaling parameters for earthquakes. The larger the earthquake, the larger the fault and amount of slip, depending on the stress drop in a particular earthquake. Observational data indicate that earthquake stress drops range between 0.1 and 10 MPa (modified after Zoback and Gorelick, [8]).

event, or stress drop. For the case modeled here, the 0.8 mm of average slip produced corre‐ sponds to an event with a magnitude equal to 2.45, which is in the range shown in Figure 6.

[2] Ohtake, M. Seismic Activity Induced by Water Injection at Matsushiro, Japan. Jour‐

A Numerical Investigation of Fault Slip Triggered by Hydraulic Fracturing

http://dx.doi.org/10.5772/56191

487

[3] Fletcher, J. B, & Sykes, L. R. Earthquakes Related to Hydraulic Mining and Natural Seismic Activity in Western New York State. Journal of Geophysical Research

[4] Pearson, C. The Relationship Between Microseismicty and High Pore Pressures Dur‐ ing Hydraulic Stimulation Experiments in Low Permeability Granitic Rocks. Journal

[5] Talwani, P, & Acree, S. Pore Pressure Diffusion and the Mechanism of Reservoir-in‐

[6] Simpson, D. W, Leith, W. S, & Scholz, C. H. Two Types of Reservoir Induced Seis‐ micity, Bulletin of the Seismological Society of America (1988). , 78, 2025-2040.

[7] Zoeback, M, & Harjes, H. P. Injection Induced Earthquakes and the Crustal Stress at 9 km Depth at the KTB Deep Drilling Site, Germany. Journal of Geophysical Re‐

[8] Zoback, M. D, & Gorelick, S. M. Earthquake Triggering and Large-scale Geologic Storage of Carbon Dioxide. Proceeding of the National Academy of Science of the

[9] Healy, J. H, Rubey, W. W, Griggs, D. T, & Raleigh, C. B. The Denver Earthquakes,

[10] Evans, D. M. The Denver Area Earthquakes and the Rocky Mountain Arsenal Dis‐

[11] Pakiser, L. C, Eaton, J. P, Healy, J. H, & Releigh, C. B. Earthquake Prediction and

[12] Healy, J. H, Lee, W. H. K, Pakiser, L. C, Raleigh, C. B, & Wood, M. D. Prospects for

[13] Gibbs, J. F, Healy, J. H, Raleigh, C. B, & Coakley, J. Seismicity in the Rangely, Colora‐ do Area: 1962-1970. Bulletin of Seismological Society of America (1973). , 63,

[14] Teng, T. L, Real, C. R, & Henyey, T. L. Microearthquakes and Water Flooding in Los Angeles, Bulletin of Seismological Society of America (1973). , 63, 859-875.

[15] Charlety, J, Cuenot, N, Dorbath, L, Dorbath, C, Haessler, H, & Frogneux, M. Large Earthquakes During Hydraulic Stimulations at the Geothermal Site of Soultz-sous-Forets. International Journal of Rock Mechanics and Mining Sciences (2007). , 44,

Earthquake Prediction and Control, Techtonophys (1972). , 319-332.

United States of America (PNAS) (2012). , 109(26), 10164-10168.

nal of Physics of the Earth (1974). , 22, 163-176.

of Geophysical Research (1981). , 86, 7855-7864.

duced Seismicity, Pageoph (1985). , 122, 947-965.

posal Well, Mountain Geologist (1966). , 3, 23-36.

Control, Science (1969). , 166, 1467-1474.

(1977). , 82(26), 13767-3780.

search (1997). B8): 18477-18491.

Science (1968). , 161, 1301-1310.

1557-1570.

1091-1105.
