Shock and Vibration
Volume 2021 (2021), Article ID 9997597, 21 pages
https://doi.org/10.1155/2021/9997597

Research Article

A Numerical Investigation of the Hazardous Injection Area of Induced Earthquake during Hydraulic Fracturing

Zhiyong Niu,¹ Shiquan Wang,¹ Hongrui Ma,¹ Songbao Feng,¹ Hengjie Luan,² and Yibin Yang³

¹School of Resources and Civil Engineering, Suzhou University, Suzhou 234000, China
²College of Energy and Mining Engineering, Shandong University of Science and Technology, Qingdao 266590, China
³Wuhu Natural Resources and Planning Bureau, Wanzhi District Branch, Wuhu 241100, China

Correspondence should be addressed to Shiquan Wang; wsq@ahszu.edu.cn

Received 4 April 2021; Accepted 20 May 2021; Published 28 May 2021

Academic Editor: Feng Du

Copyright © 2021 Zhiyong Niu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Hot dry rock (HDR) geothermal energy has many advantages, such as being renewable, clean, widely distributed, and without time and weather limitations. Hydraulic fracturing is usually needed for the exploitation of HDR geothermal energy. It has many hidden faults in the reservoir/caprock sequences. Injecting fluid into underground formations during hydraulic fracturing often induces fault slip and leads to earthquakes. Therefore, to well understand the induced fault slip and earthquakes is important for the applications and development of HDR geothermal exploitation. In this study, we investigated the hazardous injection area of the induced earthquakes during hydraulic fracturing. The study was based on a hydraulic fracturing test in Qiabuqia geothermal field in China. According to the field, a fault-surrounding rock-fracturing region system was developed to study the influences of fluid injection on the stability of the specific fault. A total of 60 hydraulic fracturing regions and 180 numerical experiments were designed. The results revealed that the hazardous injection regions that threaten the fault’s stability were near to the fault and concentrated on the following four areas: (a) above the top of the fault in underlying strata; (b) above the top of the hanging wall of the fault in underlying strata; (c) near to the fault planes in both footwall and hanging wall; (d) at the bottom of the footwall of the fault in underlying strata. The hazardous injection area can be controlled effectively by adjusting the injection pressure.

1. Introduction

Geothermal energy has great potential to simultaneously solve the energy crisis and climate problems due to its significant advantages. Firstly, the geothermal resource is rich. According to the estimation of the World Energy Council, it has about of geothermal energy (more than 10,000 times the available energy from all total fossil fuels worldwide), which can meet the energy demand for thousands of years according to the present global energy demand [1]. Secondly, it is renewable. It is mainly generated from the radioactive decay of materials within the Earth [2]. Thirdly, it is clean and green. The generation of geothermal energy produces no pollutants. Fourthly, it has no time and weather limitation. It can supply energy throughout the year due to the constant flow of heat from the Earth [3]. In addition, geothermal is widely distributed. Therefore, full exploitation of geothermal resources plays an important role in energy supply and environmental protection for all countries in the world.

According to the temperature and storage depth, geothermal energy can be divided into three categories: shallow geothermal energy, underground thermal water resources, and HDR resources [4]. Among them, HDR geothermal energy has the greatest potential. It is abundant and available almost in any area up to a certain depth in the world. It has about  JHD energy just within an accessible zone less than 10 km, becoming a strategic choice to solve the energy crisis and reduce environmental pollution [5]. Exploiting geothermal energy usually needs heat carrier fluids to bring the underground heat to the surface. Cold water can be injected into the hot dry rocks and then be heated. Heat energy can be then obtained through the heated water. For developing HDR energy economically on a commercial scale, the adequate fluid always needs to be supplied in the targeted HDR masses. However, the natural HDR mass (as shown in Figure 1(a)) is usually compact and impermeable. There is little or no water or steam in it. In addition, due to its low permeability, it does not allow fluid injecting or flowing in it. In order to solve these problems, an Enhanced Geothermal System (EGS) is usually needed [6]. In such systems, the permeability of rocks of the target reservoir is enhanced by using artificial methods (usually hydraulic fracturing). Cold fluid then can be injected into the reservoir and be heated by the hot rocks (Figure 1(b)).The hot water can be extracted to the surface to release the captured heat for power generation or other uses (Figure 1(c)). In the course of the fluid being injected into underground formations, seismicity is frequently induced. This phenomenon has been reported in many countries and regions, such as Australia [7], Germany [8], America [9], France [10], South Korea [11], Philippines [12], Switzerland [13], New Zealand [14], and Japan [15]. Large-scale, big, destructive earthquakes have exerted a great influence on the geothermal industry worldwide [16]. Induced seismicity by industrial activities often raises public concern and worry. Massive fluid injection during hydraulic fracturing and the subsequent heat mining also often induce seismicity, which may delay or hinder the acceptance of EGS. A better understanding of fluid-injection-induced fault reactivation is necessary for improving the public acceptance of EGS.

(a)

(b)

(c)

Figure 1: The EGS of HDR geothermal exploitation. (a) Natural HDR. (b) Hydrofracturing. (c) Heat exploitation.

The massive fluid injected into underground formations can often induce fault reactivation and slip. Several studies have explored this effect, and several different mechanisms have been proposed to explain it [17–21].The pore-pressure increase is considered to be the main mechanism to explain the occurrence of induced seismicity in underground fluid-injection projects. Injecting fluid into reservoir decreases the effective stress of faults which will cause them to to slip. [22]. However, the induced fault slip can also occur when the fault does not exist within the subsurface targeted zone (no longer the problem of effective stress) [23]. When the fluid is injected, the injection zone is expanded which may cause stress perturbation. The stress perturbation can be transmitted to the surrounding faults and leads them to slip. In Oklahoma, USA, the injection-linked seismicity occurs up to 35 km away from the injection wells [24]. It revealed that induced earthquakes may extend far away from the hydraulically connected region. In summary, both pore-pressure increase (direct pressure effect) and stress perturbation (indirect pressure effect) due to fluid injection can cause fault slip and induce earthquakes. Therefore, a full assessment of the stability of the fault beyond the targeted injection zones is also very important.

Numerical simulators have become important tools to solve the problems that arise in fluid-injection engineering in the deep subsurface [22, 25, 26]. In this paper, the influence of fluid injection during hydraulic fracturing on the fault’s stability was fully investigated by ABAQUS. The study took the geothermal potential area, i.e., Qiabuqia Geothermal Field (located in the east of Gonghe Basin in Northwest China) as a case. The numerical model was based on the field measurements (e.g. length, orientation, and inclination of faults) and the experimental data in field test and laboratory test (Density, Young’s modulus, Poisson’s ratio, and Porosity, etc.). In this paper, we mainly investigated the indirect pressure effect on the stability of the fault. Therefore, we only explored the stability of faults that exist not within the targeted injection zone (as shown in Figure 2). Our research focused on the spatial extent of induced seismicity as the following three aspects. (1) The influence of the spatial distribution of injection zone on induced fault slip and magnitude of seismicity. For a practical EGS project, faults naturally exist, but the injection zones are completely determined by human-made factors. For a given fault in the site, it is necessary to find out which injection zone is more likely to lead the fault to slip and induce seismicity. (2) The influence of the change of distance between the injection zone and fault on fault slip. (3) The dangerous injection region significantly affects the fault. (4) How to reduce the EGS-induced fault slip and the magnitude of seismicity.

Figure 2: Schematic of injection zone and its surrounding faults.

2. Theoretical Background

2.1. Linear Elastic Material Model

In the Qiabuqia Geothermal Field, the HDR formations mainly consist of hard granite. Granite is an elastic material with high elastic modulus and low permeability. The linear elastic material model can be used to describe its mechanical properties and deformation characteristics. Therefore, the simulation of the granite reservoir and fault was based on the linear theory of elasticity. The total stress is defined from the total elastic strain as [27]:

, where is the total stress; is the fourth-order elasticity tensor; and is the total elastic strain.

2.2. Porous Elastic Material Model

Hydraulic fracturing is usually used for establishing EGS in hot dry rocks. Artificial fracture networks can be created in the targeted region with high porosity and permeability. The porous elastic material model can be used to describe the fracturing region as shown in the flowing equations [28, 29]:

, where is the shear modulus; is the strain; is the stress; and are the pore-pressure and first stress invariant from the prefaulting state, respectively; and are the undrained and drained Poisson’s ratios, respectively; is the Kronecker delta; is the change in fluid mass per unit volume, kg/m³; is the density of the pore fluid in the reference state, kg/m³; is Skempton’s coefficient; and are drained bulk modulus and solid grain’s bulk modulus; and are the unjacketed pore compressibility and pore fluid bulk modulus, respectively; and is porosity.

The hydromechanical coupling is accomplished through the governing equation of mass conservation for the infiltrating pore fluid:

, where is the divergence of the water flux; is the direction of fluid mass flow rate; is fluid mass; and is time.

2.3. Contact Formulations for Fault Modeling

The fault and its surrounding rock are two deformable bodies. Finite-sliding formulation can be used to describe the interaction between them. Master-slave contact method was used to accomplish the interaction analysis of the fault and its surrounding rock. In our study, the surface of the fault-surrounding rock was the master surface and the surface of the fault was the slave surface. The slip of the slave surface along the master surface can be solved by the following equation:

, where is the slip; is the coordinate of the point on the contact segment; and are the interpolation functions of the point on the segment; and are the tangent and normal, respectively; and is the curvature of the contact segment. For the detailed derivation of this equation, refer to ABAQUS theory manual [27].

2.4. Moment Magnitude Scaling

The energy released by fault slip determines the magnitude of the earthquake. Seismic moment M₀ proposed by Aki was used to quantify the released energy, as follows [30]:

, where μ is rigidity of the fault, Pa; A is fault rupture area, m²; and d is the average slip of the rupture area, m.

Then, can be used to calculate the moment magnitude scale of the earthquake [31–33]:

3. Model Setup

3.1. Conceptual Model

ABAQUS was used as the numerical tool for simulating our fault-surrounding rock-fracturing region system (Figure 2) [27]. A 2D plane strain model was developed to simulate the system. As shown in Figure 3, the size of the basic model was 6,000 m × 3000 m. The length of the fault was 1000 m. Its dip angle was 60°. The distance between the top of the fault and the top boundary of the model was 1000 m. The hydraulic fracturing region was simulated by an elliptic zone with the major axis of 300 m and minor axis of 100 m. In the basic model, we designed 12 fracturing regions (TML1, TM, TMR1, FM, FT, FB, HM, HT, HB, BML1, BM, and BMR1) as shown in Figure 3. Among them, regions TML1, TM, and TMR1 were in the overlying strata of the fault. The horizontal distance between region TM and TML1 or TMR1 was 100 m. Among them, region TM was just above the fault. The distance between the bottom of region TM and the top of the fault was 300 m. Regions TML1, TM, and TMR1 were in underlying strata of the fault. The horizontal distance between region BM and BML1 or BMR1 was 100 m. Among them, region BM was just under the fault. The distance between the top of region TM and the bottom of the fault was 300 m. Regions FM, FT, and FB were in the footwall of the fault. Among them, region HM was in the middle of the footwall. Regions HM, HT, and HB were in the hanging wall of the fault. Among them, region HM was in the middle of the hanging wall. The horizontal distance between region FM and the fault was 500 m. The aim of the design was to investigate the influence of spatial distribution of hydraulic fracturing region on the stability of the fault. The material properties of the model were mainly obtained from the logs and well data from the field as well as laboratory experiments, as shown in Table 1.

Figure 3: The basic model of hydraulic fracturing in EGS engineering.

Table 1: Physical parameters of the rock mass.

Property
Density (kg/m³)
Young’s modulus (GPa)
Poisson’s ratio (—)
Porosity (—)
Friction coefficient (—)
Permeability (m/s)

Granite
2700
46
0.14
—
—
—

Hydraulic fracturing region
2100
6
0.35
0.36
—
1 × 10⁻⁷

Fault
2100
10
0.3
—
—
—

The left and right sides of the model were fixed in the horizontal direction, and the bottom was fixed in the vertical and horizontal direction. The top was allowed to deform freely.

The stress within the fault-surrounding rock-fracturing region system should be equal to the lithostatic stress before fluid injection for hydraulic fracturing. Pore pressure within fracturing region should be hydrostatic, with the initial value of 0 MPa. The initial stresses should be geostatic by the following steps: (1) the gravity was applied as an external force to obtain the geostatic stresses under the applied boundary conditions of the model; (2) the geostatic stresses were then inputted to the model as initial internal force before gravity; and (3) repeated step 2 until a mechanical equilibrium was reached between the initial stress fields and the gravity.

3.2. Numerical Testing

3.2.1. Fluid Injecting into the Hydraulic Fracturing Regions in the Basic Model

After the geostatic stresses balance, fluid was injected into the ellipse region to simulate hydraulic fracturing. According to the hydraulic fracturing test in Qiabuqia Geothermal Field, we injected fluid at a constant pressure of 100 MPa to the model plane over 216000 s (60 hours or 2.5days).

3.2.2. Location Change of the Fracturing Regions in the Vertical and Horizontal Directions

In order to investigate the influence of the location of the fracturing region on the fault’s stability, we designed more fracturing regions by changing the regions, as shown in Figure 4.

Figure 4: Comprehensive design of the fracturing regions in EGS engineering.

Firstly, we investigated fracturing regions in the overlying strata. We first changed the location of the regions in a vertical direction. By moving the regions TM, TML1, and TMR1 down 200 m, regions TB, TBL1, and TBR1 were obtained, respectively. By moving the regions TM, TML1, and TMR1 up 200 m, regions TT, TTL1, and TTR1 were obtained, respectively. We then changed the location of the regions in a horizontal direction. By moving region TBL1 400 m to the left each time, regions TBL2 and TBL3 were obtained. By moving TML1 400 m to the left each time, regions TML2 and TML3 were obtained. By moving TTL1 400 m to the left each time, regions TTL2 and TTL3 were obtained. Similarly, regions TBR2, TBR3, TMR2, TMR3, TTR2, and TTR3 can be obtained by moving regions TBR1, TMR1, and TBR1 400 m to the right each time.

Secondly, we investigated the location change of the fracturing regions in underlying strata. Using the similar method above, we obtained the regions BBL3, BBL2, BBL1, BB, BBR3, BBR2, BBR1, BTL3, BTL2, BTL1, BT, BTR3, BTR2, BTR1, BML3, BML2, BMR3, and BMR2.

Thirdly, we investigated the location change of the fracturing regions in the footwall and hanging wall of the fault. By moving regions FT, FM, and FB 400 m to the left each time, the regions FTL, FML, and FBL were obtained, respectively. By moving regions FT, FM, and FB 400 m to the right each time, the regions FTR, FMR, and FBR were obtained, respectively. Among them, regions FTR, FMR, and FBR were closer to the fault. By moving regions HT, HM and HB 400 m to the left each time, the regions HTL, HML, and HBL were obtained, respectively. By moving regions HT, HM, and HB 400 m to the right each time, the regions HTR1, HMR1, and HBR1 were obtained, respectively. Among them, regions HTL, HML, and HBLwere closer to the fault.

Above all, a total of 60 hydraulic fracturing regions were deigned (Figure 4).

3.2.3. Injection Pressure Change of the Fracturing Regions

Injection pressure is an important engineering factor that affects the stress perturbation in the fracturing region and adjacent rocks. It directly affected the induced fault slippage and the seismic moment magnitude . Based on the design of the fracturing regions in Figure 3, we change the injection pressure to 120 MPa and 80 MPa to study its influences on fault stability. The injection time remained constant at 60 hours.

Above all, a total of 180 numerical experiments were designed. After the finish of the injection, the slippage of the fault along the interface can be obtained. Based on the slippage, the seismic moment magnitude can be calculated.

4. Results and Discussion

4.1. Fracturing Region in the Overlying Strata of the Fault

4.1.1. Fracturing regions of TM, TML1, and TMR1

Figure 5 shows the fluid-injection-induced fault sliding distance when the fracturing regions were in overlying strata. As we can see, the induced slips of the three injection schemes were different. When the fracturing region was TM or TMR1, the value of induced sliding distance was positive. It revealed that the fault slid upward along the fault plane. When the fracturing region was TML1, however, the value of induced sliding distance was negative. The fault slid downward along the plane. The detailed results are summarized as follows. (1) When the fracturing region was TM (the black line with the square symbol), the fault slid upward along the whole plane. The upper part slid more significantly. At 150 m of the plane, the sliding distance was the largest, about 23.00 mm. From 50 m to 350 m, the sliding distances were all larger than 20 mm with an average sliding distance of 21.67 mm. Along the whole fault plane (0 m–1000 m), the average sliding distance was 10.14 mm. According to equations (5) and (6), an 2.8047 earthquake can be induced. (2) When the fracturing region was TMR1 (the blue line with the triangle symbol), the induced fault slip displayed a similar tendency (upward along the whole plane). However, the significant slip region moved down to the middle part. At 370 m of the plane, the sliding distance was the largest, about 26.42 mm. From 270 m to 470 m, the sliding distances were all larger than 22.93 mm, with an average sliding distance of 22.02 mm. Along the whole plane, the average sliding distance was 14.53 mm, and the of the induced earthquake was 2.9089. (3) When the fracturing region was TML1 (the red line with the circle symbol), the induced fault sliding distance displayed a contrary tendency (downward along the whole plane). The significant slip region still occurred at the upper part of the fault. At 70 m of the plane, the sliding distance was the largest, about 11.26 mm. From 20 m to 120 m, the sliding distances were all larger than 8.33 mm, with an average sliding distance of 10.60 mm. Along the whole plane, the average sliding distance was 6.46 mm with the of an induced earthquake of 2.6742. Note that region TML1 was on the left of TM (over the footwall of the fault) and region TMR1 was on the right of TM (over the hanging wall of the fault). The results may reveal that injecting fluid above the hanging wall of the fault may induce a larger earthquake than injecting fluid above the footwall of the fault.

