Study of Characteristics of Fault Slip and Induced Seismicity during Hydraulic Fracturing in HDR Geothermal Exploitation in the Yishu Fault Zone in China

Hot dry rock (HDR) geothermal energy has become promising resources for relieving the energy crisis and global warming. The exploitation of HDR geothermal energy usually needs an enhanced geothermal system (EGS) with artificial fracture networks by hydraulic fracturing. Fault reactivation and seismicity induced by hydraulic fracturing raise a great challenge. In this paper, we investigated the characteristics of fault slip and seismicity by numerical simulation. The study was based on a hydraulic fracturing project in the geothermal field of Yishu fault zone in China. It revealed that fluid injection during hydraulic fracturing can cause the faults that exist beyond the fluid-pressurized region to slip and can even induce large seismic event. It was easier to cause felt earthquakes when hydraulic fracturing was carried out in different layers simultaneously. We also examined the effects of the location, permeability, and area of the fracturing region on fault slip and magnitude of the resulting events. The results of the study can provide some useful references for establishing HDR EGS in Yishu fault zone.


Introduction
With the global economic development, the increasing depletion of fossil resources such as petroleum, coal, and natural gas is very serious, accompanied by increasingly serious environmental problems such as environmental pollution and global warming. The energy crisis and global warming have become two major problems of the world today. Accelerating the exploitation of clean and renewable energy has become the goal of many countries around the world. Geothermal resources are renewable and low carbon, becoming promising energy resources for relieving the energy crisis and global warming. Among them, geothermal energy from hot dry rock (HDR), being regarded as an almost inexhaustible source of clean energy, has the biggest development potential. HDR resources refer to the rock mass with temperature above 150°C buried 3-10 km underground, having no or only a small amount of water [1]. The total heat energy within HDR worldwide is very huge. Just the HDR geothermal energy within the most easily accessible zone (less than 10 km deep) is up to 40-400 M quads (1 quad = 1:0551 × 10 18 J), which is 100-1000 times than the available energy from all total fossil fuels worldwide [2]. HDR geothermal energy can be developed for heating generation which can be directly applied in domestic heating, aquaculture and agricultural heating, business and tourism heating (e.g., bathing and swimming), and other industrial uses. It can also be used for power generation. In addition to be clean and rich reserves, HDR geothermal is widely distributed and independent from variations in weather. Therefore, full exploitation of HDR heat resource plays an important role in energy supply and environmental protection for all countries in the world.
Cold fluid can be injected via injection wells into hot dry rocks and travels through them to capture the heat. The heated fluid then can be extracted to the surface via production to obtain heat energy. However, the natural HDR formations usually consist of hard and low-permeability granite.
Fluid cannot flow in them. Therefore, it is necessary to create artificial fracture networks in the HDR formations to improve their permeability to help the exploitation of the heat energy, which is called the enhanced geothermal system (EGS). The concept of EGS was first proposed in USA by Potter et al. (1974) as a patent in 1974 [3]. Hydraulic fracturing is one of the most significant technologies for establishing EGS in hot dry rocks [4]. Firstly, a well is drilled into the target region (Figure 1(a)). Secondly, high-pressure fracturing fluid (usually water) is pumped into rocks via the well to form cracks or open the natural joints/fractures in the rocks (Figure 1(b)). Thirdly, fracturing proppant (e.g., quartz sand) is injected into the fractures to avoid them reclosing and achieve sufficient and stable flow paths (Figure 1(c)) [5,6]. Finally, cold fluid is injected via injection well and heated by circulation through the fractured region and then extracted via production well to the surface to deliver the captured heat for power generation or other uses (Figure 1(d)).
In the course of underground fluid injection, fault reactivation and induced seismicity often raise a great challenge for its development. It is usually considered as a result of the effective stress principle [7], where the increase in pore pressure reduces the effective stress on preexisting fault which may lead to shear failure and seismic events [8]. This means that the induced seismicity by fluid injection is regarded as a fluid structure interaction (FSI) problem. Only when the injected fluid flows into or through the fault zone (as fault I shown in Figure 2) can reduce its effective stress and induce seismicity [9,10]. Influenced by the "principle of effective stress", the fluid-induced seismicity is usually interpreted as the overpore pressure reaching faults and triggering them slip [11]. Most previous studies have been based on this interpretation, no matter theory study [12][13][14], lab experimental study [15], or numerical simulation study [16][17][18]. Earlier still, in 1959, Hubbert and Rubey concluded the over pore pressure would reduce the "effective strength of rock" and thus weaken the preexisting fault [19]. In 1978, M. Lee Bell and Amos Nur suggested the "rapid unloading may cause instantaneous crustal weakening owing to excess pore pressure" [20]. Over a long period in the past, the stability of faults that are not in the pore-fluid migration region (as fault II shown in Figure 2) has been neglected to a certain extent. However, more and more researches suggest that fluid injection can also induce fault slip [21] or even trigger seismicity [22,23] outpaces pore-fluid migration. The overpressure decreased the effective stress of the injection region, leading it to an unloading state. The stress field of the region is changed, with the balance of initial stress in the surrounding rock mass being destroyed. Elasticity of rocks is an effective means    2 Geofluids of transmitting stresses to great distances [22]. The perturbation of stress field can spread to the faults and may cause them slip. Therefore, fluid injection induced fault slip, and seismicity outpaces pore-fluid migration which is a very important factor to be considered in a practical geothermal engineering. Induced seismicity has also become a problem that must be considered in EGS which may cause fear and panic among local residents, especially after the November 2017 M w 5.5 Pohang earthquake (South Korea) that has injured about 70 people and caused extensive damage in and around the city [24]. It will happen not only in the process of hydraulic circulation for heat recovery (as shown in Figure 1(d)) but also in the process of hydraulic fracturing (as shown in Figure 1(c)) [25].
China has a huge amount of HDR geothermal resources, about 23.69 M quads which is equivalent to 8.56 trillion tons of standard coal. It has become a strategic choice to ensure energy supply and social sustainable development [26]. In September 2017, dry hot rocks with temperature up to 236°C were drilled at 3705 m depth in Qinghai Gonghe Basin of China [25]. It is a major breakthrough in the exploration of dry hot rock, which directly confirmed that China has high quality dry hot rock resources. Favorable HDR resource development zones of China include Southern Tibet, Western Yunnan (Tengchong), Southeast Coast (Zhejiang, Fujian and Guangdong Provinces), North China (Bohai Bay Basin), on the southeast margin of Ordos Basin (Fenwei graben), and Northeast China (Songliao Basin), [27]. The east of the Yishu fault zone (Figure 3) is one of the best geothermal areas in Shandong Province [28]. In July 2019, a large number of dry hot rock rich areas (Ju County, Wulian County of Rizhao City and Wendeng District of Weihai City) have been found, having resources equivalent to 18.779 billion tons of standard coal. To establish the EGS in the Yishu fault zone and reduce the damage caused by induced seismicity, fault slip and induced seismicity characteristics during HDR geothermal exploitation in the Yishu fault zone were investigated in this paper by numerical simulation. Numerical simulation is an effective tool for solving problems of geotechnical engineering [29][30][31]. We used the finite-element software, ABAQUS, to estimate the induced fault slippage during hydraulic

