The crack initiation and propagation process is a dynamic coupling of hydraulic fracturing fluid flow and rock deformation. Based on Biot theory (Charlez, 1991; Chen et al., 1995), the linear elastic strain of reservoir rock and pore pressure variation can be represented by Equation 1 and Equation 2: φ c t ∂ p 0 ∂ t − α ∂ ∇ · u ∂ t = ∇ · k μ · ∇ p 0 (1) λ + G ∂ ∇ · u ∂ x i + G ∇ 2 u i − α ∂ p 0 ∂ x i = 0 (2) where φ is porosity, c t is overall compression coefficient, t is time, α is the Biot coefficient, u is the displacement vector, k is the permeability matrix, μ is the fluid viscosity, p 0 is the pore pressure, ∇ is the vector operator, λ and G are the Lamé constant, x i is the position vector.

The stress balance equation, derived from the principle of virtual work, can be expressed by Equation 3: ∫ V σ ′ + P 0 I : δ ε υ d V = ∫ S γ δ V υ d S + ∫ V f δ V d V (3) where δ ε is the virtual strain rate; γ is the surface stress per unit area; f is the body force per unit volume; I is the identity matrix; V is the accumulated deformation volume; δ V is the velocity field of the reservoir matrix; S is the loading area; and v is the strain rate of the reservoir matrix.

Assuming the fracturing fluid is an incompressible Newtonian fluid, the flow control equation of Newtonian fluid in hydraulic fractures can be obtained through the fracture pressure loss equation and the law of conservation of mass (Charlez, 1991), as shown in Equation 4: ∂ ∂ x w 3 12 μ ∂ p 0 ∂ x = ∂ w ∂ t + V 1 x , t (4) where w is the fracture width; μ is the viscosity of the fracturing fluid; and V 1 x , t is the fracturing fluid loss rate.

The leakage of fracturing fluid from the fracture into the rock matrix can be characterized by (Carter, 1957), as shown in Equation 5: V L = − k y μ ∂ p ∂ y (5) where k y is the permeability of the rock along the normal direction of the fracture surface, ∂ p ∂ y is the pressure gradient along the normal direction of the racture surface.

The interaction between contact surfaces includes normal and tangential action. As shown in Equations 6, 7, when the clearance h between surfaces satisfies h < 0, the surfaces are considered to be in contact with each other with contact pressure p > 0 (Figure 2). Conversely, for h ≥ 0, the contact surfaces are separated from each other with p = 0. p > 0 h < 0 close (6) p = 0 h ≥ 0 open (7)

The slip type includes smooth straight slip and smooth curve slip. The node X₁ in the linear slip segment is the node on the slave surface. The node X is the closest point to X₁ on the master surface, and n and t are the normal and tangent unit vectors at node X, respectively. If the model boundary is a smooth straight line, the shear displacement segment trajectory can be determined using two nodes, X₂ and X₃ (Figure 3a). If the model boundary is a smooth curve, the shear displacement trajectory can be determined by establishing a quadratic equation using three nodes, X₂, X₃ and X₄. If the contact occurs at the (convex) vertex of two segments, only a single node will enter the equations (Figure 3b).

An interpolation function N i g is defined to get the position of node X_i , as shown in Equation 8: n h = N i g X i (8)

The tangent line t of node X could be expressed by Equation 9: t def = d X d s = d N i d g X i / d X d g (9)

Linearize Equation 10 yields the contact slip equation: δ n h + n δ h = d N i d g X i δ g + N i δ X i = t δ s + N i δ X i (10) where δ s and δ h are the slip distance in the tangential and normal directions of the slip surface, respectively, expressed by Equations 11, 12: δ h = N i n · δ X i (11) δ s = − N i t · δ X i (12)

The stress stiffness of point X is obtained by derivation of functions for h and s , as shown in Equations 13, 14: d δ h = − δ X i n d N i d s N j t + t N i d N j d s n + t N i ρ n N j t · d X j (13) d δ s = δ X i t d N i d s N j t − n N i d N j d s n − n N i ρ n N j t · d X j (14) where ρ n is the curvature of the slip segment.

The simulation results were compared with experimental results by Wang et al. to validate the applicability of the proposed method (Wang et al., 2020). A cylindrical sandstone sample with a length of 1 m and a diameter of 0.5 m was prepared, a sawcut plane was located at the center of the cylinder with a dip angle α = 30°. The elastic modulus of sandstone was 15 GPa, Poisson’s ratio was 0.26, and the friction coefficient was 0.7.