Figure 5: Fault slip along the fault plane when the fracturing region was TM, or TML1, or TMR1.

4.1.2. Fracturing Regions Changed in the Vertical Direction

Based on the regions of TM, TML1, and TMR1, we studied the influence of location change of the fracturing region in the vertical direction on the induced fault sliding distances and earthquakes.

From Figure 6, we can see the induced fault sliding distance when the fracturing region changed from TM to TB and TT. When the fracturing region was changed to TB, the induced fault sliding distance increased. As was shown in the red line with the hollow circle symbol, the largest sliding distance was up to 78.24 mm. From 50 m to 350 m, the average sliding distance was 62.94 mm, over three times that of region TM. It showed that the upper part of the fault would slide violently. Along the whole plane, the average sliding distance was 30.94 mm, and the of the induced earthquake was up to 3.1277. When the fracturing region was directly over the top of the fault, it was likely to induce a felt earthquake. When the fracturing region was changed to TT, the induced fault sliding distance decreased. As we can see from the blue line with the solid circle symbol, the largest sliding distance was reduced to 15.25 mm. Along the whole plane, the average sliding distance was only 5.80 mm, and the of the induced earthquake was 2.6430. Note that the distance between the bottom of region TB and the top of the fault was 100 m, and the distance between the bottom of region TT and the top of the fault was 500 m. When the fracturing region changed from TB to TT, the average induced sliding distance reduced by more than 80 percent (from 30.94 mm to 5.80 mm). It revealed that the induced fault sliding distance and the moment magnitude of earthquake decrease rapidly with an increase of the distance between the region and fault.

Figure 6: Fault slip along the fault plane when the fracturing region was TT, or TM, or TB.

From Figure 7, we can see the induced fault sliding distance when the fracturing region changed from TML1 to TBL1 and TTL1. When the fracturing region changed to TBL1 (the red line with the hollow circle symbol), the induced fault sliding distance also increased. The largest sliding distance was 16.24 mm, and the average sliding distance was 9.74 mm. The of the induced earthquake was 2.6742. When the fracturing region changed to TTL1 (the blue line with the solid circle symbol), the induced fault sliding distance decreased. The largest sliding distance was only 5.52 mm, with an average sliding distance of 3.44 mm. The of the induced earthquake was 2.4917. When the fracturing region changed from TBL1 to TTL1, the average induced sliding distance reduced by about 65 percent (from 9.74 mm to 3.44 mm). The decreased range was less than that when the fracturing region changed from TB to TT.

Figure 7: Fault slip along the fault plane when the fracturing region was TTl1, or TMl1, or TBl1.

From Figure 8, we can see the induced fault sliding distance when the fracturing region changed from TMR1 to TBR1 and TTR1. When the fracturing region changed from TMR1 to TBR1, the sliding mode became different. As shown from the red line with the hollow circle symbol in Figure 8, the fault slip along the whole plane was presented as two parts: sliding downward from 0 m to 140 m and upward from 140 m to 1000 m. It caused the fault to be compressed at 140 m. Besides the change in the sliding mode, the amount of sliding distance changed obviously. The largest sliding distance was up to 44.37 mm, occurring at 470 m. The average sliding distance along the whole plane was up to 22.12 and the of the induced earthquake was up to 3.0305. When the fracturing region changed from TMR1 to TBR1, the sliding mode remained unchanged. The fault slid upward along the whole plane. The largest sliding distance was 10.92 mm, and the average sliding distance was 6.85 mm. The of the induced earthquake was 2.6911. When the fracturing region changed from TBR1 to TTR1, the average induced sliding distance reduced by about 70 percent (from 22.12 mm to 6.85 mm). The decreased range was also less than that when the fracturing region changed from TB to TT.

Figure 8: Fault slip along the fault plane when the fracturing region was TTR1, or TMR1, or TBR1.

The above results revealed that the induced sliding distance decreased with the vertical distance between the fault and fracturing region. The influence of vertical variation of the fracturing region on the fault slip was largest when the regions were above the top of the fault.

4.1.3. Fracturing Regions Changed in the Horizontal Direction

In this section, we studied the influence of horizontal variation of the fracturing region’s location on induced fault slip when the fracturing regions were in overlying strata of the fault.

Figure 9 shows the induced fault sliding distance of the fracturing region of TBL1-TBL3 and TBR1-TBR3. When the fracturing region shifted away from TBL1 (black dotted line with the square symbol) to TBL2 (red line with the circle symbol) and TBL3 (blue line with the triangle symbol), the tendency of induced fault slip remained similar, but the sliding distance significantly reduced. When the fracturing region was TBL2, the largest sliding distance was reduced to 3.15 mm, and the average sliding distance was only 2.52 mm. The of the induced earthquake was 2.4016. When the fracturing region was TBL3, the largest sliding distance was reduced to 0.64 mm, and the average sliding distance was 0.49 mm. The of the induced earthquake was only 1.9275. When the fracturing region shifted away from TBR1 (green dotted line with the hollow square symbol) to TBR2 (purple line with the hollow circle symbol) and to TBR3 (orange line with the hollow triangle symbol), however, the induced fault slip presented different tendencies. When the fracturing region shifted away from TBR1 to TBR2, the slip along plane was also presented in two sections. However, the zone of sliding downward of the fault became larger. As we can see, the fault slid downward from 0 m to 600 m. It resulted in the fault to be compressed at 600 m. When the fracturing region was shifted to TBR3, however, the fault slid downward along the whole plane.When the fracturing region was TBR2, the largest sliding distance was 5.60 mm, and the average sliding distance was 3.22 mm. The of the induced earthquake was only 2.4726. When the fracturing region was TBR3, the largest sliding distance was 3.50 mm, and the average sliding distance was 2.48 mm. The of the induced earthquake was only 2.3970. We knew that when the fracturing region TBL1 and TBR1 shifted away from the fault in the horizontal direction, the induced fault slip and seismicity both decreased.The horizontal variation of the fracturing region can affect both the slip pattern and the amount of sliding distance. Comparing the results of TBR2 with TBL2 (or TBR3 with TBR3), we knew that injecting fluid over the hanging wall of the fault was more likely to cause a larger fault slip than injecting fluid over the footwall under the same conditions.

Figure 9: Fault slip along the fault plane when the region of TB changed in the horizontal direction.

Figure 10 shows the induced fault sliding distance of the fracturing region of TML1-TML3 and TMR1-TMR3. When the fracturing region shifted away from TML1 (black dotted line with the square symbol) to TML2 (red line with the circle symbol) and TML3 (blue line with the triangle symbol), the variation characteristics of the induced fault slip were similar to that when the fracturing region shifted away from TBL1 to TBL2 and TBL3. When the fracturing region was TML2, the largest sliding distance was 3.39 mm, and the average sliding distance was 2.51 mm. The of the induced earthquake was only 2.4005. When the fracturing region was TML3, the largest sliding distance was 0.62 mm, and the average sliding distance was 0.50 mm. The of the induced earthquake was 1.9333. When the fracturing region shifted away from TMR1 (green dotted line with the hollow square symbol) to TMR2 (purple line with the hollow circle symbol) and TMR3 (orange line with the hollow triangle symbol), the variation characteristics of the induced fault slip were also similar to that when the fracturing region shifted away from TBR1 to TBR2 and TBR3. When the fracturing region was TMR2, the largest sliding distance was 4.10 mm, and the average sliding distance was 2.65 mm. The of the induced earthquake was 2.4162. When the fracturing region was TMR3, the largest sliding distance was 2.72 mm, and the average sliding distance was 1.88 mm. The of the induced earthquake was 2.3168. We knew that when the fracturing region TML1 and TMR1 shifted away from the fault in the horizontal direction, the induced sliding distance and seismicity decreased.The horizontal variation of the fracturing region would affect both the slip pattern and the amount of sliding distance. Comparing the results of TMR2 with TML2 (or TMR3 with TMR3), we knew that injecting fluid over the hanging wall of the fault was more likely to cause a larger fault slip than injecting fluid over the footwall under the same conditions.

Figure 10: Fault slip along the fault plane when region TM changed in the horizontal direction.