Geofluids
fracturing. ABAQUS has advantages in solving the inelastic, nonlinear, anisotropic, and heterogeneous problems.
Study on fault stability is vital for improving seismic hazard assessment and mitigation. In this paper, the study was focused on fault slip induced by fluid injection during the process of hydraulic fracturing (Figure 1(c)). In addition, we focused on the faults that exist beyond the hydraulic fracturing region. Considering a hydraulic fracturing case, we investigated the effects of fluid injection on stability of the selected faults through the following several aspects: (1) the influence of relative position between fracturing region and fault on induced fault ship and seismicity, (2) the influence of fracturing region's permeability on stability of the fault, (3) the influence of fluid injection pressure during hydraulic fracturing on stability of the fault, and (4) the influence of the area of the fracturing region on stability of the fault.

Numerical Method
2.1. Interaction between Fault and Its Surrounding Rock. The fault and its surrounding rock mass were treated as two parts (Figure 4(a)). Estimation of the induced fault slippage should well solve the interaction between the fault and surrounding rock. Contact analysis can be used to simulate the interaction between them. Contact is a special discontinuous constraint with nonlinear behavior of severe discontinuity. ABAQUS was selected to solve these problems due to its great advantages in dealing with nonlinear calculation [32].
2.1.1. Defining Contact Pair. Contact pair in ABAQUS was used to complete the contact analysis between the two deformable surfaces (fault and its surrounding rock). In the contact pair, force can be transferred through their interface. Normal force is transmitted along the vertical direction. When the contact pressure becomes zero or negative, the two surfaces will separate. Shear force can be transmitted along the tangential direction when there is friction between the two surfaces. The master-slave contact method was used to define the contact pair of the fault and surrounding rock. Before the definition of contact pair, it was needed to deter-mine which one was the master surface and which one was the slave. According to the guidelines of "defining a surface-based contact simulation" in ABAQUS [32], the surface of THE surrounding rock was chosen as the master surface ( Figure 4(b)), and the surface of fault was chosen as THE slave surface (Figure 4(c)). The reason was that the stiffness of the surrounding rock mass was larger than that of fault [30].