The main test procedures were as follows. First, the sample was loaded until the required confining pressure of 35 MPa was reached, while the pore pressure was kept constant at 5 MPa. Subsequently, the shear strength (τ_ss) was obtained by axial loading at an axial displacement rate of 1 μm/s. Next, the axial load was reduced with a displacement rate of 0.05 μm/s until the shear stress reached 0.92 × τ_ss. The axial pressure was kept constant, and water was pumped into the sample to trigger fault slip. Notably, the confining pressure remained constant at 35 MPa throughout the test.

The top of the sample was sealed to create an undrained boundary condition and pore pressure was applied through Pump A at the bottom of the sample. In order to study the correlation between fault slip and fluid pressure, in the test the fluid pressure increased from an initial value of 5 MPa at a rate of 2 MPa/min, with each increment lasting 2 min. Subsequently, the pore pressure was kept constant for 8 min. Changes in pore pressure and slip character--ristics on the sawcut plane were recorded by acoustic emission (AE) sensors. The experimental setup is shown in Figure 4a.

A 3D hydro-mechanical coupling model was established referring to the experimental settings (Figure 4b). The model geometry, material parameters and stress loading strategy were consistent with those of the test. The friction coefficient was set to 0.7, the void ratio was 0.29. The specific parameters were listed in Table 1 (Wang et al., 2020).

The fluid pressures and the slip distance on the sawcut plane obtained by numerical model is compared with the experimental results (blue bar in Figure 5). In general, the results obtained by numerical model is basically consistent with the experimental results, the regression coefficient reaches 0.97, indicating the proposed method could be applied in further research.

A numerical model of a vertical well was constructed based on the lithological parameters and fracturing strategy for a representative well in the Weiyuan gas field, Sichuan Basin. The buried depth of mudstone-sandstone interface was 3,200 m, the sandstone layer and mudstone layer were in angular unconformity, with a stratum dip angle of 40°. The friction coefficient was 0.7, and the void ratios of mudstone and sandstone were 0.1 and 0.2, respectively. The lithological parameters are shown in Table 2. The model size was set as 10 × 10 × 20 m, and the diameter of the well was 0.18 m (Figure 6a), and the mesh was finer near the wellbore to accurately capture the response characteristics of the interface (Figure 6b). The stress regime around the well was strike-slip type. The model employed fixed-displacement boundary conditions with an initial stress field defining in-situ stresses: maximum horizontal stress σ _H = 75 MPa, minimum horizontal stress σ _h = 55 MPa, and vertical stress σ _v = 60 MPa. A constant pore pressure boundary of 40 MPa was applied.

According to field reports, the drilling fluid density was 1,500 kg/m³ and the bottom hole pressure during the drilling stage was 60 MPa. During the multistage fracturing, the bottom hole pressure linearly increased from 60 MPa to 80 MPa (breaking pressure) over 1,200 s, then decreased to 70 MPa during 1,200 s–2,400 s, and finally stabilized at 70 MPa until 7,200 s (Figure 6c).

The pore pressure distribution and principal stress along the formation interface during drilling and fracturing process are shown in Figure 7. The pore pressure at the wellbore rose to 60 MPa after drilling, while the normal stresses σ ₁ and σ ₃ were reduced by 14.4 MPa and 6.97 MPa, respectively (Figure 7c). During fracturing, the pore pressure at the wellbore increased sharply and diffused outward. The permeability difference between the two sides of the formation interface and the dip angle hindered the lateral diffusion of the injected fluid, thus, the injected fluid to concentrate near the interface on the left side of the wellbore (Figure 7a), resulting in a larger pressure change on this side (Figure 7b). The pore pressure diffusion ranges on the left and right sides of the wellbore were 4.8 m and 3.5 m, respectively. The overpressure generated during the fracturing process led to a decrease in the contact stress at the interface, and the range of contact stress reduction increased with time. As shown in Figure 7c, the impact of pore pressure on the maximum principal stress was more significant than the minimum principal stress, thus reducing the stress difference around the wellbore. The contact stress near the wellbore reached 0 MPa, indicating that the tensile stress generated at the wellbore may open the interface.