Figure 11 shows the induced fault sliding distance of the fracturing region of TTL1-TTL3 and TTR1-TTR3. When the fracturing region shifted away from TTL1 (black dotted line with the square symbol) to TTL2 (red line with the circle symbol) and TTL3 (blue line with the triangle symbol), the variation characteristics of the induced fault slip were also similar to that when the fracturing region shifted away from TBL1 to TBL2 and TBL3. When the fracturing region was TTL2, the largest sliding distance, average sliding distance, and of the induced earthquake were 2.64 mm, 1.80 mm, and 2.3042, respectively.When the fracturing region was TTL3, the largest sliding distance, average sliding distance, and of the induced earthquake were 0.44 mm, 0.35 mm, and 1.8301, respectively. When the fracturing region was TTR2, the largest sliding distance, average sliding distance, and of the induced earthquake were 2.45 mm, 1.85 mm, and 2.3121, respectively. When the fracturing region was TTR3, the largest sliding distance, average sliding distance, and of the induced earthquake were 1.71 mm, 1.15 mm, and 2.1745, respectively. We knew that when the fracturing region TTL1 and TTR1 shifted away from the fault in the horizontal direction, the induced sliding distance and seismicity decreased. The horizontal variation of the fracturing region would affect both the slip pattern and the amount of sliding distance. Comparing the results of TTR2 with TTL2 (or TTR3 with TTR3), we knew that injecting fluid over the hanging wall of the fault was more likely to cause a larger fault slip than injecting fluid over the footwall under the same conditions.

Figure 11: Fault slip along the fault plane when the region TT changed in the horizontal direction.

The above results revealed that, when the fracturing regions were in overlying strata of the fault, the horizontal variation of the fracturing region may affect both the slip pattern and the amount of the sliding distance of the fault. The sliding distance decreased along with the increase of the horizontal distance between the fault and fracturing region. Injecting fluid over the hanging wall of the fault was more likely to cause a larger fault slip than injecting fluid over the footwall under the same conditions. Therefore, for the overlying strata of fault, the area close to the top of the hanging wall of the fault was the more hazardous injection area.

4.2. Fracturing Region in the Footwall and Hanging Wall of the Fault

Figure 12 shows the induced fault sliding distance when the fracturing region was in the footwall and hanging wall of the fault. The horizontal distances between the regions (FT, FM, FB, HT, HM, or HB) and the fault were 500 m. When the fracturing region was FT (black line with the square symbol), the fault slid downward along the whole plane. The largest sliding distance was 5.66 mm, and the average sliding distance was 4.37 mm. The of the induced earthquake was 2.5610. When the fracturing region was HT (green line with the hollow square symbol), the induced fault slip along the whole plane presented in two sections: from 0 m to 810 m, fault slid downward; from 810 m to 100 m, fault slid upward. It caused the fault to be compressed at 810 m. The largest sliding distance, average sliding distance, and of the induced earthquake were 8.62 mm, 5.50 mm, and 2.6276, respectively. When the fracturing region was FM (red line with the circle symbol), the fault also slid downward along the whole plane. The largest sliding distance, average sliding distance, and of the induced earthquake were 7.50 mm, 5.19 mm, and 2.6108, respectively. When the fracturing region was HM (purple line with the hollow circle symbol), the fault slid downward along the whole plane. The largest sliding distance, average sliding distance, and of the induced earthquake were 8.51 mm, 5.91 mm, and 2.6484, respectively. When the fracturing region was FB, the fault slip along the whole plane became complicated. As shown in the blue line with the triangle symbol in Figure 12, the fault slid upward from 0 m to 240 m and downward from 240 m to 1000 m. This caused the fault to be tensioned at the point of 240 m, which was prone to failure. The largest fault sliding distance was 6.04 mm, and the average sliding distance was 3.94 mm. The of the induced earthquake was 2.5310. When the fracturing region was HB, the fault also slid downward along the whole plane (orange line with the hollow triangle symbol in Figure 12). The largest sliding distance, average sliding distance, and of the induced earthquake were 5.38 mm, 4.81 mm, and 2.5888, respectively.

Figure 12: Fault slip along the fault plane when the fracturing region was FT, TB, FM, HT, HM, or HB.

When the horizontal distances between the regions and the fault were reduced to 100 m, the induced fault sliding distances were shown in Figure 13. The sliding pattern and the amount of sliding distance of fault both changed remarkably. When the fracturing region was FTR (located on the right of region FT), the induced fault slip along the whole plane (black line with the square symbol) presented in two sections: from 0 m to 60 m, the fault slid upward; from 60 m to 1000 m, the fault slid downward. It caused the fault to be tensioned at the point of 60 m and prone to failure. The largest fault sliding distance was 33.81 mm, and the average sliding distance was 21.43 mm. The of the induced earthquake was 3.0214. When the fracturing region was FMR (red line with the circle symbol), the fault slip also presented in two sections, but the range of sliding upward increased (0 m–380 m). The largest sliding distance, average sliding distance, and of the induced earthquake were 75.35 mm, 32.38 mm, and 3.1409, respectively. When the fracturing region was FBR (blue line with the triangle symbol), the range of fault that slid upward continuously increased (0 m–660 m). The largest sliding distance, average sliding distance, and of the induced earthquake were 123.32 mm, 60.01 mm, and 3.3195, respectively. When the fracturing region was HTL (located on the left of region FT), the induced fault slip showed an opposite trend (green line with the hollow square symbol). The upper part of the fault slid downward (0 m–370 m) and the lower part of fault slid upward (370 m–1000 m). It caused the fault to be compressed at the point of 370 m. The largest sliding distance, average sliding distance, and of the induced earthquake were 115.25 mm, 45.77 mm, and 3.2411, respectively. When the fracturing region was HML (purple line with the hollow circle symbol), the range of sliding upward of the fault increased to 0 m–370 m. The largest sliding distance, average sliding distance, and of the induced earthquake were 52.67 mm, 22.75 mm, and 3.0387, respectively. When the fracturing region was HBL (the orange line with the hollow triangle symbol), the induced fault slip presented only one section (downward). The largest sliding distance, average sliding distance, and of the induced earthquake were 32.84 mm, 19.92 mm, and 3.0002, respectively.

Figure 13: Fault slip along the fault plane when the fracturing region was FTR, TBR, FMR, HTL, HML, or HBL.

From the above results, we knew that when the fracturing region was close to the fault, no matter it was in the footwall or hanging wall, injecting fluid into it can induce a felt earthquake. For the fracturing region in the footwall of the fault, the induced slip became larger when the region was closer to the lower part of the fault. For the fracturing region in the hanging wall of the fault, on the contrary, the induced slip became larger when the region was closer to the upper part of the fault. It may reveal that the regions that were close to the lower part of the fault in the footwall and to the upper part of the fault in the hanging wall were more dangerous areas.

When the horizontal distances between the fracturing regions and the fault were increased to 900 m, the induced fault sliding distances were shown in Figure 14. With the fracturing region becoming far away from the fault, the sliding pattern of the fault became simple. When the fluid was injected into the six different regions (FTL, FML, FBL, HTR, HMR, or HBR), the fault all slid downward along the whole plane. When the fracturing region was FTL, the largest sliding distance, average sliding distance, and of the induced earthquake were 1.35 mm, 1.12 mm, and 2.1668, respectively. When the fracturing region was FML, the largest sliding distance, average sliding distance, and of the induced earthquake were 2.56 mm, 11.74 mm, and 2.2944, respectively. When the fracturing region was FBL, the largest sliding distance, average sliding distance, and of the induced earthquake were 3.08 mm, 1.86 mm, and 2.3137, respectively. When the fracturing region was HTR, the largest sliding distance, average sliding distance, and of the induced earthquake were 4.50 mm, 3.21 mm and 2.4717, respectively. When the fracturing region was HMR, the largest sliding distance, average sliding distance, and of the induced earthquake were 3.41 mm, 2.60 mm, and 2.4107, respectively. When the fracturing region was HBR, the largest sliding distance, average sliding distance, and of the induced earthquake were 2.40 mm, 2.16 mm, and 2.3570, respectively. The average sliding distances of six different fracturing regions were all smaller than 3.21 mm. From the above results, we knew that when the fracturing region was far away from the fault, the induced fault slip and seismicity were both small. When the fracturing region was in the footwall of the fault, the induced slip also became larger if the region was closer to the lower part of the fault. When the fracturing region was in the hanging wall of the fault, the induced slip became larger if the region was closer to the upper part of the fault. The same conclusion can be obtained as in the previous paragraph.

Figure 14: Fault slip along the fault plane when the fracturing region was FTL, TBL, FML, HTR, HMR, or HBR.

4.3. Fracturing Region in Underlying Strata of Fault

4.3.1. Fracturing Regions of BM, BML1, and BMR1

