Abstract
Shale gas hydraulic fracturing usually activates nearby faults and makes them slip. In horizontal wells, fault slip can result in serious casing deformation. Casing deformation slows the fracturing process, lowers production, and raises the cost of a shale gas well. It is challenging to obtain underground data on fault activation because of the deep shale reservoirs. As a result, the current study needed to indicate how hydraulic fracturing affects fault activation length. This made it challenging to control casing deformation. The fluid-structure coupled finite element method was used in this study to create a coupled seepage-stress model for heterogeneous shale formation. With microseismic signs and hydraulic fracking for shale gas, this model examined the variation law of pore pressure and ground stress. The fault activation coefficient was created to assess the fault activation duration and the impact of hydraulic fracturing. The model was verified by the microseismic signal. The outcomes of the numerical simulation demonstrate how the rapid rise in formation pore pressure during hydraulic fracturing affected the ground stress at the fault interface. The influence of ground stress variation at the fault interface on fault activation could not be ignored. Increased fault elastic modulus, fracture pressure, fracture time, and the fault Poisson ratio result in longer fault activation lengths. The length of the fault’s activation was decreased by the increase in fracture stage, distance from the fault, friction angle within the fault, and fault angle. Finally, a shale gas horizontal well with casing deformation in block C was analyzed. The results showed that reducing the fracturing duration can reduce the activation length of the fault by 68.25%, resulting in a 9.1 mm fault slide and a 0.86 mm casing deformation, respectively. This study offers theoretical guidelines for preventing fault activation during hydraulic fracturing in horizontal shale gas wells.
1. Introduction
Multistage hydraulic fracturing of horizontal wells is a critical technology in shale gas development. It is a widely used technology [1–3]. However, hydraulic fracturing also brings some problems for shale gas development. Fault activation and casing shear deformation occurred during shale gas hydraulic fracturing [4, 5]. The casing shear deformation hindered the running of the fracturing string, which seriously reduced the well construction efficiency of shale gas wells and increased the well construction cost [6]. The occurrence of casing shear deformation in shale gas horizontal wells was not an accidental phenomenon. It occurred in every shale gas block in Sichuan, China [7]. Casing shear deformation accounted for over 60% of the total casing deformation. It was the main form of casing deformation in shale gas horizontal wells [8]. In the current study, many scholars believed that the cause of casing shear deformation was fault slip [9–11]. Chen et al. [12] counted the casing deformation points in block C and found that casing shear deformation was related to fault slip. Li et al. [13] presented the relationship between casing shear deformation points and faults in two wells in block W, Sichuan, China. All casing shear deformation points were found to coincide with faults. We also counted 16 casing deformation points in the newly developed block L, Sichuan, China. Twelve of the casing deformation points coincided with faults or natural fractures. Casing deformation points correlated with faults up to 75%. Fault slip caused by hydraulic fracturing is the leading cause of casing shear deformation in horizontal shale gas wells. Therefore, evaluating fault activation during shale gas hydraulic fracturing was the basis for preventing fault slip and controlling casing shear deformation.
The current research reached a consensus on the fault slip mechanism. Many scholars believed that fracturing fluid entered the fault during multistage fracturing, which caused fault activation. Many scholars have done much valuable work on the mechanism of shale fault slip and the relationship between fault slip and casing deformation. Chen et al. [14] analyzed the mechanism of fault slip based on the focal mechanism principle and determined the relationship between fault slip and moment magnitude. Lu et al. [15] analyzed the mechanism of fault slip based on multiarm caliper logging data and gave the relationship between fault slip and casing deformation. Chen et al. [16] studied the relationship between fault (fracture zone) and casing deformation based on field data. The relationship between the fracture zone or fault and casing deformation was analyzed. The result showed that the fault (fracture zone) activation was the main controlling factor of casing shear deformation. Suling and Lei [17] analyzed the mechanism of casing deformation caused by fault slip based on the focal mechanism relationship. The results showed a positive correlation between fault slip and casing deformation. Meng et al. [18] established a numerical fault slip model and analyzed the casing shear deformation law under different fault dips. It showed that the casing deformation was maximum when the fault plane was perpendicular to the borehole trajectory. Li et al. [19] analyzed the casing deformation pattern and force when the casing was subjected to shear stress using laboratory experiments and proposed a new formula for calculating the casing strength. Lian et al. [20] studied the influence of different cement sheath properties and casing wall thickness on casing deformation. It was found that changing cement performance and casing wall thickness had little effect on casing deformation. Zoback et al. [21] and Daneshy [22] believed that asymmetric fracturing usually occurred in the fracturing process due to the heterogeneity of the formation, which resulted in shear dislocation. A calculation model of fault slip established by Kui et al. [23] was used to calculate the pore pressure variation in the formation during fracturing. Zhang et al. [24] established a pore pressure variation model during multistage fracturing to simulate the fault slip. Tan et al. [25, 26] developed a numerical model of hydraulic fracture extension using the XFEM-based CZM method and conducted large-scale fracturing experiments on anisotropic rocks. The fracture extension mechanism during hydraulic fracturing of deep shale was analyzed. Huang et al. [27] investigated the fracture extension mechanism of low-permeability reservoirs under rock anisotropy conditions based on a 2D discrete element method. Yoon et al. [28] developed a numerical model to study the relationship between high-rate fluid injection-induced earthquakes and fault systems. Taghipour et al. [29] used a numerical simulation method to evaluate the geometric mechanical properties of the CO2-injected formation and analyzed the mechanical state of the fault five years after CO2 injection. Vilarrasa et al. [30] proposed that the injection well should be far away from the fault to avoid fault activation. Zhang [31] developed a numerical model to analyze the effect of faults on casing damage and determine the maximum injection pressure allowed in the formation.
However, the above part of the study analyzed the mechanism of fault slip, and the other part only considered the influence of formation pore pressure changes on fault activation. Few studies considered the change of ground stress on the formation during hydraulic fracturing and the effect of this change on fault activation. The long-time and low displacement fluid injection affected by ground stress should be addressed. For example, CO2 injection into the stratum usually takes several years to decades, so the variation of ground stress by fluid injection can be neglected. However, shale gas hydraulic fracturing is a process in which large amounts of liquid are injected into a formation over a short period. This process causes the pore pressure of the formation to increase dramatically in a short time. At this time, the influence of pore pressure on ground stress cannot be ignored, and the change of ground stress can also impact fault activation. These changes also directly affect the fault activation length calculation and slip distance.
There are few studies on the effect of hydraulic fracturing on fault activation. Most of these studies are qualitative analyses and lack research on fault activation, formation characteristics, fracturing, and other construction parameters. The current research results can only provide quantitative guidance for fracturing construction. How to change the fracturing parameters to avoid fault activation caused by hydraulic fracturing is still a significant problem.
A coupled seepage-stress numerical model was established with the microseismic monitoring data in the fracturing process. This model considered shale heterogeneity and the effect of hydraulic fracturing on the ground stress and pore pressure at the fault interface. At the same time, the fault activation coefficient is established to judge the fault activation degree and fault activation length according to the Mohr-Coulomb criterion. A method to determine the fault activation length is given. Secondly, in the sensitivity analysis of fracture parameters, fault dip angle, fault rock properties, and the variation of fault activation length under different conditions are quantified. Finally, a case study of a shale gas horizontal well with casing deformation was carried out. The effect of varying fracture parameters on fault activation length was analyzed.
2. Physical Model and Fault Activation Determination Method
2.1. The Physical Model
Figure 1 shows the casing deformation pattern of well L1-3 at a depth of 5008 m in L block, Sichuan, China, which was formed using multiarm diameter logging data. The maximum deformation of the casing was 31.52 mm, and the deformation form of the casing was shear deformation. Casing shear deformation was usually caused by fault slip. This characteristic was common in casing deformation caused by fault slip in block Changning and block Weirong.
Figure 2 illustrates well trajectory and fault location relationships on the L1 platform. In Figure 2, the colored bands represent faults, and the black bands represent fracture zones. The position of the hole trajectory for well L1-3 is also indicated. We extracted the relationship between the trajectory and faults in L1-3 in Figure 2 to form a simplified geometric model. The geometric model is shown in Figure 3(a). Meanwhile, the hydraulic fracturing microseismic signal at the casing deformation position is extracted, as shown in Figure 3(b). Microseismic signals are often used to indicate the range of influence of fracturing fluids. In Figure 3(b), the microseismic signals are distributed on the right side of the fault. Therefore, the cause of casing deformation in this section is fault activation caused by hydraulic fracturing. It should also be noted that the microseismic signal did not intersect well L1-4. So, the effect of well L1-4 on well L1-3 could be ignored.
In this paper, the geological and fracturing parameters of well L1-3 in L block at 5008 m are taken as the basis to establish the physical model, as shown in Figure 4. This model is a coupled seepage-stress model for hydraulic fracturing of shale gas horizontal wells with faults. In the model, the shale gas formation is represented as the blue part, which is 2000 m long and 1000 m wide. The fault lies in the middle of the formation and is shown in orange. The horizontal line in the middle of the model is the location of the horizontal well for shale gas. The horizontal well traverses the fault, and the angle between the horizontal well and the fault is . The colored box to the right of the fault is the fracture stage. Each fracture stage is 50 m in length [32]. The horizontal well fracturing direction is from right to left. The distance between the purple fracture stage and the center point of the fault interface is .
2.2. Seepage and Stress Coupling Mathematical Model
In hydraulic fracturing, the pore pressure change in the reservoir causes formation deformation. Due to the impermeable nature of the fault, the stress state of the formation near the fault is different from that of the homogeneous layer. Therefore, the model should consider the effect of hydraulic fracturing on the formation pore pressure and the effect of pore pressure variation on ground stress. We generally adopt the effective stress method to calculate the seepage and stress coupling. The specific equation is as follows: (a)Effective stress equation
Based on the principle of effective stress [33, 34], the total stress is composed of effective stress , pore water pressure , and pore gas pressure . The expression is as follows: where is a function of saturation and is the identity matrix. When the material is completely saturated, . When the material is not saturated, .
During hydraulic fracturing, the water pressure generated by fracturing is much higher than the initial pore pressure of the formation. Therefore, we assume that the analytical process of shale gas does not occur in hydraulic fracturing. Only the variation law of formation pore pressure is considered. Only the variation of formation pore water pressure is considered. At this point, the formation is considered completely saturated. So, we can simplify the above formula as follows: (b)Formation stress balance equation
In this paper, the stress balance equation of strata and fault is expressed by the virtual work principle commonly used in finite element method [35]. where is the virtual displacement field (m); is the virtual strain, ; is the surface force per unit area (MPa); is the volume force per unit volume (excluding fluid weight)(MPa); is the porosity; is the volume (m3); is the area boundary; is the density of the fluid (kg/m3); and is the acceleration due to gravity (m/s2). (c)Equation of seepage continuity
The continuity equation should be satisfied for the fluid in each model element. That is, the fluid flow rate of a finite element equals the increased rate of fluid volume. The specific formula is as follows [36]: where is the average inflow velocity from boundary S (m/s), is the external normal direction vector of boundary S, and is the fluid reference density (kg/m3). (d)Relationship between formation permeability and permeability coefficient
Darcy’s law is used to describe the formation permeability. In ABAQUS, the software uses the permeability coefficient to represent the material’s permeability [36]. The relationship between the permeability coefficient and permeability is as follows: where is the permeability coefficient (dimensionless), is the permeability of a fully saturated medium (D), is a velocity coefficient that depends on the porosity of the material, is the flow rate of the fluid (m/s), and is the relationship between permeability and saturation of the wetting liquid; when saturation is equal to 1, . In this paper, we set the formation as fully saturated.
2.3. Fault Activation Determination Method
In the current research, scholars commonly used the Mohr-Coulomb criterion to determine the activation of faults [37, 38]. In some previous studies, only pore pressure changes in the formation were considered. As shown in Figure 5(a), the reaction size in the Mohr circle is unchanged. Only the Mohr circle is shifted to the left from the right. When the Mohr circle formed by subtracting the pore pressure from the ground stress is tangent to or intersects the Mohr envelope of the fault, it could be judged that the fault was activated. However, the impermeable characteristics of the fault interface rapidly increase the pore pressure at the fault in hydraulic fracturing. This leads to a change of ground stress in the fracturing area of the formation and near the fault. Zhu et al. [39] calculated the variation of ground stress during the development of unconventional reservoirs due to the reduction of formation pore pressure. Similarly, the ground stress of the shale formation changed due to the increase of the formation pore pressure during the fracturing process. When the fracturing fluid enters the fault, the ground stress state on the fault interface also changes with the pore pressure change. Then, the change of the Mohr circle appears as shown in Figure 5(b). The change of maximum and minimum ground stress is not equal, which causes the radius of the Mohr circle to increase or decrease. This change directly affects the judgment of fault activation length. Therefore, a new fault activation coefficient is defined in this paper to determine the activation length of the fault. This coefficient considers the effects of ground stress and formation pore pressure changes on fault activation during shale hydraulic fracturing.
(a)
(b)
According to Mohr-Coulomb’s criterion, we set as the distance between the center of Mohr-coulomb’s circle and the envelope and as the radius of Mohr-Coulomb’s circle, as shown in Figure 6. We define the fault activation coefficient as to describe the fault activation state, and is the ratio of and , as shown in where is the internal friction angle of fault rock (°), is the maximum in situ stress at the fault interface (MPa), is the minimum in situ stress at the fault interface (MPa), and is the pore pressure at the fault interface (MPa).
Thus, appeared in three forms, as shown in Figure 7. When (Figure 7(a)), the Mohr circle is located on the right side of the envelope, and the fault is in an inactive state. When (Figure 7(b)), the Mohr circle is tangent to the envelope, and the fault is in the critical activation state. When (Figure 7(c)), the Mohr circle is partially or entirely to the left of the envelope, and the fault is activated.
(a)
(b)
(c)
3. Mesh Model and Parameter Setting
3.1. Mesh Model
According to the physical model in Figure 4 and the actual situation of well L1-3, we established a numerical model and meshed it in ABAQUS. Figure 8 is a schematic diagram of the numerical model obtained after meshing. The size of the model is . The horizontal well of shale gas is set at the center line of the model, 500 m away from the upper and lower boundaries. The orange part is the fault that intersects the shale gas horizontal well in the model’s center. The angle between the fault and the horizontal well is . In the validation of the model, this angle was set as the angle between the well trajectory of well L1-3 and fault 1#, which was 60°. In the parameter sensitivity analysis, the angle between the fault and the horizontal well was set to 90° for computational convenience without specific instructions. The size of the formation model is much larger than the flow size of the fracturing fluid. This is to eliminate the influence of the model boundary on the variation of pore pressure and ground stress in the fractured part. The fracturing stage is positioned in the horizontal well to the right of the fault. The length of the fracturing stage is set to 50 m according to the actual working conditions. The fracturing sequence is from right to left, with the first stage on the far left and the third stage closest to the fault. The third stage is 25 m away from the fault unless otherwise specified. The mesh size has been optimized before the calculation. The current mesh size met the calculation accuracy requirements and ensured the calculation efficiency.
The mesh elements in the model are quadrilateral elements with fluid-structure interaction. The mesh elements are CPE4P pore pressure-stress coupling elements. The mesh element is a 2D 4-node planar stress element. At the same time, this kind of element is coupled with the calculation of pore pressure. The integration point is at the center point of the element. The coupled solution of pore pressure and stress in the element is introduced in Section 3.4. To ensure the accuracy of the calculation, a large mesh density is used in the fault and near the shale gas horizontal wells. Under the premise of ensuring calculation accuracy, the mesh density is gradually reduced in the vicinity of faults and horizontal wells on both sides. This way, the total number of model meshes can be reduced as much as possible, and the computation efficiency can be improved.
3.2. Boundary Conditions and Load Application
The displacement boundary conditions of the numerical model are defined by imposing displacement constraints on the model. The four edges of the model were set as fixed boundaries, and no displacement or rotation occurred in the three directions.
Faults are generally considered impermeable [40]. There is minimal seepage from the left boundary of the fault into the shale formation on the left. The contact between the formation and the fault is a binding contact mode. This contact mode ensures that the left boundary of the fault only conducts stress without fluid diffusion.
The initial pore pressure of the formation is applied to the four sides of the model, respectively. This means that the pore pressure at the formation boundary is constant. At the same time, the initial pore pressure in the formation is consistent with the initial pore pressure at the boundary before fracturing begins. Through the trial calculation, we found that the increase in pore pressure had almost no effect on the pore pressure and ground stress at the fault interface when fracturing 150 m before the fault.
Fracturing operations are typically performed from stage 1 to stage 3. At the start of fracturing, the pressure in the wellbore is the fracturing pressure. When the stage is fractured, the well is no longer subjected to fracturing pressure. Therefore, the bottom-hole pressure of the fracturing stage is applied through a pore pressure load in this model, which is used to simulate the fracturing process of horizontal wells of shale gas. The pore pressure load is applied in the order of stage 1 to stage 3. After the completion of one stage, the pore pressure load in that stage is deactivated, and the pore pressure boundary load in the next stage is activated.
3.3. Parameter Settings
Ground stresses and initial pore pressure were determined based on the exploration wells in block L, where well L1-3 were located. We determine the initial pore pressure and ground stress in a shale formation based on measurements taken during the drilling history. The initial pore pressure is 30 MPa, the minimum horizontal ground stress is 52 MPa, the maximum horizontal ground stress is 62 MPa, and the vertical ground stress is 65 MPa.
The formation and fault rocks meet the Mohr-Coulomb criterion. The rock mechanical parameters were derived from the field-measured values. In the article, the mechanical properties of the fault were set to be the same as those of the formation without cohesive force [29]. This is because the fault contains fault mud inside, which is highly hydrated when encountering fluids such as fracturing fluid, resulting in damage to the fault when subjected to slight shear force. Specific parameters are shown in Table 1.
The shale permeability during multistage fracturing could be divided into matrix permeability and hydraulic fracture permeability. In this paper, these two components were converted uniformly using equivalent permeability. The formation’s permeability was set to the fracturing fluid’s influence range. Section 2.1 mentions that microseismic signals were commonly used in the field to calibrate fracturing fluids’ influence range. Therefore, the permeability coefficient of the model is determined according to the microseismic signal of well L1-3. This is shown in Section 4.1. The permeability of shale was orthotropic, and the permeability in the direction is about 1.45 times that in the direction. The permeability of the fault was set to match that of the formation without being specified otherwise.
At the time of calculation, the distance from the nearest stage to the fault was 25 m without special instructions. According to Wei’s [32] paper, the length of each fracturing stage is 50 m. The bottom hole fracturing pressure was set as 100 MPa according to the pump pressure and well depth in the fracturing design.
3.4. Solving Method
In ABAQUS, a sequential iterative coupling calculation method is adopted for seepage-stress coupling analysis [41]. The seepage-stress coupling model is divided into two parts. The seepage model and the mechanical model are established and solved, respectively. The calculation process is as follows: firstly, the seepage model is solved. The seepage model is solved in time order. The change in pore pressure is solved at each time step; secondly, the results of the seepage model are transferred to the mechanical model to solve the stress balance equation. The deformation and stress changes of the inland layer are calculated at this time step; finally, the updated pore permeability parameters are transferred to the seepage model to calculate the next time step. In ABAQUS, the time step can be automatically adjusted according to the calculation results to improve the model’s computational efficiency. If there is no convergence within this time step, the software reduces the time step to 1/4 of the original time step and recalculates.
4. Result and Parameter Sensitivity Analysis
4.1. Model Validation
4.1.1. Verification of Calculation Results
We used microseismic signals to verify the model’s accuracy. The microseismic signal generally represents the influence of the fracturing fluid. In previous studies, the profile of microseismic signals was used to model and modify the entire microseismic influence area [42]. This paper assigns the equivalent formation permeability to simulate the entire fracturing process. The calculation results are shown in Figure 9. The computed pore pressure cloud image was used to compare and verify the results with the microseismic signals. Figure 9(a) shows the complete picture of the pore pressure cloud map of the model. Figure 9(b) shows a magnified view near the fault. Figure 9(c) is the microseismic signal diagram of fracturing. The range of pore pressure in the model and the microseismic signals were calculated separately. In the pore pressure cloud image, the area above the initial pore pressure in the model was selected as the area of hydraulic fracturing influence. The influence range of fracturing fluid in the direction is 470 m and in the direction is 324 m. In the microseismic signal diagram, the influence range of fracturing fluid in the direction is 480 m, and the influence range of fracturing fluid in the direction is 330 m. The calculation error is below 5%. The fracturing fluid influence range calculated by the model is consistent with the microseismic signal. Moreover, the influence range is also concentrated on the right side of the fault. The calculation results of this model are close to the engineering practice and meet the engineering requirements.
Furthermore, the fault activation length calculated by seepage-stress coupling and the fault activation length considering only pore pressure changes are analyzed in Figure 10. In Figure 10, the -axis is the distance from the fault plane, and the point A in the mesh model in Figure 8 was set as the point 0 m on the -axis. The -axis on the left shows the calculation results of the fault activation coefficient of the seepage-stress coupling method. The seepage-stress coupling method was calculated using the method in Section 2.3 to obtain the fault activation state at each point along the fault surface. The result is shown as a black line in Figure 10.
The -axis on the right shows the value of fault activation, considering only pore pressure changes. This method is used by Chen et al.’s [43] paper. In this method, the increase of pore pressure causes the overall left shift of Mohr’s circle (Figure 5(a)). In this method, we used the difference between the distance from the circle’s center to the envelope and the radius of the Mohr circle to determine the fault activation state. When the calculation result is less than 0, it indicates that the fault is activated. The blue line represents it in Figure 10.
The calculated curves in Figure 10 have a U-shaped distribution, with the tip of the U-shape representing a greater activation of the fault here and a lesser activation of the fault on both sides. The two different background color regions represent the fault activation lengths under the two methods. The green background represents the calculation results of Chen’s method, and the orange represents the calculation results of this paper’s method. According to the calculation results, the fault activation length obtained by the seepage and stress coupling method is 428 m. When only the seepage method is considered, the activation length of the fault is only 298 m. The difference between the two calculations was 43%. According to Zhao et al.’s [44] article, the activated fault length of the example well is 337.5 m, the fault slip distance can be 23~43 mm, and the casing deformation can be between 10 mm and 20 mm. The casing deformation in well L1-3 was 31.25 mm, indicating that the fault activation length should be greater than 298 mm. According to the formula in Chen et al.’s [12] article, it is calculated that when the fault activation length is 428 m, the fault slip is between 33 mm and 62.3 mm. According to Zhang et al.’s [45] calculation, the casing deformation can reach more than 30 mm when the fault slip distance exceeds 50 mm. This calculation is consistent with the multiarm diameter measurements. The results show that the coupled model is more suitable for engineering practice.
According to the calculation results, since a large amount of liquid is injected into the hydraulic fracturing process of shale gas in a short time, the change of ground stress at the fault interface significantly influences the fault activation. Therefore, the influence of ground stress variation at the fault interface on fault activation under hydraulic fracturing conditions must be addressed.
4.1.2. Analysis of Calculation Results
Generally, the angle between the wellbore trajectory and the fault is more than 45°. Furthermore, the angle of 90° is the majority. The angle between the fault and horizontal well trajectory of 26 wells in block L and block W was calculated. The ant body and the well trajectory of some wells are shown in Figure 11. Figure 12 shows the statistical results. The 85° to 90° between the fault and horizontal well is 65.39%. Most of the faults are perpendicular to the horizontal well trajectory. Therefore, in the subsequent calculation, we adopt the condition that the angle between the fault and the horizontal trajectory is 90° to conduct the rule analysis for convenience. The advantage is that the graph is more symmetrical, which can better show the influence law of parameters on fault activation.
(a) Ant body diagram of block W
(b) Ant body diagram of block L
Based on the method described in Section 2.3, the pore pressure, formation stress, and fault activation factor were calculated for different fracturing stages. The calculation results are shown in Figures 13–17.
(a)
(b)
(c)
Figure 13 shows the fracturing fluid diffusion cloud diagram before and after the 1~3 fracturing stage. It can be seen that almost no fluid enters the fault in the first two fracturing stages, and the pore pressure in the fault does not change much. After the third fracturing stage, the fracturing fluid enters the fault. This is also illustrated in Figure 14. Figure 14 represents the pore pressure at the fault interface. There was a slight increase in the formation pore pressure during the first two stages. This is due to the fault’s impermeability and the compression of the fault caused by the increase of pore pressure in the right formation, leading to a slight rise of pore pressure at the fault interface. After the third stage is fractured, the flow of the fault fluid into the fault increases the pore pressure at the fault interface.
However, we find that the variation of ground stress and pore pressure is different in Figure 15. This is due to the formation being squeezed and displaced by the fracturing fluid as the fracturing fluid enters the formation. Due to the existence of faults, the displacement of the formation leads to a change in the ground stress at the fault interface. In Figure 15(a), the minimum horizontal ground stress increases first and then decreases at the fault center. After the second fracturing stage, the minimum horizontal ground stress in the fault center is more significant than in the first stage. The minimum horizontal ground stress in the center of the fault is greater than that in the first fracturing stage. However, the variation is not very large, and the variation range is about 1 MPa. After the third fracturing stage, the fracturing fluid enters the fault, and the impermeable nature of the fault allows the fracturing fluid to flow only inside the fault. This causes a decrease in the normal stress on the fault interface, which leads to a decrease in the minimum horizontal ground stress on the fault interface. This rule can also be seen in the minimum ground stress cloud image in Figure 16. Similar rules are also reflected in the maximum horizontal and vertical crustal stresses, as shown in Figures 15(b) and 15(c). The impact of the first two stages was also minimal. Both of these stresses decreased significantly after the third stage was fractured. This is the same result as previous studies have shown. When the fracturing fluid enters the fault, the pore pressure on the formation increases, and the normal stress decreases, resulting in lower friction and fault slip. The difference is that the maximum horizontal and vertical crustal stresses at the fault interface gradually decrease with the progress of fracturing instead of increasing first and then decreasing. This is because the fault is not squeezed in the and directions, but the stress decreases due to increased pore pressure.
Finally, the activation coefficients of each stage were calculated, and the results are shown in Figure 17. In Figure 17, the -axis is the position of the fault interface. The zero point is the position of the fault midpoint, the negative number is above the fault center in the model, and the positive number is below the fault center in the model. The -axis shows the fault activation coefficient. It can be seen that the activation coefficient at the fault interface is greater than 1 during the first and second fracturing stages, and the fault is in a state of no activation at this time. After fracturing the third stage, the activation coefficient on the fault surface starts to fall below 1. This indicates fault activation after the third fracturing stage.
We noticed that the activation factor was higher in the center of the fault than on the sides during the first two stages. This phenomenon is related to the variation of pore pressure and ground stress. This indicates that the fracturing fluid did not enter the fault during the first two fracturing stages. However, the fracturing fluid began to squeeze the formation at the fault interface. This increases the normal stress on the fault and makes it less likely to activate it. This further indicates that the influence of ground stress at the fault interface on fault activation during hydraulic fracturing cannot be ignored. It also indicates that fault activation occurs only when the fracturing fluid enters the fault and causes a significant increase in pore pressure inside the fault. This is consistent with Ji et al.’s [46] point. Fracture fluid entry into the fault plane is the most critical factor for fault activation.
4.2. Sensitivity Analysis
4.2.1. Influence of Fracturing Parameters
Two critical parameters in the fracturing process are fracturing pump pressure and fracturing fluid injection volume. The fracturing pump pressure response in the model is the fracturing pressure at the bottom of the well. The variation of bottom hole pressure can change the influence of fracturing fluid. The fluid injection volume is represented by the fracturing duration in the model. The variation of fracturing time affects the volume of fracturing fluid injected into the formation, that is, the injection volume.
In this section, we set the bottom hole fracturing pressure as 80 MPa, 90 MPa, and 100 MPa, and the fracturing duration as 2 h, 3 h, and 4 h, respectively. The influence of downhole fracturing pressure and fracturing duration on fault activation length was analyzed.
Figure 18 shows the influence of bottom hole pressure on the fault activation length. As the bottom hole pressure decreases, the fault activation length decreases. When the bottom hole pressure is 100 MPa, the fault activation length is 340 m. Moreover, when the bottom hole pressure is reduced to 80 MPa, the fault activation length is 200 m. The activation length of the fault decreased by 41.2%. According to Zhao and Zhang’s paper, we can obtain that the fault slip decreases by about 18 mm. According to Zhang’s paper, a 10 mm reduction in fault slip can effectively reduce casing deformation. However, the bottom hole pressure depends on the well’s vertical depth and the formation’s mechanical properties in the actual fracturing process, such as the fracture pressure. All these factors bring engineering difficulties to controlling the bottom hole pressure.
Figure 19 shows the influence of fracturing duration on fault activation length. As the fracturing duration decreases, the activation length of the fault decreases significantly. When the fracturing duration is 4 h, the fault activation length is 340 m. When the fracturing duration is reduced to 2 hours, the fault activation length is only 100 m. The reduction can reach 70.6%. The fault slip distance decreases to 9.1 mm with a fault activation length of 100 mm, and the casing deformation at 9.1 mm is only 0.86 mm. The reduction in fracturing duration can effectively reduce the fault activation length and fault slip distance. The variation in fracturing duration also reduces the amount of fluid injected. The control of fault slip could start from this aspect.
4.2.2. Influence of Fault Rock Mechanical Parameters
The rock mechanics of the fault determines the type of fracture parameters we use for different formations. However, the properties of fault rocks vary significantly with different burial depths and stratigraphic origins. Shale reservoirs are currently buried between 1500 m and 4000 m, which can significantly affect the mechanical parameters of fault rocks. Therefore, we set the different elastic modulus, Poisson’s ratio, and internal friction angle of fault rocks to analyze the variation rule of fault activation length.
The fault elastic modulus was set as 20 GPa, 30 GPa, and 40 GPa, and the calculation results are shown in Figure 20. The activation length of the fault increases with the increase of the elastic modulus of the fault. When the elastic modulus is 20 GPa, the fault activation length is 340 m. When the elastic modulus is 40 GPa, the fault activation length increases to 540 m. With the development of shale gas, the shale reservoir is buried deeper, and the elastic modulus of the formation also increases. Shale depth increases, resulting in longer fault activation lengths and more casing deformation problems. This is also consistent with current field knowledge. The buried depth of the Changning block is about 2500 m, and the casing deformation ratio is 23% (49/211). The buried depth of the Weirong block is 2800~3000 m, and the occurrence probability of casing deformation is 41.8% (23/55) [47]. The buried depth of the Luzhou block is 3500~3800 m. As of June 2022, the probability of casing deformation in fractured wells was 59.02% (13/22). Even 23 unfractured wells had casing deformation. With the increased buried depth, the casing deformation becomes increasingly serious. Therefore, more attention should be paid to the fault activation caused by hydraulic fracturing for shale reservoirs with deep burial depth or high elastic modulus.
The Poisson ratio of the fault was set as 0.20, 0.25, and 0.30, and the calculation results are shown in Figure 21. With the increase of fault Poisson’s ratio, the activation length of the fault increases and the activation coefficient at the fault center increases. When Poisson’s ratio is 0.2, the fault activation length is 400 m, while when Poisson’s ratio increases to 0.3, the fault activation length is 480 m. The increase is only 20%. The fault Poisson’s ratio is not the main factor affecting the fault activation length.
One of the essential properties affecting rock’s shear failure is the internal friction angle. We set the internal friction angles of the fault as 20°, 30°, and 40°, respectively. Moreover, the calculation results are shown in Figure 22. With the fault rock’s internal friction angle increasing, the activation length of the fault decreases gradually. This conclusion is also by the Mohr-Coulomb criterion. That is, the larger the internal friction angle of the rock, the greater the shear strength of the rock, and the more difficult the fault activation. The fracture displacement can be appropriately increased to maximize the fracturing effect of the fault with a larger internal friction angle.
4.2.3. The Influence of the Distance from the Injection Section to the Fault
The distance of the fracturing stage from the fault was set to be 25 m, 50 m, and 75 m. After the third stage fractured, the pore pressure and ground stress were selected to calculate the fault activation coefficient. The influence of the fracturing stage distance from the fault on the fault activation length is analyzed. The calculation results are shown in Figure 23. When the distance between the fracturing stage and the fault is greater than 50 m, the activation coefficient of the entire fault section is greater than 1, indicating that the fault is inactive. Moreover, when the distance equals 50 m, the activation coefficient at the fault center is close to the critical activation state. Therefore, during the fracturing process, the distance between the fracturing stage and the fault should be controlled according to the formation characteristics to avoid extensive fault activation.
At the same time, the fault activation coefficient is abnormal in the middle of the fault. The calculation results of ground stress in two directions at the fault section are shown in Figure 24. It can be found that the minimum horizontal ground stress and the maximum horizontal ground stress are reversed in the center of the fault. This is due to the impermeable nature of the fault interface, where the fracturing fluid can only percolate within the fault. The transfer of pressure triggers a change in the ground stress at the fault interface. The reduction of the minimum horizontal ground stress is less than the reduction of the maximum horizontal ground stress.
4.2.4. Influence of Fault Dip Angle
In Section 4.1, we discussed the angle between the fault and horizontal well trajectory and found that most of the angles were above 75°, but a few were 60° or even smaller. Therefore, the included angle between the horizontal well trajectory and the fault is set as 45°, 60°, and 90° in this section. The influence of different fault dip conditions on the fault activation length was calculated. The results are shown in Figure 25.
The calculation results show that the smaller the angle between the fault and the horizontal well, the longer the activation length of the fault. This is because when the angle between the fault and the horizontal well is less than 90°, the fracturing fluid can communicate with the fault not only at the point where the fault meets the horizontal well but also below the fault. This leads to greater variation in pore pressure within the fault, as shown in Figure 26. However, when the angle is 60°, the activation length of the fault is close to that of the angle is 45°. This indicates that the increment of fault activation length gradually decreases and eventually becomes stable under constant fracture parameters. However, the degree of fault activation becomes larger. The activation coefficient becomes smaller, leading to greater fault slip distance.
5. Case Study
This paper selects an example of Wei’s [32] paper for analysis. The well is the N2-5 well located in block C, Sichuan Basin, China. It is a shale gas horizontal well with a well depth of 4450 m and a vertical depth of 2500 m. The well depth of point A in the horizontal well is 2950 m. The horizontal section is 1500 m long. Casing shear deformation occurred in well N2-5 at 3802 m, which was judged to be caused by fault activation and slip. The positional relationship between the wellbore trajectory and the ant body of the fault in well N2-5 is shown in Figure 27. In this section, the casing deformation point at 3802 m is selected for the case study.
It can be seen from Figure 27 that the well trajectory at casing deformation point 1 is in the eccentric position of the fault. The fault and the well trajectory are perpendicular. Therefore, the physical model of well N2-5 is established according to Figure 27. The fault length is 500 m long based on the scale. The horizontal well trajectory is eccentric with the fault section, and the angle between the fault and the horizontal well is 90°. It is shown in Figure 28.
In the numerical model, a mesh model was established to reduce the influence of the boundary effect at the model’s boundary on the stress and pore pressure at the fault, as shown in Figure 29. Shale formations were added around the fault. The effect of the boundary on pore pressure and ground stress variation at the fault boundary is reduced. The size of the numerical model is . Calculation parameter settings are shown in Table 1.
According to the numerical simulation results, the data were extracted and processed by the method in Section 2.3. The result of the fault activation coefficient of well N2-5 was obtained. The calculation results are shown in Figure 30. It can be seen from the figure that the fault activation coefficient before fracturing is 2.1, which is larger than the fault activation coefficient 1. The fault is overall stable. After fracturing, the activation coefficient of the fault decreases due to the change in pore pressure and ground stress. The fault activation length is about 312 m. According to the formula in Chen et al.’s [14] paper, the fault slip distance caused by this fault activation length is 31 mm. A 31 mm fault slip distance is about 8 mm with a 139.7 mm casing deformation. The amount of deformation is sufficient to prevent the plug from passing through. Subsequent fracturing operations are affected by casing deformation.
To investigate which method can effectively reduce the activation length of the fault, we propose three improvements of fracturing parameters for the fracturing stage closest to the fault. (1) reducing the bottom hole fracturing pressure from 100 MPa to 80 MPa; (2) reducing the duration of fracturing by half, that is, reducing the amount of liquid injected by half; (3) increasing the distance between the stage and the fault. The distance from the fracture stage to the fault is doubled. The calculation results are shown in Figure 31. The results show that all three methods are effective. All three methods can reduce the fault activation length. Thus, the fault slip and casing deformation can be reduced. The activation lengths obtained by the three methods are 198 m, 150 m, and 100 m, respectively. Reducing fracturing duration (or fluid injection) is most effective in reducing the fault activation length. The activation length of the fault was reduced by 68.25%.
The field engineer’s primary concern is that the fault slip distance created does not interfere with subsequent fracturing operations. The fault slip distance was calculated under three fault activation lengths. The fault slip distances were 18.05 mm, 13.6 mm, and 9.1 mm, respectively. According to the fault slip distance, the casing deformation was calculated using the numerical model in Meng et al. [48]. The casing deformation under the three conditions was calculated as 4.5 mm, 1.2 mm, and 0.86 mm. The outer diameter of the bridge plug is generally 2 to 4 mm smaller than the outer diameter of the casing. Therefore, when the casing deformation is less than 2 mm, the running of the bridge plug is not affected [32]. So increasing the distance between the fracture stage and the fault by more than 50 m or reducing the fracturing duration by 2 h can reduce the casing deformation to less than 2 mm.
Therefore, when faults are encountered in shale gas hydraulic fracturing, it is recommended to increase the distance between the fracturing stage and the fault or reduce the fracturing duration (i.e., reduce the amount of fluid injected). This way, the activation length of the fault can be better reduced, and the amount of fault slip can be reduced. This results in less casing deformation and less stage loss. Meanwhile, reducing fracture pressure to reduce the fault slip is not recommended. This makes it difficult to control the fracturing pressure. Further, reducing fracture pressure reduces fracturing effectiveness and the production of horizontal wells. This is the opposite of the goal of fracturing.
In the subsequent fracturing of block C, hydraulic fracturing near the fault was performed in two ways: (1)Reducing fracturing duration reduces the amount of fluid entering the formation, which ensures that the activation length in the stage near the fault is reduced. The probability of casing deformation was effectively reduced, with the casing deformation ratio decreasing from 44% in 2018 to 16.7% in 2020(2)The fracture stage across the fault is fractured with the next stage. This reduced the casing loss rate. The lost stage rate of horizontal wells decreased from 4.58% in 2018 to 1.72% in 2020
The two methods in this paper effectively reduced the proportion of casing deformation and the stage loss rate. These optimization methods of fracturing parameters have been proven to improve the shale gas wells’ construction efficiency and have good popularization value.
6. Conclusion
This paper studied the effect of hydraulic fracturing on fault activation in shale gas horizontal wells. A numerical model of shale formation was developed to simulate the fault activation state during hydraulic fracturing. At the same time, the fault activation coefficient was established to evaluate the activation state of each point on the fault surface according to the Mohr-Coulomb criterion. Moreover, the fault activation length under fracturing conditions was calculated. The influence of different engineering and geological conditions on fault activation length was analyzed. The following conclusions were drawn: (1)A seepage-stress-coupled hydraulic fracturing model for shale gas horizontal wells was established. The variation law of pore pressure and ground stress at the fault interface was analyzed. The fault activation coefficient was established to evaluate the fault activation state and calculate the fault activation length during shale gas hydraulic fracturing. According to the calculation results, the rapid increase of formation pore pressure in hydraulic fracturing resulted in the change of ground stress at the fault interface. The influence of ground stress variation on fault activation could not be ignored(2)A sensitivity analysis was conducted. The effects of fracturing parameters, mechanical parameters of the rock, fracturing stage, distance from the fault, and fault dip on the activation length of the fault were analyzed. The numerical results showed that fracture pressure, fracturing duration, fault elastic modulus, and fault Poisson’s ratio positively correlated with the fault activation length. The distance between the fracturing stage and fault, friction angle, and fault dip angle negatively correlated with the fault activation length(3)A case study was conducted using the data from C block, Sichuan Basin, China. The casing deformation data was used for backstepping verification. Three methods to reduce the fault activation length were proposed. The results showed that reducing the fracturing duration and increasing the distance between the fracture stage and the fault were the most effective methods for reducing the fault activation length. These two methods have good application and popularization values
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant No. U19B6003, Grant No. 52004013, and Grant No. 52204018). This paper is also supported by the talent introduction project of the Karamay Campus of China University of Petroleum, Beijing (no. XQZX20220019).