XSite, a hydraulic fracturing simulation software developed based on discrete lattice theory and synthetic rock mass technology, has gained widespread academic recognition for its precision and adaptability in simulating fracture propagation (Zhao et al., 2021; Bakhshi et al., 2019). The 3D discrete lattice method integrates synthetic rock mass (SRM) technology and discrete lattice theory, utilizing a simplified bonded particle model (BPM) to accurately and intuitively describe the initiation, propagation, and evolution of hydraulic fractures in rock masses (Zhao et al., 2021), as shown in Figure 1.

Within this theoretical framework, the rock medium is discretized into a system of interacting particles, where mechanical responses between particles are simulated via a spring-mass system. The tensile and shear failure of springs correspond to the tensile and shear failure mechanisms of rock materials. Fluid units are located at the centers of springs, forming a fluid transport network that enables fluid flow between discrete units. Additionally, to accurately characterize pre-existing discontinuities in geological structures, XSite adopts the smooth joint model (SJM) to simulate the mechanical behaviors of faults, joints, and other geological weak planes, including complex responses such as sliding, opening, and closing.

The core advantage of XSite lies in its ability to directly simulate highly nonlinear mechanical behaviors during rock deformation, such as fracturing, sliding, and joint opening/closing. The dynamic evolution of these behaviors is solved using a central difference algorithm at each node. The translational motion of nodes is governed by Equations 1, 2: u ˙ i t + Δ t 2 = u ˙ i t − Δ t 2 + ∑ F i t Δ t m (1) u i t + Δ t = u i t + u ˙ i t + Δ t / 2 Δ t (2) Where u ˙ i t and u i t represent the velocity and position of component i (i = 1,3) at time t, ∑ F i denotes the sum of all forces acting on the node of mass m, and Δt is the time step. This explicit integration scheme ensures computational efficiency and numerical stability.

Rock deformation and fracture are simulated through the relative displacements between lattice nodes, with mechanical responses following the incremental relationship calculated by Equations 3, 4: F N ← F N + u ˙ N k N Δ t (3) F i S ← F i S + u ˙ S k S Δ t (4) Where k N , k S are the normal and shear stiffness of the spring; F N and F i S are the normal and shear forces; u ˙ N and u ˙ S are the normal and shear relative displacements.

XSite employs a coupled hydro-mechanical model to simulate fluid flow within fractures and the matrix. Fracture flow, based on Darcy’s law, is solved through a network of fluid nodes and pipes constructed dynamically during simulation. Matrix flow represents permeability, storage, and leakage in porous media, characterized by pore pressure fields stored in spring elements. Fluid exchange between these mechanisms is driven by pressure gradients, achieving fully coupled simulations of fracture-matrix systems.

This methodology excels in capturing complex physical phenomena during hydraulic fracturing, including fracture initiation, propagation, branching, and interactions with natural fractures, providing a reliable theoretical foundation and numerical tool for designing and optimizing hydraulic fracturing in unconventional reservoirs.

A numerical simulation method for hydraulic fracture propagation in continental shale was established using the 3D discrete lattice algorithm. By constructing a 3D mechanical model that incorporates reservoir-interlayer structural characteristics, we systematically studied the control mechanisms of in-situ stress differences and rock mechanical parameters on hydraulic fracture height extension.

A 3D discrete lattice method is used to construct a hydraulic fracturing model for interlayer interstibed shale, which can accurately simulate the complex mechanical behavior of fracture propagation in heterogeneous layers. The model dimensions are 25 m × 25 m × 5 m, comprising a central reservoir layer and upper/lower interlayers. A horizontal wellbore is positioned at the model center, aligned parallel to the minimum horizontal stress σ h direction to simulate field conditions. Boundary conditions apply triaxial in-situ stresses: minimum horizontal stress σ h , maximum horizontal stress σ H , and vertical stress σ v , reflecting realistic subsurface stress environments, as shown in Figure 2.