Figure 15 shows the fluid-injection-induced fault sliding distance when the fracturing regions were in underlying strata. When the fracturing region was BM (the black line with the square symbol), the fault slid upward along the whole plane, and the sliding distance along the lower part was large. From 500 m to 1000 m, the average sliding distance was 11.30 mm. From 0 m to 500 m, the average sliding distance was only 5.65 mm. Along the whole fault plane, the average sliding distance was 8.33 mm. The of the induced earthquake was 2.7478. When the fracturing region was BML1 (the red line with the square symbol), the fault also slid upward along the whole plane, and the sliding distance became larger than that of BM. The sliding distance along the middle section of the fault plane was much larger. The largest fault sliding distance was 39.08 mm, occurring at 510 m. The average sliding distance was 24.64 mm, and the of the induced earthquake was 3.0618. When the fracturing region was BMR1 (the blue line with the triangle symbol), the fault slid downward along the whole plane. The largest sliding distance, average sliding distance, and of the induced earthquake were 8.82 mm, 5.75 mm, and 2.6405, respectively. The results revealed that injecting fluid under the footwall of the fault may induce a much larger fault sliding distance than injecting fluid under hanging wall of the fault.

Figure 15: Fault slip along the fault plane when the fracturing region was BM, BML1, or BMR1.

4.3.2. Fracturing Regions Changed in the Vertical Direction

We then studied the influence of location change of the fracturing regions in the vertical direction on the induced fault sliding distances and earthquakes.

Figure 16 shows the induced fault sliding distance when the fracturing region changed from BM to BT and BB. When the fracturing region changed to BT (closer to the fault), the induced fault sliding distance increased. As shown from the blue line with the solid square symbol, the largest sliding distance was 23.68 mm (at 970 m). The average sliding distance was 9.84 mm, and the of the induced earthquake was 2.7960. When the fracturing region changed to BB (farther to the fault), the induced fault sliding distance decreased (red line with the hollow square symbol). The largest sliding distance was reduced to 10.10 mm. The average sliding distance was reduced to 6.27 mm, and the of the induced earthquake was 2.6656. When the fracturing region changed from BT to BB, the average induced sliding distance reduced by 36 percent (from 9.84 mm to 6.27 mm).

Figure 16: Fault slip along the fault plane when the fracturing region was BT, BM, or BB.

Figure 17 shows the induced fault sliding distance when the fracturing region changed from BML1 to BTL1 and BBL1. When the region changed to BTL1 (closer to the fault), the induced fault sliding distance increased (the blue line with the solid circle symbol). As shown, the largest sliding distance was increased to 59.72 mm (at 490 m). The average sliding distance was 35 mm, and the of the induced earthquake was 3.1634. When the fracturing region changed to BBL1 (farther to the fault), the induced fault sliding distance decreased (red line with the hollow circle symbol). The largest sliding distance was reduced to 25.40 mm. The average sliding distance was 15.97 mm, and the of the induced earthquake was reduced to 2.9362. When the fracturing region changed from BTL1 to BBL1, the average induced sliding distance reduced by 54 percent (from 9.84 mm to 6.27 mm).

Figure 17: Fault slip along the fault plane when the fracturing region was BTL1, BML1, or BBL1.

Figure 18 shows the induced fault sliding distance when the fracturing region changed from BMR1 to BTR1 and to BBR1.When the fracturing region changed to BTR1 (closer to the fault), the induced fault sliding distance increased. As shown from the blue line with the solid triangle symbol, the largest sliding distance was increased to 12.53 mm (at 830 m). The average sliding distance was increased to 8.46 mm, and the of the induced earthquake was 2.7523. When the fracturing region changed from BMR1 to BBR1 (becoming farther to the fault), the induced fault sliding distance decreased (red line with the hollow triangle symbol). The largest sliding distance was reduced to 4.31 mm. The average sliding distance was reduced to 3.03 mm, and the of the induced earthquake was 2.4550. When the fracturing region changed from BTR1 to BBR1, the average induced sliding distance reduced by 52 percent (from 9.84 mm to 6.27 mm).

Figure 18: Fault slip along the fault plane when the fracturing region was BTR1, BMR1, or BBR1.

From the above results, we knew that when region BM, BML1, and BMR1 moved up 200 m to BT, BTL1, and BTL1, respectively, the induced fault slip and seismicity all increased. When region BML1 moved up to BTL1, the increase was the largest. Taking the average induced sliding distance, for example, it increased by 10.36 mm (from 24.64 mm to 35.00 mm). It revealed that the influence of vertical variation of the fracturing region’s location on fault slip was largest when the regions were under the footwall of the fault. For underlying strata of the fault, the area that was close to the bottom of the footwall of the fault was a hazardous injection area.

4.3.3. Fracturing Regions Changed in the Horizontal Direction

In this section, we studied the influence of the fracturing region changing in the horizontal direction on induced fault slip when the fracturing regions were in underlying strata of the fault.

Figure 19 shows the induced fault sliding distance of the fracturing region of BBL1-BBL3 and BBR1- BBR3. As shown, when the fracturing region shifted away from BBL1 (black dotted line with the square symbol) to BBL2 (red line with the circle symbol), the induced fault slip presented a similar tendency (upward). However, the amount of the sliding distance changed significantly. The largest sliding distance was reduced to 10.05 mm, occurring at 370 m. The average sliding distance was 6.00 mm. The of the induced earthquake was 2.6528. When the fracturing region changed to BBL3 (blue line with the triangle symbol), the induced fault slip turned into the opposite tendency (downward). The largest sliding distance was only 3.24 mm, and the average sliding distance was 2.67 mm. The of the induced earthquake was 2.4184. When the fracturing region shifted away from BBR1 (green dotted line with the hollow square symbol) to BBR2 (purple line with the hollow circle symbol), the induced fault slip also presented a similar tendency (downward). The largest sliding distance was 4.61 mm, and the average sliding distance was 3.40 mm. The of the induced earthquake was 2.4883. When the fracturing region changed to BBR3 (orange line with the hollow triangle symbol), the induced fault slip became complicated and presented two sections: from 0 m to 600 m, the fault slid upward; from 600 m to 1000 m, it slid downward. The largest sliding distance was 2.69 mm, and the average sliding distance was 1.14 mm. The of the induced earthquake was 2.1720. We knew that when the fracturing region BBL1 and BBR1 shifted away from the fault in the horizontal direction, the induced sliding distance and seismicity decreased. The horizontal variation of the fracturing region would affect both the slip pattern and the amount of sliding distance. Comparing the results of BBR2 with BBL2 (or BBR3 with BBR3), we knew that injecting fluid under the footwall was more likely to cause a larger fault sliding distance than injecting fluid under the hanging wall of the fault under the same conditions.

Figure 19: Fault slip along the fault plane when region BB changed in the horizontal direction.

Figure 20 shows the induced fault sliding distance of the fracturing region of BML1-BML3 and BMR1-BMR3. When the fracturing region shifted away from BML1 (black dotted line with the square symbol) to BML2 (red line with the circle symbol), the induced fault slip became complicated and presented two sections: from 0 m to 800 m, the fault slid upward; from 800 m to 1000 m, the fault slid downward. The largest sliding distance, average sliding distance, and of the induced earthquake were 8.70 mm, 4.56 mm, and 2.5733, respectively. When the fracturing region was BML3, the induced fault slip (blue line with the triangle symbol) presented an opposite tendency with that when the region was BML1. The fault slid downward along the whole plane and the maximum sliding distance occurred at 610 m, about 2.93 mm. The average sliding distance was 2.45 mm, and the of the induced earthquake was 2.3935. When the fracturing region shifted away from BMR1 (green dotted line with the hollow square symbol) to BMR2 (purple line with the hollow circle symbol), the induced fault slip presented a similar tendency (downward). The largest sliding distance was 4.85 mm, and the average sliding distance was 3.86 mm. The of the induced earthquake was 2.5251. When the fracturing region changed to BMR3 (orange line with the hollow triangle symbol), the induced fault slip also presented two sections: from 0 m to 380 m, the fault slid upward; from 380 m to 1000 m, the fault slid downward. The largest sliding distance was 1.60 mm, and the average sliding distance was 0.90 mm. The of the induced earthquake was 2.1350. We knew that when the fracturing region BML1 and BMR1 shifted away from the fault in the horizontal direction, the induced sliding distance and seismicity decreased.The horizontal variation of the fracturing region would affect both the slip pattern and the amount of sliding distance. Comparing the results of BMR2 with BBL2 (or BMR3 withBML3), We knew that injecting fluid under the footwall was more likely to cause a larger fault sliding distance than injecting fluid under the hanging wall of the fault under the same conditions.

Figure 20: Fault slip along the fault plane when region BM changed in the horizontal direction.