Defining Contact Properties.
After the definition of contact pair, the properties of contact were needed to be defined. To avoid one surface penetrating to another, the "hard" contact was assigned in the normal direction of the interface of contact. The Coulomb friction model was used in the tangential direction. In this model, the two contacting surfaces can carry shear stresses up to a certain magnitude across their interface before they start sliding relative to one another. The Coulomb friction model defines this critical shear stress, τ cirt , at which sliding of the surfaces starts as a fraction of the contact pressure, p, between the surfaces (τ cirt = μp, μ is the coefficient of friction) [32].

Contact Tracking
Approach. It has two tracking approaches to account for the relative motion of two contacting surfaces in mechanical contact simulations: finite sliding and small sliding. In small-sliding formulation, only relatively small sliding relative to each other of the interacting surfaces is allowed. In the finite sliding formulation, however, any size of sliding is permitted. Small-sliding contact is computationally less expensive than finite-sliding contact. However, the degree of the induced fault slip along the surrounding rock often cannot predict in advance. Large slippage may also occur. Therefore, the finite sliding formulation was used for this study.
The contact is described by a node n 1 on the slave surface with a quadratic segment with nodes n 2 , n 3 , and n 4 ( Figure 5). To derive the equations governing these elements, the coordinates in the plane of the slide line is considered. Firstly, the point x on the segment closest to the point x 1 on the slave surface is determined. Secondly, the normal n and tangent t to the segment at that point is determined.   4 Geofluids Then, point x and the normal n have the following relationship [32]: where h is the overclosure. The slip can be solved by the following equations: where s is the slip. Detailed derivation of the equations can refer to ABA-QUS theory manual [32].

Hydromechanical Coupling of the Hydraulic Fracturing
Region. The process of fluid injection into the fracturing region was regarded as a problem of fluid-structure coupling. In this study, the simulation of the hydraulic fracturing region is based on porous elastic theory. The porous elastic medium is described by the flowing equations [33,34]: where G is the shear modulus, ε ij is the strain, σ ij is the stress, σ kk is the mean stress deviation, Pa, p is the pore pressure, Pa, v u and v are undrained and drained Poisson's ratios, respectively, δ ij is the Kronecker delta, Δm is the change in fluid mass per unit volume, kg/m 3 , ρ 0 is the density of the pore fluid in the reference state, kg/m 3 , B is Skempton's B coefficient, K and K s are drained bulk modulus and solid grain's bulk modulus, and 1/K n and K f are the unjacketed pore compressibility and pore fluid bulk modulus, respectively.
The flow is described by Darcy's law and a mass conservation equation. Darcy's law describes the pore fluid flowing in a porous medium: where f is the volumetric flux vector, m/s, κ is the intrinsic permeability of the medium, m 2 , μ is the viscosity, Pa·s, and ∇p is the pressure gradient vector, Pa/m. The mass conservation equation of pore fluid is [35]: where K is the fluid mass.  Aki (1967) can be used to quantify the released energy, as follows [36]: where μ is the rigidity of fault, Pa, A is the fault rupture area and A = lw (l and w are length and width of fault rupture, respectively), and d is the average slip of the rupture area.
On the base of seismic moment M 0 , moment magnitude scale M S proposed by Kanamori and Anderson (1975) can be calculated to measure the strength of the induced seismic event [37][38][39], as follows:

Conceptual Model.
A two-dimension model was designed to simulate the EGS-Fault system by ABAQUS [32]. A numerical model was simplified from a hydraulic fracturing field test in Yishu fault zone. We mainly investigated the characteristics of the fault slippage and seismicity induced by fluid injection during hydraulic fracturing of EGS engineering. According to geological exploration data, two major faults (named F1 and F2) that adjacent to but beyond the hydraulic fracturing region were selected. Contact analysis was used to describe the interaction behavior between the faults and the surrounding rock. Mohr-Coulomb failure criterion was used to control the slip of fault. A plane strain model with dimensions of 8,000 m × 400 m ( Figure 6) was developed to simulate the geothermal system. The thicknesses of caprock and granite were 1000 m and 3000 m, respectively. Actually, the total caprock above granite was about 2000 m. In order to improve the simulation efficiency, the caprock in our model was only 1000 m. The other 1000 m caprock was replaced by the uniformly distributed