The key parameters required for model calculation are shown in Table 1. In order to study the influence of mechanical properties of interlayer on fracture height extension, the elastic modulus of upper interlayer was set in a range of 6000–30000 MPa and Poisson’s ratio was set to 0.1–0.3 to reflect the mechanical properties of different lithologic interlayer. The elastic modulus of the reservoir and the lower interlayer is fixed at 12,000 MPa and Poisson’s ratio is 0.2. In the model, the same strength parameters were used for both reservoir and interlayer: uniaxial compressive strength of 65 MPa, tensile strength of 4 MPa and fracture toughness of 4 MPa·m^0.5 to highlight the influence of elastic parameters and stress differences. In the setting of ground stress conditions, the foundation stress state of reservoir and interlayer is the same: the minimum horizontal principal stress is 10 MPa, the maximum horizontal principal stress is 20 MPa, and the vertical ground stress is 20 MPa. On this basis, by adjusting the minimum horizontal principal stress value of the interlayer, the effect of stress difference between reservoir and interlayer on fracture height extension is studied. In terms of fracturing construction parameters, the injection displacement was set at 0.05 m3/s and the fracturing fluid viscosity was set at 10 mPa·s to simulate the actual hydraulic fracturing conditions.

Based on the 3D discrete lattice theory mentioned above and the numerical simulation framework of XSite software, the accuracy of the model is verified in this section, which provides a reliable theoretical basis for the subsequent research on hydraulic fracturing simulation of laminated shale. By comparing the numerical simulation results with the classical theoretical analytical solutions, the accuracy and reliability of the model in predicting hydraulic fracture propagation are evaluated.

Under ideal conditions (no fracturing fluid loss, homogeneous isotropic medium, constant injection rate), the expansion of circular hydraulic fracture meets the GeerSMA -de Klerk model, and the theoretical analytical solution of its radius can be expressed as Equation 5: R t = 0.52 E Q c 3 t 4 μ 1 − v 2 1 / 9 (5) Where, R(t) is the radius of hydraulic fracture, m; E is the elastic modulus, Pa; Q c is Fluid injection rate, m/s; t is the injection time, s; v is Poisson’s ratio, dimensionless; μ is the fluid viscosity, Pa·s. The analytical solution is based on the theory of linear elastic fracture mechanics, assuming that the fracture remains circular during expansion and the fluid flow satisfies laminar flow conditions.

In order to verify the calculation accuracy of XSite software based on 3D discrete lattice method, a standard verification model is constructed. The model adopts the same discrete lattice theoretical framework as Section 2.1, and the rock medium is discretized into a particle system, and the interaction between particles is simulated by a spring-particle system. The size of the verified model is 25 m × 25 m × 5 m, and the uniform mesh is used to ensure sufficient calculation accuracy.

The material parameter Settings of the model are consistent with the basic parameters in Table 1. The initial fracture is placed in the center of the model to ensure effective injection of fracturing fluid. The simulation process uses the explicit time integration algorithm described in Section 2.1, and the time step is automatically adjusted according to the computational stability conditions to ensure the convergence and accuracy of the numerical solution. The coupling of fluid flow and rock deformation is realized by iterative solution in each time step to ensure the coupling accuracy.

Figure 3 shows the comparison between XSite numerical simulation results and theoretical analytical solutions under different injection times. The abscissa represents the injection time (s) and the ordinate represents the fracture radius (m). It can be clearly seen from the figure that in the whole simulation process, the numerical solution is in good agreement with the theoretical solution, which verifies the accuracy of XSite software in simulating hydraulic fracture propagation.

It is worth noting that the numerical simulation results were able to capture some small fluctuations that were not reflected in the theoretical solution, mainly due to the fact that the discrete lattice method can more accurately simulate the stress concentration and energy release process near the crack tip. In addition, the simulation results also show that the distribution of crack opening width is basically consistent with the theoretical prediction, which further confirms the reliability of the model.

In this study, the influence mechanism of elastic modulus of different interlayer on the height growth of hydraulic fracturing was quantitatively analyzed. The unified initial geostress field conditions are set up in the experiment: the maximum horizontal principal stress is 20 MPa, the vertical geostress is 20 MPa, and the minimum horizontal principal stress is 10 MPa. The elastic modulus of the reservoir and the lower interlayer is maintained at a constant value of 12,000 MPa, while the elastic modulus of the upper interlayer is controlled as a variable, which is set at 6000, 9000, 12,000, 20,000 and 30,000 MPa respectively. The fracture longitudinal propagation characteristics obtained through numerical simulation are shown in Figure 4.

The results show that the elastic modulus of the upper interlayer is positively correlated with the longitudinal propagation of cracks. When the overlying elastic modulus is low (6000 MPa), vertical fracture propagation is significantly inhibited. With the gradual increase of elastic modulus, the height expansion of fractures is gradually enhanced, and the volume of fractures in overlying strata is also increased.