Figure 21 shows the induced fault sliding distance of the fracturing region of BTL1-BTL3 and BTR1-BTR3. When the fracturing region shifted away from BTL1 (black dotted line with the square symbol) to BTL2 (red line with the circle symbol), the induced fault slip also became complicated and presented two sections: from 0 m to 570 m, the fault slid upward; from 570 m to 1000 m, the fault slid downward. The largest sliding distance, average sliding distance, and of the induced earthquake were 9.13 mm, 4.01 mm, and 2.5361, respectively. When the fracturing region was BML3 (blue line with the triangle symbol), the fault slid downward along the whole plane. The maximum sliding distance was 2.47 mm. The average sliding distance was 2.10 mm, and the of the induced earthquake was 2.3488. When the fracturing region shifted away from BTR1 (green dotted line with the hollow square symbol) to BMR2 (purple line with the hollow circle symbol), the induced fault slip presented a similar tendency (downward). The largest sliding distance was 4.41 mm, and the average sliding distance was 3.64 mm. The of the induced earthquake was 2.5081. When the fracturing region changed to BTR3 (orange line with the hollow triangle symbol), the induced fault slip also presented two sections: from 0 m to 140 m, the fault slid upward; from 140 m to 1000 m, the fault slid downward. The largest sliding distance was 1.83 mm, and the average sliding distance was 1.10 mm. The of the induced earthquake was 2.1616. We knew that when the fracturing regions BTL1 and BTR1 shifted away from the fault in the horizontal direction, the induced sliding distance and seismicity decreased. The horizontal variation of the fracturing region would affect both the slip pattern and the amount of sliding distance. Comparing the results of BTR2 with BBL2 (or BTR3 withBML3), we knew that injecting fluid under the footwall was more likely to cause a larger fault sliding distance than injecting fluid under the hanging wall of the fault under the same conditions.

Figure 21: Fault slip along the fault plane when region BT changed in the horizontal direction.

The above results revealed that, when the fracturing regions were in underlying strata of the fault, the horizontal variation of the fracturing region may affect both the slip pattern and the amount of sliding distance. The sliding distance decreased along with the increase of the horizontal distance between the fault and fracturing region. Injecting fluid under the footwall of the fault was more likely to cause a larger fault slip than injecting fluid under the hanging wall of the fault under the same conditions. Therefore, for underlying strata of fault, the area close to the bottom of footwall of the fault was the more hazardous injection area.

4.4. Hazardous Injection Area

Risk assessment of the induced earthquakes in the process of hydraulic fracturing for establishing EGS is essential. Therefore, systematic research on seismic activities and spatial evolution characteristics is needed. Hazardous injection areas of induced felt earthquakes can be identified by numerical simulation. It is helpful for the seismic disaster prevention of such underground fluid-injection engineering. Based on the above analysis, the results of induced fault slip and seismicity were given in a tabular form, as shown in Table 2. The induced earthquakes with larger than 3.0 can be felt by people and may cause public panic. Therefore, we selected the fracturing region that induced earthquake as hazardous injection region. The hazardous injection area of the specific fault can be obtained by connecting all the hazardous regions. When the injection pressure was 100 MPa, the hazardous injection regions were TB, TBR1, FTR1, FMR1, RBR1, HTL1, HML1, HBL1, BTL1, and BML1. According to the ten hazardous injection regions, the hazardous injection area can be drawn, as shown in Figure 22. It revealed that when the injection pressure was 100 MPa, the hazardous injection regions were near the fault and concentrated on four areas: (1) above the top of the fault in underlying strata; (2) above the top of the hanging wall of the fault in underlying strata; (3) near the fault planes in both footwall and hanging wall; and (4) the bottom of the footwall of the fault in underlying strata.

Table 2: The induced maximal slippage, average slippage, and moment magnitude of fault when the pressure was 100 MPa.

Fracturing region
Maximal slippage (mm)
Average slippage (mm)
Moment magnitude ()

TTL1
5.52
3.44
2.4917

TTL2
2.64
1.80
2.3042

TTL3
0.44
0.35
1.8301

TT
15.25
5.80
2.6430

TTR1
10.92
6.85
2.6911

TTR2
2.45
1.85
2.3121

TTR3
1.71
1.15
2.1745

TML1
11.26
6.46
2.6742

TML2
3.39
2.51
2.4005

TML3
0.62
0.50
1.9333

TM
22.72
10.14
2.8047

TMR1
26.42
14.53
2.9089

TMR2
4.10
2.65
2.4162

TMR3
2.72
1.88
2.3168

TBL1
16.24
9.74
2.7931

TBL2
3.15
2.52
2.4016

TBL3
0.64
0.49
1.9275

TB
78.24
30.94
3.1277

TBR1
44.37
22.12
3.0305

TBR2
5.60
3.22
2.4726

TBR3
3.50
2.48
2.3970

FTL1
1.35
1.12
2.1668

FT
5.66
4.37
2.5610

FTR1
33.81
21.43
3.0214

FML1
2.56
1.74
2.2944

FM
7.50
5.19
2.6108

FMR1
75.35
32.38
3.1409

FBL1
3.08
1.74
2.2944

FB
6.04
3.94
2.5310

FBR1
123.32
60.01
3.3195

HTL1
115.25
45.77
3.2411

HT
8.62
5.50
2.6276

HTR1
4.50
3.21
2.4717

HML1
52.67
22.75
3.0387

HM
8.51
5.91
2.6484

HMR1
3.41
2.60
2.4107

HBL1
32.84
19.92
3.0002

HB
5.38
4.81
2.5888

HBR1
2.40
2.16
2.3570

BTL1
59.72
35.00
3.1634

BTL2
9.13
4.01
2.5361

BTL3
2.47
2.10
2.3488

BT
23.68
9.84
2.7960

BTR1
12.53
7.15
2.7036

BTR2
4.41
3.64
2.5081

BTR3
1.83
1.10
2.1616

BML1
39.08
24.64
3.0618

BML2
8.70
4.56
2.5733

BML3
2.93
2.45
2.3935

BM
13.41
8.33
2.7478

BMR1
8.82
5.75
2.6405

BMR2
4.85
3.86
2.5251

BMR3
1.60
0.90
2.1035

BBL1
25.40
15.97
2.9362

BBL2
10.05
6.00
2.6528

BBL3
3.24
2.67
2.4184

BB
10.10
6.27
2.6655

BBR1
4.31
3.43
2.4909

BBR2
4.61
3.40
2.4883

BBR3
2.69
1.14
2.1720

Figure 22: Hazardous injection area when injection pressure was 100 MPa.

When the injection pressure was increased to 120 MPa, the induced maximal slippage, average slippage, and moment magnitude of all fracturing regions were as shown in Table 3. According to the results, the hazardous injection area can be obtained (Figure 23). As shown, it contained 12 hazardous injection regions (TB, TBR1, TMR1, FTR1, FMR1, RBR1, HTL1, HML1, HBL1, BTL1, BML1, and BBL1). The hazardous injection area was enlarged with the increase of injection pressure. The spatial distribution of the regions was almost the same as that when the injection pressure was 100 MPa.

Table 3: The induced maximal slippage, average slippage, and moment magnitude of fault when the pressure was 120 MPa.

Fracturing region
Maximal slippage (mm)
Average slippage (mm)
Moment magnitude ()

TTL1
6.59 (80 m)
4.10
2.5425

TTL2
3.16 (110 m)
2.15
2.3556

TTL3
0.52 (420 m)
0.42
1.8829

TT
18.30 (120 m)
6.96
2.6958

TTR1
14.05 (400 m)
8.54
2.7550

TTR2
2.93 (90 m)
1.45
2.2416

TTR3
2.05 (320 m)
1.37
2.2252

TML1
13.56 (70 m)
7.77
2.7276

TML2
4.06 (160 m)
3.00
2.4521

TML3
0.74 (490 m)
0.59
1.9813

TM
34.05 (140 m)
14.75
2.9132

TMR1
34.70 (370 m)
19.92
3.0002

TMR2
4.90 (130 m)
2.67
2.4184

TMR3
3.25 (370 m)
2.25
2.3688

TBL1
19.51 (85 m)
11.70
2.8461

TBL2
3.78 (310 m)
3.02
2.4540

TBL3
0.76 (630 m)
0.59
1.9813

TB
102.62 (70 m)
40.97
3.2090

TBR1
57.25 (470 m)
28.95
3.1085

TBR2
6.74 (230 m)
3.86
2.5251

TBR3
4.20 (460 m)
2.98
2.4502

FTL1
1.62 (60 m)
1.34
2.2188

FT
6.80 (65 m)
5.24
2.6136

FTR1
40.77 (330 m)
25.78
3.0749

FML1
3.07 (110 m)
2.08
2.3461

FM
8.99 (240 m)
6.23
2.6637

FMR1
90.63 (220 m)
38.94
3.1943

FBL1
3.69 (310 m)
2.60
2.4107

FB
7.29 (510 m)
4.76
2.5858