Geofluids
pressure that equals to its weight. The length of fault F1 was 2000 m, with the length in granite and caprock that was 1500 m and 500 m, respectively. The length of fault F2 is 1000 m, with the length in granite and caprock that is 500 m and 500 m, respectively. The thicknesses and inclinations of fault F1 and F2 were 10 m and 60°, respectively. In the field test, hydraulic fracturing was carried out in three different layers. As shown in Figure 6, three elliptic regions (regions N1, N2, and N3) were used to simulate the fracturing region in each layer. The long axis and short axis of the ellipse were 600 m and 200 m, respectively. The vertical distances from the center of regions N1, N2, and N3 to bottom of fault F1 were 0 m, 300 m and 600 m, respectively. The distance from the midpoint of the fracturing region N1 to the ground surface was 3600 m. The horizontal distance between the center of region N1 and the bottom of fault F1 was 1220 m. Injection well was a vertical well along the short axis of the ellipse region. 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. The parameters that could not be measured were mainly determined from other public literature [40,41].
The boundary conditions of model are (1) the bottom of the model was fixed, in which both vertical and horizontal displacement are not allowed, (2) the left and right sides of the model are constrained by rollers, only vertical displace-ment was allowed, and (3) the top of the model is allowed to deform freely.
The stress within the fault-stratum system should be equal to the lithostatic stress prior to hydraulic fracturing for geothermal exploitation. The initial stresses should be geostatic. Water pressure is hydrostatic within the hydraulic fracture region, with the initial value to be 0 MPa. Firstly, the gravity is considered as external force applied to the model. After the application of gravity, the stresses of each node are obtained. Secondly, the stresses are inputted to the model as initial internal force before gravity to reach a mechanical equilibrium between the initial stress fields and the gravity under the applied boundary conditions of the model. Then, fluid is injected into the ellipse region to simulate the process of hydraulic fracturing.

Numerical Experiments.
Fault reactivation and slip are a complex mechanical behavior of fault and its surrounding rocks involving the interaction of total stresses, pore pressures, and temperature stresses [42]. In this study, we focused on the additional influence of pore pressure changes within the fracturing region on induced fault slippage and seismicity.  fracture is very weak. The injected fluid can only flow in the fracturing region. However, the fluid pressurization in the fracturing region can transmit solid stresses to the adjacent faults that may trigger seismic event. The location of the fracturing region will directly control the stress perturbation of fault zone and consequently affect its stability. Therefore, the location of the fracturing region is a very important engineering parameter that influences the stability of faults. Regions N1, N2, and N3 were the three actual fracturing regions in the field test of hydraulic fracturing. In order to fully study the influence of the location of the fracturing region on fault's stability, we designed three hypothetical fracturing regions (regions P1, P2, and P3) that were closer to fault F1 (Figure 7). Regions P1, P2, and P3 can be obtained by shifting N1, N2, and N3 to the left with 820 m, horizontally. The horizontal distance between the left endpoint of region P1 and the bottom of fault F1 was 100 m. The distance between region P1 and N1 was 220 m. In summary, a total of six hydraulic fracturing regions were designed. As shown in Table 2 In these numerical experiments, all injection wells are vertical with the injection pressure of 80 MPa. After 120 hours of fluid injection, we obtained the slippage along the fault F1 and F2 (Figures 9-14). Based on them, the seismic moment magnitude M w can be then calculated (Tables 3-5). 3.2.2. Permeability of the Fracturing Region. The engineering effect of an enhanced geothermal system is largely dependent on the permeability of the fracturing region. Too low or too high permeability will both reduce the thermal recovery effect. The permeability directly affects fluid flowing in the fracturing region. If the permeability is too low, the fluid flows very slowly which affects the heat recovery rate. If the permeability is too high, however, the fluid flows very fast which affects the heat absorption from HDR. The permeability will also significantly affect the pressure build-up in the fracturing region and then affect the stress perturbation of the adjacent faults. When the injection pressure and time are kept constant, the fracturing region permeability will become the key factor that affects the stability of fault. Therefore, the influence of fracturing region permeability on fault slip and induced seismicity needs to be investigated. Based on experiments No. C1 (Figure 8(a)), a total of 5 numerical experiments (Nos. C1K1-C1K5) of this section were designed by changing the fracturing region permeability from 10 -10 m/s to 10 -6 m/s ( Table 6).
In these numerical experiments, all injection wells are vertical with the injection pressure of 80 MPa. After 120 hours of hydraulic fracturing by fluid injection, we obtained the slippage along the faults F1 and F2 (Figures 15 and 16). Based on fault's slippage, the seismic moment magnitude can be then calculated (Table 7).