The stress variation, interface opening and slip distance at the formation interface during drilling and fracturing process were shown in Figure 8. In the drilling process, the changes of contact stress and shear stress were only observed near the wellbore, and no slip occurred at the formation interface. After fracturing, the stress change on the left side of wellbore was more significant due to the uneven fluid pressure distribution (Figures 8a,b). The tension stress led to partial opening of the formation interface, and the opening near the wellbore first increased and then decreased with time. At 1,200 s the opening was 0.015 mm. Between 2,400 s and 7,200 s, after a linear reduction the pore pressure remained constant. Fluid leakage at the formation interface reduced the opening to 0.013 mm, but enlarged its range (Figure 8c). The slip occurred mainly during the fracturing stage, and the slip distance increases with the fracturing time. Influenced by the dip angle and the difference in rock permeability, the injected fluid concentrated on the left side of the wellbore, resulting in a more significant reduction in principal stress. The interface opening and contact pressure decrease led to significant slip at the formation interface (Figure 8d). Fluid aggregation near the wellbore progressively expanded the stress-affected zone, displacing the peak slip location along the formation interface farther from the injection point over time.

As shown in Figure 9, both contact pressure and shear stress at the wellbore decreased significantly by 2,400 s, but the contact pressure decreased more sharply than the shear stress. For this reason, slip occurred mainly in the first 2,400 s. After 2,400 s, as the contact stress and shear stress remained basically unchanged, the interface slip remained constant.

To study the mechanical response characteristics of different interface dip angles during fluid injection, the formation interface dip angle α was set to 0°–60° while other parameters kept constant. According to Equations 15,16 the contact stress and shear stress first increase and then decrease as α increases (Figure 10a,b). In addition, as shown in Figure 11b, when α = 10°, the direction of shear stress changed due to fluid injection, and it reached a maximum when α = 40°–50°. σ = 1 2 σ 1 + σ 3 + 1 2 σ 1 − σ 3 cos ⁡ 2 α (15) τ = 1 2 σ 1 − σ 3 sin ⁡ 2 α (16)

Figure 10c presents the formation interface opening after fracturing. The opening first increases and then decreases as α increases. The slip distance at the formation interface after the fracturing follows a similar trend, reaching a high value for α between 40° and 50° (Figure 10d). The maximum slip distance occurred at the wellbore when the dip angle of the formation interface was between 0° and 20°. However, when the dip angle was between 30° and 60°, the injected fluid concentrated in the sandstone on the left side of the wellbore, resulting in larger induced stress, and the maximum slip occurred about 1 m to the left of the wellbore.

The permeability coefficient of sandstone formation increased from 2.0 × 10⁻¹¹ m/s to 1.0 × 10⁻¹⁰ m/s with other factors kept constant. The contact pressure of formation interface decreased as the increase of permeability coefficient (Figure 11a). As shown in Figure 11b, the shear stress decreased with the increase of sandstone permeability coefficient. Meanwhile, both the magnitude and range of the interface opening increased, extending from 4.5 m to 5.1 m (Figure 11c). The slip distance at the formation interface increased with permeability coefficient. The maximum slip distance was 0.45 mm, and the slip distance at the wellbore was 0.35 mm (Figure 11d). In addition, the maximum slip position at the formation interface occurred about 1 m to the left of the wellbore, which did not change with increasing permeability coefficient, indicating that the formation dip angle was the key control factor for the location of the maximum interface slip.

The stress difference of the formation was adjusted by changing the maximum horizontal in-situ stress. The variation range of the horizontal maximum principal stress was set to 55 MPa–75 MPa, and the variation range of the horizontal stress difference was 0 MPa ∼ 20 MPa. As shown in Figure 12a, the contact pressure near the wellbore (within 1 m from the wellbore) decreased with the increase of the stress difference, indicating that the induced stress increased with stress difference, and shear stress similarly increased with stress difference. When the stress difference was 10 MPa–20 MPa, the model was in strike-slip state, and when the stress difference was reduced to between 0 MPa and 5 MPa, the model was in normal stress state, the direction of shear stress near the wellbore was reversed after fracturing (Figure 12b), leading to a change in slip direction in Figure 12d. As shown in Figure 12c, the opening of the formation interface increased as the stress difference rosed from 0 MPa to 20 MPa, and the opening range on the left side of the wellbore increased from 3.8 m to 4.9 m. Notably, when the stress difference was 10 MPa, fluid injection changed the local stress state near the wellbore from strike-slip to normal, causing slip direction vary along the interface and resulting in a smaller slip distance. Therefore, the slip distance around the wellbore first decreased and then increased with the increase of the stress difference, and the maximum slip distance was 0.36 mm when the stress difference rose to 20 MPa.