FBR1
146.67 (410 m)
72.06
3.3725

HTL1
138.94 (590 m)
55.84
3.2987

HT
10.33 (505 m)
6.21
2.6627

HTR1
5.40 (685 m)
3.85
2.5243

HML1
62.65 (800 m)
27.32
3.0917

HM
10.21 (760 m)
7.09
2.7011

HMR1
4.10 (815 m)
3.12
2.4635

HBL1
39.29 (830 m)
23.61
3.0494

HB
6.41 (875 m)
5.79
2.6425

HBR1
2.89 (500 m)
2.60
2.4107

BTL1
74.91 (490 m)
44.24
3.2312

BTL2
11.20 (150 m)
4.97
2.5983

BTL3
2.20 (620 m)
1.32
2.2144

BT
27.94 (970 m)
12.37
2.8623

BTR1
15.06 (830 m)
10.07
2.8027

BTR2
5.30 (600 m)
4.38
2.5617

BTR3
2.96 (535 m)
2.52
2.4016

BML1
50.91 (540 m)
32.39
3.1410

BML2
11.27 (310 m)
5.86
2.6460

BML3
1.48 (770 m)
1.08
2.1563

BM
17.13 (630 m)
10.28
2.8087

BMR1
10.63 (895 m)
6.89
2.6928

BMR2
5.83 (750 m)
4.64
2.5784

BMR3
3.51 (615 m)
2.94
2.4463

BBL1
34.34 (560 m)
21.74
3.0255

BBL2
12.71 (280 m)
7.63
2.7224

BBL3
3.23 (150 m)
1.37
2.2252

BB
12.16 (570 m)
7.47
2.7162

BBR1
5.19 (910 m)
4.32
2.5577

BBR2
5.54 (770 m)
3.65
2.5089

BBR3
3.89 (730 m)
3.21
2.4717

Figure 23: Hazardous injection area when the pressure was 120 MPa.

When the injection pressure was decreased to 80 MPa, the induced maximal slippage, average slippage, and moment magnitude of all fracturing regions were as shown in Table 4. According to the results, the hazardous injection area can be obtained (Figure 24). As shown, it contained only five hazardous injection regions (TB, FMR1, RBR1, HTL1, and BTL1). The area was reduced with the decrease of injection pressure. When the injection pressure was 80 MPa, the hazardous injection regions were all near the fault and concentrated on the following areas:(1) above the top of the fault in underlying strata; (2) near the lower part of the fault in the footwall; (3) near the upper part of the fault in the hanging wall; and (4) under the bottom of the footwall of the fault in underlying strata.

Table 4: The induced maximal slippage, average slippage, and moment magnitude of fault when the pressure was 80 MPa.

Experiment no.
Fault F1

Maximal slippage (mm)
Average slippage (mm)
Moment magnitude ()

TTL1
4.46 (75 m)
2.78
2.4300

TTL2
2.12 (110 m)
1.44
2.2396

TTL3
0.36 (420 m)
0.29
1.7756

TT
12.20 (120 m)
4.64
2.5784

TTR1
8.77 (400 m)
5.50
2.6276

TTR2
1.96 (90 m)
0.97
2.1252

TTR3
2.05 (320 m)
0.92
2.1099

TML1
8.97 (70 m)
5.15
2.6086

TML2
2.72 (160 m)
2.01
2.3362

TML3
0.50 (490 m)
0.40
1.8687

TM
2058 (140 m)
8.70
2.7604

TMR1
19.15 (370 m)
10.60
2.8176

TMR2
3.29 (130 m)
1.80
2.3042

TMR3
2.18 (370 m)
1.51
2.2533

TBL1
12.98 (85 m)
7.79
2.7284

TBL2
2.52 (310 m)
2.01
2.3362

TBL3
0.51(630 m)
0.39
1.8614

TB
54.66 (70 m)
21.32
3.0199

TBR1
31.63 (470 m)
15.48
2.9272

TBR2
4.48 (230 m)
2.58
2.4084

TBR3
2.80 (460 m)
1.99
2.3333

FTL1
1.08 (60 m)
0.89
2.1003

FT
4.53 (65 m)
3.41
2.4892

FTR1
26.88 (330 m)
17.08
2.9557

FML1
2.05 (110 m)
1.39
2.2294

FM
6.01 (240 m)
4.16
2.5467

FMR1
59.91 (220 m)
25.82
3.0753

FBL1
2.47 (310 m)
1.74
2.2944

FB
4.80 (510 m)
3.06
2.4578

FBR1
97.51 (410 m)
47.05
3.2491

HTL1
88.70 (590 m)
34.97
3.1632

HT
6.91 (505 m)
4.16
2.5467

HTR1
3.60 (685 m)
2.56
2.4062

HML1
38.93 (800 m)
17.48
2.9624

HM
6.81 (760 m)
4.73
2.5839

HMR1
2.73 (815 m)
2.08
2.3461

HBL1
26.36 (830 m)
15.71
2.9315

HB
4.30 (875 m)
3.84
2.5236

HBR1
1.92 (500 m)
1.72
2.2910

BTL1
44.58 (490 m)
25.85
3.0757

BTL2
7.34 (150 m)
3.20
2.4708

BTL3
1.47 (620 m)
0.88
2.0970

BT
19.51 (970 m)
7.93
2.7335

BTR1
10.01 (830 m)
6.76
2.6873

BTR2
3.52 (600 m)
2.91
2.4433

BTR3
1.97 (535 m)
1.68
2.2842

BML1
27.56 (540 m)
17.11
2.9562

BML2
6.82 (150 m)
3.46
2.4934

BML3
1.28 (770 m)
0.72
2.0389

BM
10.30 (630 m)
6.62
2.6813

BMR1
7.03 (895 m)
4.56
2.5733

BMR2
3.88 (750 m)
3.09
2.4607

BMR3
2.34 (615 m)
1.96
2.3289

BBL1
16.71 (560 m)
10.58
2.8170

BBL2
7.35 (280 m)
4.56
2.5733

BBL3
2.15 (150 m)
0.91
2.1067

BB
8.05 (570 m)
5.10
2.6057

BBR1
3.44 (910 m)
2.42
2.3899

BBR2
3.68 (770 m)
2.87
2.4393

BBR3
2.59 (730 m)
2.14
2.3543

Figure 24: Hazardous injection areas when pressure was 80 MPa.

5. Conclusions

Injecting fluid into underground formations can exercise a great influence on the stability of the natural faults. When the fault is connected to the injection area, the main cause of fault reactivation and slip is the reduction of effective stress and the weakening of the physical and mechanical properties of the fault induced by fluid injection. When the fault is not connected to the injection area, the fault can also be reactivated by the solid stress perturbation transmitted from the injection area. Therefore, for an EGS project, it is necessary to estimate the stability of the fault that is not connected to the injection area during hydraulic fracturing. In this paper, the influence of fluid injection during hydraulic fracturing on fault’s stability was investigated by a 2D numerical model. The study took Qiabuqia Geothermal Field as a case. The research was focused on the following aspects: (1) the effect of the spatial distribution of the injection zone on induced fault slip and the magnitude of seismicity; (2) the influence of the change of distance between the injection zone and the fault on fault slip; (3) the dangerous injection area that threatens the fault’s stability; and (4) the methods to reduce the fluid-injection-induced fault slip and the seismicity. The results revealed that the spatial distribution of the fracturing region had a significant influence on the induced fault slip and of an induced earthquake. Delimiting hazardous injection area is very helpful for the selection of the location of the hydraulic fracturing region. It is useful for the prevention of fluid-injection-induced earthquakes. This study produced the following conclusions:

(1)When the fracturing region was in the overlying strata of fault, the horizontal variation of the region may affect both the slip pattern and the amount of the sliding distance. The induced sliding distance decreased with the increase of the horizontal distance between the fault and fracturing region. It also decreased with the vertical distance. The influence of vertical variation of the fracturing region’s location on fault slip was the largest when the regions were above the top of the fault. Injecting fluid over the hanging wall of the fault was more likely to cause a larger fault slip than injecting fluid over the footwall under the same conditions. For overlying strata of the fault, the area that closes to the top of the hanging wall of the fault was the more hazardous injection area.(2)When the fracturing regions were in underlying strata of the fault, the horizontal variation of the fracturing region may also affect both the slip pattern and the value of the sliding distance of the fault. The sliding distance decreased along with the horizontal distance (or vertical distance) between the fault and fracturing region. Injecting fluid under the footwall of the fault was more likely to cause a larger fault slip than injecting fluid under the hanging wall of the fault under the same conditions. For underlying strata, the area that closes to the bottom of footwall of the fault was the more hazardous injection area.(3)When the fracturing region was in the footwall or hanging wall of the fault, injecting fluid into the region can induce felt earthquake if it was close to the fault. For the fracturing region in the footwall of the fault, the induced slip became larger when the region was close to the lower part of the fault. For the fracturing region in the hanging wall of the fault, on the contrary, the induced slip became larger when the region was close to the upper part of the fault. It revealed that the region that was close to the lower part of the fault in the footwall and the region that was close to the upper part of the fault in the hanging wall were more dangerous areas. The hazardous injection area was significantly influenced by the injection pressure. The area would be increased with the increase of the pressure.(4)The hazardous injection regions that threaten the fault’s stability were near the fault and concentrated on four areas: (a) above the top of the fault in underlying strata; (b) above the top of the hanging wall of the fault in underlying strata; (c) near to the fault planes in both footwall and hanging wall; and (d) the bottom of the footwall of the fault in underlying strata.