Injection Pressure.
Injection pressure is another important factor that affects stress perturbation in the fracturing region and adjacent rocks. The internal pore pressure in the region increases rapidly when fluid is injected, reducing the effective stress of the region and leading it to an unloading state. The region will expand and cause the surrounding rock to lose stress balance. The stress perturbation may transfer to the nearby faults, which may induce them to loss stability. The influence of injection pressure on the induced fault slip and seismicity needs to be investigated. Based on the experiment No. C1, a total of 5 numerical experiments (Nos. C1P1-C1P5) of this section were designed by changing the injection pressure from 60 MPa to 70 MPa, 80 MPa, 90 MPa, and 100 MPa (Table 8).
In these numerical experiments, the permeability of the fracturing region is 10 -10 m/s. All injection wells are vertical and after 120 hours of hydraulic fracturing by fluid injection, we obtained the slippage along the faults F1 and F2 (Figures 17 and 18). Based on fault's slippage, the seismic moment magnitude can be then calculated as shown in Table 9.    Table 10). In these experiments, the fracturing region permeability and injection pressure were, respectively, remained constant at 10 -8 m/s and 80 MPa. After 120 hours of fluid injection, we obtained the slippage along the faults F1 and F2 (Figures 19 and 20). The fault slippage and induced seismic moment magnitude M w were as shown in Table 9.    Pa. What we used in this study was a twodimension model, and the length of fault plane F1 was 2000 m, but the width was unknown. According to the relevant observations and geological data, the width of the fault plane was assumed to be 200 m. Therefore, the rupture area A of F1 was 400,000 m 2 . According to the above data and equation (10) Table 3. The results were clearly inconsistent with our intuitive judgment. It was generally though that the influence of the hydraulic fracturing region on fault was larger, when it was closer to the fault. However, the results showed that the injection-induced fault slip became larger when the fracturing region changed downward. From the table, we can see that with the fracturing region changing   9 Geofluids from N1 to N2 and then to N3, the induced slippage and moment magnitude of fault F1 both became larger. In fact, the influence of fluid injection on the adjacent faults' stability is not simply depended on the position of injection region, but on the stress disturbance region induced by it. It is the disturbance of stress and displacement field result from fluid injecting in the fracturing region that leads the fault slip. Figure 21 illustrated the distribution of vertical displacement of experiments Nos. C1 (Figure 21(a)), C2 (Figure 21(b)), and C3 (Figure 21(c)). When fluid was injected into the fracturing region, the region expanded, causing the upper surrounding rock to uplift and the lower surrounding rock to subside. The disturbance degree of the displacement field close to fault can directly affect its slid. When injection pressure and injection time remained constant, the displacement disturbance field close to the fracturing region kept the same. With the fracturing region changing from N1 to N2 and then to N3, the displacement disturbance field moved down, correspondingly. It may lead the displacement disturbance of fault F1 and its near area become much larger. Taking the displacement contour of U2 ≈ 0:1 m for example, it was approximate ellipse in the three experiments (as indicated by black dotted lines in Figure 21). When fluid was injected into regions N1, N2, and N3, the minimum distance between the contour and fault F1 was 780 m, 700 m, and 640 m, respectively. It meant that the displacement disturbance field     close to fault F1 became larger. This can explain why the moment magnitude M w became larger when the fracturing region became far away the fault F1. Since the displacement contour of U2 ≈ 0:1 m was roughly elliptical, we could foresee that the distance between the contour and fault F1 would become smaller if the region continues to move down. Then, the moment magnitude M w would not become larger. In contrast, it would become smaller. In order to verify this, we designed another numerical experiment by only moving region N3 down by 600 m to be N4 (the vertical distance from the center of N4 to the bottom of fault was then 1200 m). After 120 hours of fluid injection, we obtained the vertical displacement field and fault slippage as shown in Figures 22 and 23. When the fracturing region was N4, the maximal and average slippage of fault F1 was 8.21 mm (occurred at 1500 m of the fault plane) and 5.15 mm, respectively. The moment magnitude M w was 2.6732, and it was smaller than that when fluid was injected into regions N1, N2, or N3.