Specifically, when the elastic modulus of the upper interlayer is 6 GPa, the fracture height expansion is strictly limited and the permeability volume is very small, which is mainly concentrated in the reservoir area. When the elastic modulus of the upper interlayer increases to 9 GPa, the fracture height propagation increases slightly, but it still shows obvious limitation. When the elastic modulus of the upper interlayer reaches 12 GPa (the same as that of the reservoir), the fracture height increases significantly, and the volume in the overlying layer expands obviously, showing uniform expansion. When the elastic modulus of the upper interlayer further increases to 20 GPa and 30 GPa, the fracture height extension continues to increase, and gradually shows a trend of preferential expansion of the overlying strata, and reaches the maximum expansion degree when the elastic modulus is 30 GPa.

This phenomenon can be explained by the principle of elastic mechanics: strata with higher elastic modulus have smaller strain when subjected to the same stress, so they have weaker resistance to fracture propagation. However, the formation with lower elastic modulus is more prone to deformation, which is a strong obstacle to fracture propagation. This is consistent with previous studies on the control effect of rock mechanical properties on fracture propagation (Bakhshi et al., 2019), and this influence relationship is further quantified.

The comprehensive analysis shows that under the same in-situ stress conditions and injection parameters, the ratio of elastic modulus between upper interlayer and reservoir is the key parameter to control vertical fracture propagation: when the ratio <1, vertical fracture propagation is inhibited; When the ratio is greater than 1, vertical fracture propagation is promoted. When the ratio is 1, the fracture spreads uniformly.

Stress barriers are a key factor in controlling the height propagation of hydraulic fractures, especially in the stiff layer under compression. In order to quantify the effect of the difference of ground stress on the fracture penetration behavior, a systematic numerical simulation experiment scheme is designed. The maximum horizontal principal stress is 20MPa, the vertical ground stress is 20 MPa, and the elastic modulus of the reservoir and the upper and lower interlayers are 12000 MPa. The horizontal minimum principal stress of the reservoir is fixed at 10 MPa, while the horizontal minimum principal stress of the upper interlayer is set at 8, 10, 12, and 15 MPa, respectively, to study the longitudinal propagation characteristics of the fractures under different stress difference conditions.

Figure 5 shows the simulation results of longitudinal crack propagation under different stress difference conditions. It is found that the difference of horizontal minimum principal stress between reservoir and interlayer has a significant control effect on fracture height propagation. When the minimum horizontal geostress difference between reservoir and interlayer is less than 5MPa, the stress barrier has a weak sealing effect on the fracture, and the fracture can penetrate the upper interlayer relatively easily, forming a larger volume fracture network in the overlying formation. With the increase of the horizontal ground stress difference between reservoir and interlayer, the fracture height growth is gradually inhibited, and the fracture volume in the overlying strata is obviously reduced, and finally the fracture height growth stops.

Through quantitative analysis of the crack propagation morphology under different stress difference conditions, the control mechanism of stress difference on the crack penetration behavior can be obtained: When the crack spreads from the low stress area to the high stress area, additional stress resistance needs to be overcome, and this process requires more energy consumption. With the increase of stress difference between reservoir and interlayer, the pressure required to overcome this energy barrier also increases. When a certain threshold is reached, the fracture cannot continue to expand upward.

Through the parameter sensitivity analysis of the system, the critical conditions for fracture penetration are quantified, as shown in Figure 6: When the minimum horizontal ground stress difference between layers is greater than 7 MPa, vertical fracture propagation will be effectively inhibited. When the ratio of elastic modulus between layers is less than 0.5, the fracture height propagation is also significantly limited.

The impact of the elastic modulus difference between the mudstone interlayer and the sandstone reservoir on fracture height was investigated through numerical simulation. The results show that the difference of elastic modulus between mudstone interlayer and sandstone reservoir has significant influence on the hydraulic fracture height. The experimental results show that the mudstone interlayer has obvious limiting effect on fracture height. With the decrease of the elastic modulus of mudstone interlayer, the width of hydraulic fracture increases, and the fracture height decreases under the same fluid injection condition. In addition, the toughness characteristics of mudstone make it difficult to form a high stress concentration area at the tip of hydraulic fracture, which increases the difficulty of fracture initiation and extension.