Note that this study was only focused on one fault in the site. The above hazardous injection area was only determined by the stability of the specific fault. However, it usually has many faults in the site. It has at least two effects on the hazardous injection area. Firstly, the interactions between the faults can directly influence their stability and slippage during fluid injection. Secondly, for an EGS project, the hazardous injection area should be a superimposed area of the hazardous regions of all the faults. Therefore, it is necessary to study the hazardous injection area determined by all the major faults. This involves simulating several faults in a model. It is beyond the scope of this paper, and further research is needed.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

The authors are grateful for the financial support from the start-up fund for doctoral research of Suzhou University (no. 2019jb09 and no. 2019jb18) and Dual-ability teaching team project (Geological Engineering Dual-Energy Teaching Team, no. 2020XJSN06). The authors would like to express their sincere gratitude to Mahmoud from the University of Alberta for his help in proofreading this manuscript.

References

W. G. P. Kumari and P. G. Ranjith, “Sustainable development of enhanced geothermal systems based on geotechnical research—a review,” Earth-Science Reviews, vol. 199, Article ID 102955, 2019.

P. Bayer, L. Rybach, P. Blum, and R. Brauchler, “Review on life cycle environmental effects of geothermal power generation,” Renewable and Sustainable Energy Reviews, vol. 26, pp. 446–463, 2013.

A. Aghahosseini and C. Breyer, “From hot rock to useful energy: a global estimate of enhanced geothermal systems potential,” Applied Energy, vol. 279, Article ID 115769, 2020.

Y. Wang, Y. Liu, J. Duo, M. Liu, and M. Zeng, “Geothermal energy in China: status, challenges, and policy recommendations,” Utilities Policy, vol. 64, Article ID 101020, 2013.

Z. Feng, Y. Zhao, A. Zhou, and N. Zhang, “Development program of hot dry rock geothermal resource in the Yangbajing Basin of China,” Renewable Energy, vol. 39, no. 1, pp. 490–495, 2012.

R. M. Potter, E. S. Robinson, and M. C. Smith, Method of Extracting Heat from Dry Geothermal Reservoirs, US, Los Alamos, NM, USA, 1974.

E. L. Majer, R. Baria, M. Stark et al., “Induced seismicity associated with enhanced geothermal systems,” Geothermics, vol. 36, no. 3, pp. 185–222, 2007.

M. D. Zoback and H.-P. Harjes, “Injection-induced earthquakes and crustal stress at 9 km depth at the KTB deep drilling site, Germany,” Journal of Geophysical Research: Solid Earth, vol. 102, no. B8, pp. 18477–18491, 1997.

R. P. Denlinger and D. R. H. O’Connell, “Evolution of faulting induced by deep fluid injection, Paradox Valley, Colorado,” Bulletin of the Seismological Society of America, vol. 110, no. 5, pp. 2308–2327, 2020.

R. Baria, S. Michelet, J. Baumgärtner et al., “Creation and mapping of 5000 m deep HDR/HFR reservoir to produce electricity,” in Proceedings World Geothermal Congress 2005, vol. 8, p. 1627, Antalya, Turkey, April 2005.

K. H. Palgunadi, A. A. Gabriel, T. Ulrich, J.-N. López-Comino, and P. M. Mai, “Dynamic fault interaction during a fluid-injection-induced earthquake: the 2017 5.5 Pohang event,” Bulletin of the Seismological Society of America, vol. 110, no. 5, pp. 2328–2349, 2020.

C. J. Bromley, C. F. Pearson, and D. M. Rigor, “PNOC-EDC, microearthquakes at the Puhagan geothermal field, Philippines—a case of induced seismicity,” Journal of Volcanology & Geothermal Research, vol. 31, no. 3-4, pp. 293–311, 1987.

N. Deichmann and J. Ernst, “Earthquake focal mechanisms of the induced seismicity in 2006 and 2007 below Basel (Switzerland),” Swiss Journal of Geosciences, vol. 102, no. 3, pp. 457–466, 2009.

S. Sherburn, C. Bromley, S. Bannister, S. Sewell, and S. Bourguignon, “New Zealand geothermal induced seismicity: an overview,” in Proceedings of the World Geothermal Congress, Melbourne, Australia, April 2015.

K. y. Nishigami, K. Tadokoro, S. Nagai, T. Mizuno, Y. Kano, and Y. Hiramatsu, “Triggering characteristics of induced earthquakes: brief review of water injection experiments in the Nojima fault and other regions,” Journal of Geography (Chigaku Zasshi), vol. 111, no. 2, pp. 268–276, 2002.

F. Grigoli, S. Cesca, A. P. Rinaldi et al., “The November 2017 5.5 Pohang earthquake: a possible case of induced seismicity in South Korea,” Science, vol. 360, no. 6392, pp. 1003–1006, 2018.

M. D. Zoback and S. M. Gorelick, “Earthquake triggering and large-scale geologic storage of carbon dioxide,” Proceedings of the National Academy of Sciences, vol. 109, no. 26, pp. 10164–10168, 2012.

J. E. Streit and R. R. Hillis, “Estimating fault stability and sustainable fluid pressures for underground storage of CO₂ in porous rock,” Energy, vol. 29, no. 9-10, pp. 1445–1456, 2004.

J. Rutqvist, J. Birkholzer, F. Cappa, and C.-F. Tsang, “Estimating maximum sustainable injection pressure during geological sequestration of CO₂ using coupled fluid flow and geomechanical fault-slip analysis,” Energy Conversion and Management, vol. 48, no. 6, pp. 1798–1807, 2007.

X. C. Wei, Q. Li, X.-Y. Li, Y.-K. Sun, and X. H. Liu, “Uncertainty analysis of impact indicators for the integrity of combined caprock during CO₂ geosequestration,” Engineering Geology, vol. 196, pp. 37–46, 2015.

F. Cappa and J. Rutqvist, “Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO₂,” International Journal of Greenhouse Gas Control, vol. 5, no. 2, pp. 336–346, 2011.

J. Rutqvist, A. P. Rinaldi, F. Cappa, and G. J. Moridis, “Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs,” Journal of Petroleum Science and Engineering, vol. 107, pp. 31–44, 2013.

T. H. W. Goebel and E. E. Brodsky, “The spatial footprint of injection wells in a global compilation of induced earthquake sequences,” Science, vol. 361, no. 6405, pp. 899–904, 2018.

K. M. Keranen, M. Weingarten, G. A. Abers, B. A. Bekins, and S. Ge, “Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection,” Science, vol. 345, no. 6195, pp. 448–451, 2014.

H. Liu, H. Wang, H. Lei, L. Zhang, M. Bai, and L. Zhou, “Numerical modeling of thermal breakthrough induced by geothermal production in fractured granite,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 12, no. 4, pp. 900–916, 2020.

D. Zbinden, A. P. Rinaldi, T. Diehl, and S. Wiemer, “Hydromechanical modeling of fault reactivation in the St.Gallen deep geothermal project (Switzerland): poroelasticity or hydraulic connection?” Geophysical Research Letters, vol. 47, pp. 1–10, 2020.

Hibbitt, Karlsson, and Sorensen, ABAQUS Theory Manual, Version 6.2, Hibbitt, Karlsson, Sorensen Inc., Providence, RI, USA, 2001.

M. A. Biot, “General theory of three-dimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941.

J. R. Rice and M. P. Cleary, “Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents,” Reviews of Geophysics, vol. 14, no. 2, pp. 227–241, 1976.

K. Aki, “Scaling law of seismic spectrum,” Journal of Geophysical Research, vol. 72, no. 4, pp. 1217–1231, 1967.

H. Kanamori and D. L. Anderson, “Theoretical basis of some empirical relations in seismology,” Bulletin of the Seismological Society of America, vol. 65, no. 5, pp. 1073–1095, 1975.

H. Kanamori, “The energy release in great earthquakes,” Journal of Geophysical Research, vol. 82, no. 20, 1977.

T. C. Hanks and H. Kanamori, “A moment magnitude scale,” Journal of Geophysical Research, vol. 84, no. B5, pp. 2348–2350, 1979.