Results and Discussion
Secondly, we studied the induced slippage along fault F2 ( Figure 10). When fluid was injected into region N1, fault F2 slipped upward along the fault plane as shown in the black line with square symbol. The slippage along the upper half of the fault was obviously larger than that along the lower half. From 0 to 500 m, the maximal and average slippage was 4.70 mm and 3.80 mm, respectively. From 500 m to 1000 m, however, the maximal and average slippage was 4.70 mm and only 0.90 mm, respectively. As shown in the displacement cloud of No. C1 in Figure 21(a), the fluidinjection-induced displacement near the upper portion of the fault was larger than that near the lower portion. Along the whole fault plane of F2, the maximal slippage fault was 4.70 mm (occurred at 160 m) with the average slippage of 2.29 mm. The moment magnitude M w of the induced seismicity of fault F2 was 2.2378. When fluid was injected into region N2, the displacement disturbance field moved down along with the fracturing region (as shown in Figure 21(b)). This caused the slippage along the upper portion of the fault to become smaller and that along lower portion of the fault to become larger. As shown in the red line with circle symbol in Figure 10, the difference of the slippage along the upper and lower part of the fault became smaller. From 0 to 500 m, the maximal and average slippage was 2.84 mm and 2.52 mm, respectively. From 500 m to 1000 m, the maximal and average slippage was 1.96 mm and only 1.67 mm, respectively. Along the whole fault plane of F2, the maximal and average slippage of fault F2 was 2.84 mm (occurred at 250 m of the fault plane) and 2.14 mm, respectively. The moment magnitude M w of the induced seismicity of fault F2 was 2.2182. When fluid was injected into region N3, the induced slippage of fault F2 displayed similar tendency. The maximal and average slippage of fault F2 was 2.69 mm (occurred at 250 m of the fault plane) and 2.10 mm, respectively. The moment magnitude M w of the induced seismicity of fault F2 was 2.2128. From Table 3, we can see that when the fracturing region changing from N1, N2, to N3, the induced slippage and moment magnitude M w of fault F2 both became smaller. It was contrary to that of fault F1. That was because when the fracturing region changing from N1, N2, to N3, the displacement disturbance close to fault F2 became smaller.
The above results revealed that the influence of the hydraulic fracturing region on stability of fault was not simply depended on their relative spatial location. It was the displacement and stress disturbance induced by fluid injection in the fracturing region that ultimately caused the fault sliding instability.
(2) Hypothetical Fracturing Program. Figure 11 shows the fault slippage of fault F1 of experiment Nos. H1-H3. The induced slippage of fault F1 in the three experiments displayed similar tendency. The slippage along the upper portion of the fault was obviously smaller than that along the lower portion. Because the hypothetical fracturing region was near the bottom of fault F1, the disturbance degree of stress and displacement in the lower portion of fault F1 and its adjacent area was larger. Taking experiment No. H1 for example, as shown in Figure 24, the vertical displacement of the lower portion fault F1 and its adjacent area was obvious larger than that in upper portion. When fluid was injected into region P1 (the black line with square symbol in the Figure) 11 Geofluids fracturing region changed from P1 to P2 and then to P3, the induced slippage and moment magnitude M w of fault F1 became smaller. It was different to that when fluid was injected into regions N1, N2, and N3, respectively. This can be also explained by the vertical displacement cloud as shown in Figure 24. With the fracturing region changing from N1 to N2 and N3, the range and degree of displacement disturbance in the lower portion of fault F1 and its adjacent area became smaller. (2) If the fracturing region was too close to the fault, as the region P1, injecting fluid into this region may induce felt seismic event (M w ≥ 3). Therefore, determination of rational position of fracturing region in HDR EGS engineering is a key step for site stability.  Table 4, it can be seen that the maximal slippages of fault F2 in Experiment Nos. H1-H3 were 2.35 mm, 1.64 mm, and 1.01 mm, respectively. The average slippage was 1.00 mm, 0.78 mm, and 0.62 mm, respectively. Accordingly, the moment magnitude M w of the induced seismicity of induced seismicity was 1.99, 1.9260, and 1.8595, respectively. Because the fracturing regions P1-P3 were far from fault F2, the induced fault slippage and moment magnitude M w were very small. Figure 13 illustrates the slippage of fault F1 when fluid was simultaneously injected into three     Figure 14 illustrates the induced fault slippage of F2 of experiments No. C4 and No. H4. Fault F2 slipped downward along the whole fault plane in experiment No. C4. However, in experiment No. H4, it presented two sections with slipping upward along the upper portion and slipping downward along the lower portion. When fluid was injected into regions N1, N2, and N3 at the same time, the maximal slippage was 10.08 mm (occurred at 230 m), and the average slippage was 8.50 mm. In this case, the moment magnitude M w was 2.6176. When fluid was injected into regions P1, P2, and P3 at the same time, the maximal slippage along the whole fault plane was only 2.80 mm, and the average slippage was 1.40 mm. The moment magnitude M w was only 2.0954.   Figure 19: Pore pressure profiles along the long axis of region N1 (PP′) of different injection pressures (Nos. C1P1-C1P5).