Figure 7 quantitatively shows the mechanism of sandstone-mudstone modulus difference affecting fracture height. When the elastic modulus of the mudstone interlayer is 11 GPa lower than that of the sandstone reservoir (Figure 7a), the hydraulic fracture propagates vertically in the sandstone reservoir, but after entering the mudstone interlayer, the propagation is significantly inhibited, showing the characteristics of increasing width and limiting height. When the difference of elastic modulus decreases to 6 GPa (Figure 7b), the crack height increases correspondingly, indicating that the difference of elastic modulus is negatively correlated with the crack height.

In addition, the difference of minimum horizontal principal stress between mudstone interlayer and sandstone reservoir also has a key influence on fracture height. Figure 8 shows this regulatory effect: When the minimum horizontal principal stress of the mudstone interlayer is 4 MPa lower than that of the sandstone reservoir (Figure 8a), the fracture propagation resistance in the mudstone interlayer is small and the fracture height is larger. However, when the stress difference is reversed to a higher value of 4 MPa (Figure 8b), the crack height decreases significantly. This reveals the positive correlation between the difference of ground stress and fracture height, and provides an important theoretical basis for reservoir reconstruction design.

To explore the effects of fracturing fluid parameters on fracture penetration, numerical simulations were conducted based on the model parameters described in the preceding sections to investigate the influence of different fracturing fluid parameters on fracture penetration. The results show that the viscosity of fracturing fluid and the flow rate are the key parameters that affect the penetration ability of fractures. The high viscosity fracturing fluid significantly enhances the ability to penetrate the interlayer interface, and the injection strategy of “high viscosity leading and low viscosity following” can effectively achieve the dual goals of main fracture expansion and weak interlayer opening.

Figure 9 shows the mechanism by which fracturing fluid viscosity affects fracture penetration. When the viscosity is 5 mPa·s, the low-viscosity fracturing fluid promotes the expansion of the main fracture, but the opening effect of the weak surface between layers is not good. When the viscosity increased to 20 mPa·s, the ability of high-viscosity fracturing fluid to penetrate the interlayer interface was significantly improved, and both the main fracture and the weak interlayer were well opened. By using the compound injection scheme of 20 mPa·s followed by 5 mPa·s, the cooperative optimization of the expansion of the main seam and the opening of the weak surface between the layers was realized.

At the same time, fracturing fluid injection rate plays a significant role in regulating fracture height extension distance. The experimental results show that with the injection displacement increasing from 3.0 mL/min to 7.0 mL/min, the longitudinal penetration height of hydraulic fracture shows an obvious increasing trend. Figure 10 quantitatively describes the variation law of crack height expansion distance under different displacement: When the displacement is 3.0 mL/min, crack longitudinal expansion is significantly inhibited. When the displacement was increased to 5.0 mL/min, the expansion of seam height was improved. When the displacement reaches 7.0 mL/min, the fracture longitudinal extension effect is most significant. It is worth noting that with the increase of injection displacement, the fracture extension pressure gradually increases, especially when the artificial fracture extends to the upper and lower mudstone layers, the extension pressure increases sharply. When the injection displacement increases from 3.0 mL/min to 7.0 mL/min, the fracture extension pressure in the interlayer increases by about 10 MPa.

Under the condition of displacement of 6 m³/min, the longitudinal penetration of seam height is significantly inhibited. Figure 11 shows the fracture propagation patterns at different times under this displacement: fractures start in the sandstone layer and begin to spread to the mudstone layer in 30 s; The fracture propagation height in mudstone layer is limited at 110 s. At 180 s, the fracture propagation in the mudstone layer is completely blocked, and the fracture height propagation stops. This is mainly due to the fact that when hydraulic fractures spread from sandstone to mudstone reservoir, the lower elastic modulus of mudstone leads to the increase of fracture width, reduces the stress concentration of crack tip, and hinders the further extension of fractures.

Figure 12 Systematically compares the fracture propagation patterns under different discharge rates: when the discharge rate is 6 m³/min, fracture propagation is blocked in mudstone layer; When the displacement increased to 8 m³/min and 9 m³/min, the seam height expansion gradually increased but did not fully penetrate the layer. When the displacement reaches 10 m³/min, the fracture can completely penetrate the layer, but the fracture length decreases significantly, indicating that there is a critical relationship between the displacement and the penetration ability.

In summary, high viscosity fracturing fluid and high injection rate are conducive to hydraulic fracture penetration of the interlayer interface. The compound injection strategy of “high viscosity leading and low viscosity following” can take into account the dual needs of main seam expansion and weak surface opening between layers. With the increase of injection displacement, the fracture extension pressure and longitudinal expansion height both increase, but the tradeoff between displacement and fracture length should be paid attention to.