Geofluids
Based on the results, the induced slippage and moment magnitude M w were obtained as shown in Table 5. Comparing it with Tables 3 and 4, it can be seen that the induced fault slippage and moment magnitude M w of experiment No. C4 were much larger than that of experiment Nos. C1-C3. It revealed that it was easier to cause felt earthquakes when hydraulic fracturing was carried out in different layers simultaneously. In the current hydraulic fracturing program, three different fracturing layers between fault F1 and F2 were set. It must be treated with caution when decided whether to carry out hydraulic fracturing in different layers at the same time. Figures 15 and 16 and Table 7 show the results of the induced fault slippage and moment magnitude M w of experiment Nos. C1K1-C1K5. When permeability of the fracturing region changed from 10 -10 m/s to 10 -6 m/s, the fault slippage and moment magnitude of fault F1 and F2 both showed an increasing tendency.

Permeability of the Fracturing Region.
Firstly, we studied the induced fault slippage of F1 as shown in Figure 9 and Table 7. When the permeability was 10 -10 m/s (the black curve with the symbol of square), the maximal and average-induced slippages were 1.20 mm and 0.93 mm, respectively. The moment magnitude M w was 2.1779. When the permeability was 10 -9 m/s (the red curve with the symbol of circle), the maximal and averageinduced slippages were 3.35 mm and 2.60 mm, respectively. The moment magnitude M w was 2.1779. When the permeability was 10 -8 m/s (the blue curve with the symbol of triangle), the maximal and average-induced slippages were 6.93 mm and 5.34 mm, respectively. The moment magnitude M w was 2.6837. When the permeability was 10 -7 m/s (the green curve with the symbol of pentagram), the maximal and average-induced slippages were 7.04 mm and 5.43 mm, respectively.    14 Geofluids the permeability was 10 -6 m/s (the magenta curve with the symbol of regular pentagon), the maximal and averageinduced slippages were 7.24 mm and 5.58 mm, respectively. The moment magnitude M w was 2.6964. Secondly, we studied the induced fault slippage of F2 as shown in Figure 16 and Table 7. When the permeability was 10 -10 m/s, the maximal and average-induced slippages were 0.84 mm and 0.40 mm, respectively, with the moment magnitude M w of the induced seismicity of 1.7302. When the permeability was 10 -9 m/s, the maximal and averageinduced slippages were 2.32 mm and 1.11 mm, respectively, with the moment magnitude M w of the induced seismicity of 2.0286. When the permeability was 10 -8 m/s, the maximal and average-induced slippages were 4.70 mm and 2.29 mm, respectively, with the moment magnitude M w of the induced seismicity of 2.2380. When the permeability was 10 -7 m/s, the maximal and average-induced slippages were 4.77 mm and 2.33 mm, respectively, with the moment magnitude M w of the induced seismicity of 2.2427. When the permeability was 10 -6 m/s, the maximal and averageinduced slippages were 4.91 mm and 2.39 mm, respectively,  with the moment magnitude M w of the induced seismicity of 2.2506.
The above results revealed that the slip of fault F1 and F2 obviously increased with the increase of fracturing region permeability. When the permeability increased from 10 -10 m/s to 10 -8 m/s, the moment magnitude M w increased by about 0.5058. However, it had almost no increase when the permeability increased from 10 -8 m/s to 10 -6 m/s. The magnitude increased by only 0.0127. Figure 25 illustrates the distribution of pore pressure along the long axis of the fracturing region N1 of experiment Nos. C1K1-C1K5. The pore pressure was distributed symmetrically along the injection well (short axis of region N1) and decreased from the well to both sides. When the permeability was 10 -10 m/s, the black curve with square symbol in Figure 25, the pore pressure decreased rapidly from the injection well to both sides. The pore pressure was 80 MPa at the midpoint (at 300 m), and rapidly decreased to 7 MPa at 220 m and 380 m and then gradually decreased to 0.7 MPa at the endpoints. With the increase of fracturing region permeability, the pore pressure on both sides of the injection well obviously increased. When the permeability increased to 10 -8 m/s, the pore pressure at the endpoint was up to 76 MPa. When the permeability was 10 -7 m/s and 10 -6 m/s, the distribution of pore pressure was almost a horizontal line. It indicated that the whole fracturing region would be in high pore pressure state when its permeability was larger than 10 -7 m/s. Therefore, when the injection pressure and time kept constant, the permeability of the fracturing region became the key factor affecting the pressure build-up in the fracturing region and the stability of fault. On one hand, enhancing the permeability of granite through hydraulic fracturing is the engineering purpose for geothermal exploitation. On the other hand, too high a permeability will threaten the stability of the site fault. In addition, the heat transfer effect will be greatly affected if the fluid flows too fast in the fracturing region. Therefore, it is needed to appropriately control the permeability of the fracturing region for a successful geothermal exploitation project. Figures 17 and 18 and Table 9 show the results of the induced fault slippage and moment magnitude M w of experiments Nos. C1P1-C1P5. We can see that the induced slippage of fault F1 presented similar tendency under different injection pressures. Firstly, we studied the induced fault slippage of F1 as shown in Figure 18 and  Figure 18 and The above results revealed that the induced slippage and moment magnitude of fault F1 and F2 were all evenly increased with the injection pressure. However, the increment of fault F2 was smaller than that of F1. That was because the fault rupture area of F2 was smaller than that of F1, and F2 was farther from the fracturing region N1 than F1 was. Figure 19 illustrates the distribution of pore pressure along the major axis of region N1 under different injection pressures of experiments Nos. C1P1-C1P5. As can be seen, the pore pressure in region N1 increased with the injection pressure. The increment of pore pressure in the region was equal to the increment of injection pressure. This may explain why the induced fault slippage and moment magnitude M w evenly increased with the injection pressure. Table 11 show the results of the induced fault slippage and moment magnitude M w of experiments Nos. C1A1-C1A5. We can see that the induced slippage of fault F1 presented similar tendency under different areas of the fracturing region. Firstly, we studied the induced fault slippage of F1 as shown in Figure 20 and Table 11. As can be seen, when the area of the fracturing region was 0.36 A, 0.64 A, A, 1.44 A, and 1.96 A, the maximal fault slippage was 3.35 mm,  Figure 26 and

Conclusions
The exploitation of HDR geothermal energy usually needs an enhanced geothermal system with artificial fracture networks. Hydraulic fracturing is one of the most significant technologies for making artificial fracture networks in hot dry rocks. Fault reactivation and seismicity induced by fluid injection during hydraulic fracturing often raise a great challenge because it is easy to arouse social concern. A better understanding of injection-induced fault reactivation is important for improving disaster assessment and prevention of such underground fluid injection engineering. Whether the fault is connected with fluid region or not, underground fluid injection may induce fault reactivation. However, most studies have been focused on pore pressure increase of the fault hydraulically connected to the injection. The stability of faults that are not in the pore-fluid migration region has been neglected. In this paper, the induced fault slip and seismicity characteristics during HDR geothermal exploitation were studied by a 2D numerical model by ABAQUS. The study was based on a hydraulic fracturing project in the geothermal field of Yishu fault zone of China. In this study, we focused on the faults that exist beyond the hydraulic fracturing region. The results of the study can provide some useful references for establishing HDR EGS in the Yishu fault zone. The main conclusions were as follows: (1) Fluid injection during hydraulic fracturing can cause the fault that exist beyond the fluid-pressurized region to slip and can even induce large seismic event (M w ≥ 3 in our range of study) (2) The influence of the hydraulic fracturing region on stability of fault is not simply depended on their relative spatial location of the fault and fracturing region. It was the displacement and stress disturbance induced by fluid injection in the fracturing region that ultimately caused the fault sliding instability (3) If the fracturing region is close to the fault, as the region P1, injecting fluid into this region may induce felt seismic event. Therefore, determination of ratio-nal position of the fracturing region in an HDR EGS engineering is a key step for site stability in a hot dry rock geothermal engineering (4) It is easier to cause felt earthquakes when hydraulic fracturing is simultaneously carried out in different layers. In the current hydraulic fracturing program, three different fracturing layers between fault F1 and F2 were set. It must be treated with caution when decided whether to carry out hydraulic fracturing in different layers at the same time (5) When the injection pressure and time keep constant, the permeability of the fracturing region becomes the key factor affecting the pressure build-up in the fracturing region and the stability of fault. On one hand, enhancing the permeability of granite through hydraulic fracturing is the engineering purpose for geothermal exploitation. On the other hand, too high a permeability will threaten the stability of the site fault. In addition, the heat transfer effect will be greatly affected if the fluid flows too fast in the fracturing region. Therefore, it is needed to appropriately control the permeability of the fracturing region for a successful geothermal exploitation project (6) When the fracturing region and injection time keep constant, the increment of pore pressure in the fracturing region is equal to the increment of injection pressure. As a result, the induced fault slippage and moment magnitude evenly increase with the injection pressure (7) The induced slippage and moment magnitude of fault F1 and F2 both increase with the area of the fracturing region. In an HDR EGS engineering, the influence of the fracturing region area on the stability of the adjacent fault must be adequately assessed

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

Conflicts of Interest
The authors declare that they have no conflicts of interest.